https://sangaku0418.hatenablog.com/entry/2024/02/28/003155

福島県安達郡東和町木幡山（現二本松市）隠津島神社明治17年(1884)

街角の数学 Street Wasan　～落書き帳「○△□」～ 15．一寸
http://streetwasan.web.fc2.com/math15.5.22.html
キーワード：円2個，半円，扇
#Julia #SymPy #算額 #和算 #数学

半径が 1 寸の扇の中に，甲半円と乙円が入っている。乙円の直径が最大になるとき，甲円の直径はいかほどか。

扇の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ R - r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ R - r_1 + r_2)\)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive
eq1 = x2^2 + r2^2 - (r1 + r2)^2
eq2 = x2^2 + (R - r1 + r2)^2 - (R - r2)^2
res = solve([eq1, eq2], (r1, x2))

2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
( (-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))
( (R^(3/2)*sqrt(R - 8*r2)/2 + sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(R^(3/2)*sqrt(R - 8*r2)/2 + sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))

2 組の解が求まるが，基本的にはどちらも適解である。
どちらの解も，\(R - 8r_2\) の平方根を求める操作が必要なので，乙円の半径 \(r_2\) は扇の半径 \(R\) の 1/8 よりも大きくはなれない。すなわち \(r_2\) の最大値は \(R/8\) であり，このときの \(r_1\) が求める答えである。

扇の半径が 1 寸ならば，\(r_2\) の最大値は 1/8 = 0.125 であり，このとき \(r_1\) は 0.5 寸（直径は 1 寸）である。

res[1][1](R => 1, r2 => 0.125) |> println

0.500000000000000

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

function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = 1
r2 = R/8
(r1, x2) = ( (-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))
@printf("扇の半径 = %g; 乙円の直径 = %g; 甲円の直径 = %g\n", R, 2r2, 2r1)
r3 = R - r1
y = R - r1
x = sqrt(R^2 - y^2)
θ = atand(y/ x)
plot()
circle(0, 0, R, beginangle=θ, endangle=180-θ, n=500)
circle(0, 0, r3, :magenta, beginangle=θ, endangle=180-θ, n=500)
circle(0, y, r1, :gray80)
circle(0, y, r1, :blue, beginangle=0, endangle=180)
circle2(x2, y + r2, r2, :green)
segment(-x, y, x, y)
segment(0, 0, x, y)
segment(0, 0, -x, y)
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, :center, :bottom, delta=delta)
point(0, R - r1, "甲半円:r1,(0,R-r1)", :blue, :center, :bottom, delta=delta)
point(x2, R - r1 + r2, "乙円:r2,(x2,R-r1+r2) ", :green, :right, :vcenter)
end
end;