87 岩手県室根村折壁字天王前 弥栄神社 明治16年(1833)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円14個，外円
#Julia #SymPy #算額 #和算 #数学
外円の中に，甲円，乙円，丙円，丁円を容れる。甲円の直径は十二寸五分，「甲乙円径五分三」（注1）のとき，丙円の直径はいかほどか。
注1：これは，「乙円の直径は甲円の直径の 3/5」であろう。ほかに解釈のしようがない。乙円の直径は七寸五分である。しかしこの時点で，算額のような図は描けないことが明らかである。下図の A が直径が 12.5 寸の甲円である。
甲円がどの程度互いに交差するのかは変化しうるが。乙円 B は直径が 7.5 寸である。乙円が外円に接するように描いているが，甲円と乙円は交わってしまう。甲円の交わりが大きくなると，丙円との交わりも大きくなってしまう。逆に甲円の交わりが小さくなり，交わりが 0 になってもまだ甲円と乙円は交わったままである。甲円の交わりがなくなって（中心が離れて）は，算額の図のようにはならない。
なぜこのようなことになるかといえば，それは乙円の大きさを固定してしまったからである。実際以下のように連立方程式を立ててみると，条件が 1 個多すぎる。それは乙円を既知とするからである。そこで，以下では乙円の直径は未知として解を求めることにする。
外円の半径と中心座標を \(R, (0, 0)\)
甲円の半径と中心座標を \(r_1, (0, R - r_1)\)
乙円の半径と中心座標を \(r_2, (x_2, r_2)\)
丙円の半径と中心座標を \(r_3, (x_{31}, 0), (x_{32}, y_3)\)
丁円の半径と中心座標を \(r_4, (0, 0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive,
x2::positive, r3::positive, x31::positive,
x32::positive, y3::positive, r4::positive
eq1 = x2^2 + r2^2 - (R - r2)^2
eq2 = x32^2 + y3^2 - (R - r3)^2
eq3 = x2^2 + (R - r1 - r2)^2 - (r1 + r2)^2
eq4 = x31^2 + (R - r1)^2 - (r1 + r3)^2
eq5 = x32^2 + (y3 - R + r1)^2 - (r1 + r3)^2
eq6 = (x2 - x32)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq7 = (x2 - x31)^2 + r2^2 - (r2 + r3)^2
eq8 = 2r1 - r4 - R;
## res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (R, r2, x2, r3, x31, x32, y3, r4))
function H(u)
#(R, x2, r3, x31, x32, y3, r4) = u
(R, r2, x2, r3, x31, x32, y3, r4) = u
return [
r2^2 + x2^2 - (R - r2)^2, # eq1
x32^2 + y3^2 - (R - r3)^2, # eq2
x2^2 - (r1 + r2)^2 + (R - r1 - r2)^2, # eq3
x31^2 + (R - r1)^2 - (r1 + r3)^2, # eq4
x32^2 - (r1 + r3)^2 + (-R + r1 + y3)^2, # eq5
(-r2 + y3)^2 - (r2 + r3)^2 + (x2 - x32)^2, # eq6
r2^2 - (r2 + r3)^2 + (x2 - x31)^2, # eq7
-R + 2*r1 - r4, # eq8
]
end;
r1 = 12.5/2
iniv = BigFloat[1.9, 0.45, 1.38, 0.24, 0.85, 1.22, 1.12, 0.099] .*r1
res = nls(H, ini=iniv)
([11.881752126800865, 2.815876063400433, 8.617479375809692, 1.5081750129387437, 5.335976715996894, 7.63184686468951, 7.02609532892722, 0.6182478731991338], true)
甲円の直径が 12.5 のとき，乙円の直径は 5.63175212680087（注2），丙円の直径は 3.01635002587749（注3）である。
その他のパラメータは以下のとおりである。
r1 = 6.25; r2 = 2.81588; r3 = 1.50818; r4 = 0.618248; R = 11.8818
0.45054017014406933
注2：「問」では乙円の直径は甲円の直径の 3/5=0.6 と言っていたが，実際には約 0.45 である。先に示した図でも，乙円の直径が言っているより小さければ，甲円に外接し外円に内接するという条件を満たすことは自明であった。
注3：「答」では，「乙円径五十分十九」といっている。五十分の十九寸でないことは図を見ても明らかである。
では，別の解釈は丙円の直径はなにかの 19/50 だろうが，該当するようなものは見当たらない。「問」に続いて素直に読めば，「丙円は甲円の 19/50」と言いたいのだろうが，物差しで測ってみると「丙円径/甲円径≒1/6」のようだ（計算結果では「丙円径/甲円径≒1/4」)。
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 12.5/2
(R, r2, x2, r3, x31, x32, y3, r4) = res[1]
@printf("甲円の直径が %g のとき，乙円の直径は %.15g，丙円の直径は %.15g である。\n", 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; r3 = %g; r4 = %g; R = %g\n", r1, r2, r3, r4, R)
println(r2/r1)
plot()
circle(0, 0, R)
circle22(0, R - r1, r1, :blue)
circle4(x2, r2, r2, :green)
circle4(x32, y3, r3, :magenta)
circle2(x31, 0, r3, :magenta)
circle(0, 0, r4, :orange)
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, " R", :red, :left, :bottom, delta=delta/2)
point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :blue, :center, delta=-delta)
point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
point(x31, 0, "丙円:r3\n(x31,r3)", :magenta, :center, :vcenter)
point(x32, y3, "丙円:r3\n(x32,r3)", :magenta, :center, :vcenter)
point(0, 0, "丁円:r4,(0,0)", :black, :center, delta=-2delta)
end
end;
以下のアイコンをクリックして応援してください
