一〇八 埼玉県加須市騎西町 玉敷神社 大正4年(1915)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個
#Julia #SymPy #算額 #和算 #数学
外円の中に甲円，乙円を 2 個ずつ入れ，甲円，乙円に外接する丙円を 2 個描く。甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき，丙円の直径はいかほどか。
後述するが，上の図は外円，甲円，乙円の直径が 20, 6, 3 のときのもので，丙円の直径は 12.3435 である。
その他のパラメータは，
\(R = 10;\ r_1 = 3;\ r_2 = 1.5;\ y_1 = -6.32456;\ y_2 = 8.3666\)
\(r_3 = 6.17176;\ y_3 = 2.28132\)
外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (r_1,\ y_1)\)
乙円の半径と中心座標を \(r_2,\ (r_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (r_3,\ y_3)\)
とおき，以下の連立方程式を解く。
まとめて解くには SymPy の能力が足りない。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive,
r1::positive, y1::negative,
r2::positive, y2::positive,
r3::positive, y3::positive
eq1 = r1^2 + y1^2 - (R - r1)^2 |> expand
eq2 = r2^2 + y2^2 - (R - r2)^2 |> expand
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;
外円，甲円，乙円の直径が与えられるので，eq1, eq2 はそれぞれだけで未知数 \(y_1, y_2\) を求めることができる。
(ans_y1, ans_y2) = solve([eq1, eq2], (y1, y2))[1]
(-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))
\(y_1, y_2\) が既知として，eq3, eq4 の連立方程式を解く。
(y1, y2) = (-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;
(ans_r3, ans_y3) = solve([eq3, eq4], (r3, y3))[1] # 1 of 2
( (-R*r1 + R*r2 + (sqrt(R)*sqrt(R - 2*r1) + sqrt(R)*sqrt(R - 2*r2))*(-sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2)))/(2*(r1 - r2)), -sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2))
\(R,\ r_1,\ r_2\) を与えれば，それぞれのパラメータの数値解が得られる。
甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき，以下のようになる。
ans_y1(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y2(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_r3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
-7.21110255092798
11.4017542509914
7.73565831477412
4.58897582448691
「術」では，「1385.4 を 89.57 で割れば丙円の直径が得られる」とあるが，「甲円，乙円，外円の直径がそれぞれ 9 寸，3 寸，26 寸のとき」の解であるし 1385.4/89.57 = 15.467232332254104 で，「答」の「丙円径十五寸四分二厘」とも小数点以下第2位で不一致である。
丙円の直径の正確な値は上に示すように 2*7.73565831477412 = 15.47131662954824 である。
任意の 甲円，乙円，外円の直径に対しては，上述の数式解に定数を代入すればよい。
\(y_1,\ y_2,\ r_3,\ y_3\) の式は長く，SymPy では簡約化できないが，描画ソースプログラム内に示すように，一時変数を使えばプログラム的には短くなる。
描画関数プログラムのソースを見る
function draw(R, r1, r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
s = √(R - 2r1)
t = √(R - 2r2)
u = -√(2R*r1*r2*(R - r1 - r2 + s*t))
v = √R*(r1*t + r2*s)
w = r1 - r2
(y1, y2, r3, y3) = (
-s√R,
t√R,
(s + t)*(u + v)√R/(2w^2) - R/2,
(u + v)/w
)
@printf("外円，甲円，乙円の直径が %g, %g, %g のとき，丙円の直径は %g\n", 2R, 2r1, 2r2, 2r3)
@printf("R = %g; r1 = %g; r2 = %g; y1 = %g; y2 = %g; r3 = %g; y3 = %g\n", R, r1, r2, y1, y2, r3, y3)
plot()
circle(0, 0, R, :brown)
circle2(r1, y1, r1, :green)
circle2(r2, y2, r2)
circle2(r3, y3, r3, :blue)
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", :green, :left, :bottom, delta=delta/2)
point(r1, y1, "甲円:r1,(r1,y1)", :green, :center, delta=-delta/2)
point(r2, y2, "乙円:r2,(r2,y2)", :black, :center, :bottom, delta=delta/2)
point(r3, y3, "丙円:r3,(r3,y3)", :blue, :center, delta=-delta/2)
end
end;
draw(20/2, 6/2, 3/2, more=false)
以下のアイコンをクリックして応援してください
