https://sangaku0418.hatenablog.com/entry/2025/04/08/175026

宮城県塩竈市森山塩竈神社文化14年(1817) 昭和60年(1985)復元奉納

http://www.wasan.jp/miyagi/siogama4.html
キーワード：円3個，楕円
#Julia #SymPy #算額 #和算 #数学

楕円の中に，甲円，乙円，丙円を容れる。楕円の長径が 13 寸，短径が 5 寸，乙円の直径が 4 寸のとき，丙円の直径はいかほどか。

楕円の長半径，短半径，中心座標を \(a,\ b,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (d_1,\ 0)\)
乙円の半径と中心座標を \(r_2,\ (d_2,\ 0)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
丙円と楕円の接点座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

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

using SymPy
@syms a, b, r2, r1, r3, x3, y3, d1, d2, x0, y0
d1 = sqrt( (a^2 - b^2)*(b^2 - r2^2))/b
d2 = sqrt( (a^2 - b^2)*(b^2 - r1^2))/b
eq1 = a^4*(r1 - r2)^2 - 4b^2*(a^2*b^2 - a^2*r2*r1 - b^4) # 算法助術-公式87
eq2 = (x3 + d1)^2 + y3^2 - (r3 + r2)^2
eq3 = (x3 - d2)^2 + y3^2 - (r3 + r1)^2
eq4 = (x0 - x3)^2 + (y0 - y3)^2 - r3^2
eq5 = x0^2/a^2 + y0^2/b^2 - 1
eq6 = -b^2*x0/(a^2*y0) + (x0 - x3)/(y0 - y3);
# solve([eq1, eq2, eq3, eq4, eq5, eq6], (r1, r3, x3, y3, x0, y0))

ans_r1 = solve(eq1, r1)[2] # 2 of 2
ans_r1 |> println

(2*b*sqrt( (a - b)*(a + b)*(b - r2)*(b + r2)) + r2*(a^2 - 2*b^2))/a^2

ans_r1(a => 13/2, b => 5/2, r2 => 4/2).evalf() |> println

2.47337278106509

function H(u)
(r1, r3, x3, y3, x0, y0) = u
return [
a^4*(r1 - r2)^2 - 4*b^2*(a^2*b^2 - a^2*r1*r2 - b^4), # eq1
y3^2 - (r2 + r3)^2 + (x3 + sqrt( (a^2 - b^2)*(b^2 - r2^2))/b)^2, # eq2
y3^2 - (r1 + r3)^2 + (x3 - sqrt( (a^2 - b^2)*(b^2 - r1^2))/b)^2, # eq3
-r3^2 + (x0 - x3)^2 + (y0 - y3)^2, # eq4
-1 + y0^2/b^2 + x0^2/a^2, # eq5
(x0 - x3)/(y0 - y3) - b^2*x0/(a^2*y0), # eq6
]
end;
a = 13/2
b = 5/2
r2 = 4/2
iniv = BigFloat[2.47337, 0.62, -1.67, 1.8, -1.7, 2.42]
res = nls(H, ini=iniv)

([2.473372781065089, 0.63, -1.6666666666666667, 1.783009316358785, -1.7333333333333334, 2.4094720491334933], true)

長径 = 13, 短径 = 5, 乙円の直径 = 4 のとき，丙円の直径は 1.26 である。

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

function draw(a, b, r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
#r1 = (2*b*sqrt( (a - b)*(a + b)*(b - r2)*(b + r2)) + r2*(a^2 - 2*b^2))/a^2
(r1, r3, x3, y3, x0, y0) = [2.473372781065089, 0.63, -1.6666666666666667, 1.783009316358785, -1.7333333333333334, 2.4094720491334933]
d2 = -sqrt( (a^2 - b^2)*(b^2 - r2^2))/b
d1 = sqrt( (a^2 - b^2)*(b^2 - r1^2))/b
@printf("長径 = %g, 短径 = %g, 乙円の直径 = %g のとき，丙円の直径は %g である。\n", 2a, 2b, 2r2, 2r3)
plot()
ellipse(0, 0, a, b)
circle(d2, 0, r2)
circle(d1, 0, r1, :blue)
circle(x3, y3, r3, :green)
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(x0, y0, "(x0,y0)", :green, :center, :bottom, delta=2delta)
point(d1, 0, "甲円:r1,(d1,0)", :blue, :center, delta=-delta)
point(d2, 0, "乙円:r2,(d2,0)", :red, :center, delta=-delta)
point(x3, y3, "丙円:r3,(x3,y3)", :green, :left, delta=-delta, deltax=8delta)
point(a, 0, " a", :black, :left, :vcenter)
point(0, b, "b", :black, :center, :bottom, delta=delta)
end
end;