https://sangaku0418.hatenablog.com/entry/2024/12/16/120757

32 岩手県一関市山目字館配志和神社嘉永5年(1852)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円17個，外円
#Julia #SymPy #算額 #和算 #数学

外円の中に中円を容れ，甲円，乙円，丙円を同数個ずつ互いに外接するように容れる。甲円と乙円は外円に内接し，丙円は中円と外接する。中円の直径が与えられたとき甲円の個数から外円の直径を求める術を述べよ。

外円の半径と中心座標を \(R,\ (0,\ 0)\)
中円の半径と中心座標を \(r_0,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ R - r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
とおく。
原点と乙円の中心を結ぶ直線は丙円の中心を通る。
甲円の中心，原点，乙円の中心がなす角を \(\theta\) とすると，甲円，乙円，丙円の個数がそれぞれ \(n\) 個のとき，\(\theta = \pi/n\) である。
\(n\) は前もって与えられており，それにより \(\theta,\ s = \sin\theta，c = \cos\theta\) も前もって計算されているとする。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy

@syms n, θ, R, r0, r1, r2, x2, y2, r3, x3, y3
@syms c, s
# θ = PI/n
# c = cos(θ)
# s = sin(θ)
eq1 = x2^2 + (y2 - (R - r1))^2 - (r1 + r2)^2
eq2 = x3^2 + (y3 - (R - r1))^2 - (r1 + r3)^2
eq3 = (r0 + 2r2 + 2r3) - R
eq4 = (r0 + 2r1) - R
eq5 = x2 - (R - r2)*s
eq6 = y2 - (R - r2)*c
eq7 = x3 - (r0 + r3)*s
eq8 = y3 - (r0 + r3)*c;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (R, r1, r2, x2, y2, r3, x3, y3))[3]; # 3 of 4

例えば，\(n = 6\) のとき，R, \(r_1,\ r_2,\ x_2,\ y_2,\ r_3,\ x_3,\ y_3\) はそれぞれ以下のように計算される。

n = 6
θ = π/n
res[1](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[2](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[3](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[4](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[5](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[6](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[7](c => cos(θ), s => sin(θ), r0 => 1) |> println
res[8](c => cos(θ), s => sin(θ), r0 => 1) |> println

2.12423571211134
0.562117856055670
0.333379955279379
0.895427878415980
1.55092657993009
0.228737900776291
0.614368950388145
1.06411823666503

外円の半径は，以下のように長い式で表される。

res[1] |> println

r0*(c^6 + 3*c^4*s^2 + c^4 - 8*c^3 + 3*c^2*s^4 - 2*c^2*s^2 - 4*c^2*sqrt(c^6 - 2*c^5 + 2*c^4*s^2 - 17*c^4 - 4*c^3*s^2 + 36*c^3 + c^2*s^4 - 32*c^2*s^2 + 63*c^2 - 2*c*s^4 + 36*c*s^2 - 162*c - 15*s^4 + 126*s^2 + 81) - 77*c^2 - 8*c*s^2 + 72*c + s^6 - 3*s^4 - 4*s^2*sqrt(c^6 - 2*c^5 + 2*c^4*s^2 - 17*c^4 - 4*c^3*s^2 + 36*c^3 + c^2*s^4 - 32*c^2*s^2 + 63*c^2 - 2*c*s^4 + 36*c*s^2 - 162*c - 15*s^4 + 126*s^2 + 81) - 105*s^2 - 12*sqrt(c^6 - 2*c^5 + 2*c^4*s^2 - 17*c^4 - 4*c^3*s^2 + 36*c^3 + c^2*s^4 - 32*c^2*s^2 + 63*c^2 - 2*c*s^4 + 36*c*s^2 - 162*c - 15*s^4 + 126*s^2 + 81) - 117)/(c^6 + 3*c^4*s^2 - 7*c^4 + 8*c^3 + 3*c^2*s^4 - 10*c^2*s^2 - 13*c^2 + 8*c*s^2 - 72*c + s^6 - 3*s^4 + 15*s^2 - 45)

SymPy では簡約化できないので，SymPy の手助けを借りながら手動で以下のように簡約化できる。

ans_R = r0*(4sqrt( (1 - c)*(c + 3)) + 7 - c*(c + 2))/(c + 1)^2
ans_R(c => cos(θ)) |> println

2.12423571211134*r0

甲，乙，丙円が \(n\) 個ずつの場合，外円の直径は中円の直径の \( (4\sqrt{ (1 - c) (c + 3)} + 7 - c(c + 2))/(c + 1)^2\) 倍である。
ここで，\(\theta = \pi/n,\ s = \sin \theta,\ c = \cos\theta\) である。

中円の直径が 1 のとき，\(n = 3, \dots, 10\) に応じた外円の直径は以下のようになる。

for n = 3:10
θ = π/n
s = sin(θ)
c = cos(θ)
@printf("n = %2d, 外円の直径 = %.15g\n", n, 2(4sqrt( (1 - c)*(c + 3)) + 7 - c*(c + 2))/(c + 1)^2)
end

n = 3, 外円の直径 = 9.81466899744816
n = 4, 外円の直径 = 6.35082454034786
n = 5, 外円の直径 = 4.97417671222864
n = 6, 外円の直径 = 4.24847142422268
n = 7, 外円の直径 = 3.80359616681383
n = 8, 外円の直径 = 3.50404360012199
n = 9, 外円の直径 = 3.28902767325085
n = 10, 外円の直径 = 3.12737748258041

以下は，図を描くためのプログラムである。
\(n,\ r_0\) を与えて全ての描画パラメータを得る関数は以下のように定義できる。

function getparameters(arg_n, arg_r0)
@syms n, θ, R, r0, r1, r2, x2, y2, r3, x3, y3
n = arg_n
# θ = Sym(180)/n # 180 をシンボルにすると長精度演算になり n < 7 でないと解けない
θ = pi/n
eq1 = x2^2 + (y2 - (R - r1))^2 - (r1 + r2)^2
eq2 = x3^2 + (y3 - (R - r1))^2 - (r1 + r3)^2
eq3 = (r0 + 2r2 + 2r3) - R
eq4 = (r0 + 2r1) - R
eq5 = x2 - (R - r2)*sin(θ)
eq6 = y2 - (R - r2)*cos(θ)
eq7 = x3 - (r0 + r3)*sin(θ)
eq8 = y3 - (r0 + r3)*cos(θ);
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (R, r1, r2, x2, y2, r3, x3, y3))[4]
[float(x(r0 => arg_r0).evalf()) for x in res]
end;

2*getparameters(3, 1)[1] |> println
2*getparameters(10, 1)[1] |> println

9.814668997448164
3.127377482580407

getparameters(6, 1) |> println

[2.124235712111339, 0.5621178560556694, 0.333379955279379, 0.8954278784159799, 1.5509265799300853, 0.22873790077629041, 0.6143689503881452, 1.064118236665031]

描画関数プログラムのソースを見る

function draw(n, r0=1, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r2, x2, y2, r3, x3, y3) = getparameters(n, r0)
plot()
circle(0, 0, R, :magenta)
circle(0, 0, r0, :orange)
rotate(0, R - r1, r1, :blue, angle=360/n)
rotate(x2, y2, r2, :green, angle=360/n)
rotate(x3, y3, r3, angle=360/n)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, R - r1, "甲円:r1,(0,R-r1)", :blue, :center, delta=-delta/2)
point(x2, y2, "乙円:r2\n(x2, y2)", :green, :center, delta=-delta/2)
point(x3, y3, "丙円:r3,(x3, y3)", :red, :left, delta=-delta/2, deltax=-delta)
point(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
point(0, r0, "r0", :black, :center, delta=-delta/2)
end
end;