算額（その2110）
百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，扇
#Julia, #SymPy, #算額, #和算
扇面に大円 1 個，等円 3 個を容れる。扇長（要から先端までの長さ）が与えられたとき，大円の直径を求める術を述べよ。
扇長を R，扇の端（図参照）の座標を (x0, sqrt(R^2 - x0^2))
大円の半径と中心座標を r1, (0, R - r1)
等円の半径と中心座標を r2, (0, R - 2r1 - r2), (x2, y2)
残念ながら，答えを得るための「術」は SymPy の能力では得られないようだ（何らかの手立てはあるはず）。
以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms R::positive, x0::positive, r1::positive,
r2::positive, x2::positive, y2::positive;
sinθ = x0/R
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x2^2 + (y2 - R + r1)^2 - (r1 + r2)^2
eq3 = dist2(0, 0, x0, sqrt(R^2 - x0^2), 0, R - r1, r1)
eq4 = dist2(0, 0, x0, sqrt(R^2 - x0^2), 0, R - 2r1 - r2, r2)
eq5 = dist2(0, 0, x0, sqrt(R^2 - x0^2), x2, y2, r2);
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(r1, r2, x2, y2, x0) = u
return [
x2^2 + y2^2 - (R - r2)^2, # eq1
x2^2 - (r1 + r2)^2 + (-R + r1 + y2)^2, # eq2
-R^2*r1^2 + R^2*x0^2 - 2*R*r1*x0^2 + r1^2*x0^2, # eq3
-R^2*r2^2 + R^2*x0^2 - 4*R*r1*x0^2 - 2*R*r2*x0^2 + 4*r1^2*x0^2 + 4*r1*r2*x0^2 + r2^2*x0^2, # eq4
(R^2*(-r2^2 + x2^2) - x0^2*x2^2 + x0^2*y2^2 - 2*x0*x2*y2*sqrt(R^2 - x0^2))/R^2, # eq5
]
end;
R = 30
iniv = BigFloat[11, 3, 13, 24, 17]
res = nls(H, ini=iniv)
([10.98076211353316, 2.942286340599478, 13.126775968941422, 23.660254037844386, 17.320508075688775], true)
たとえば，扇長が 30 寸のとき，大円の直径は 2*10.98076211353316 = 21.96152422706632 寸である。
「術」は「3 の平方根から 1 を引き，扇長を掛ける」と，実に簡明な答えが示されている。
上の数値解は正しく，正確なものであることがわかる。
2*10.98076211353316 |> println
(√3 - 1)*30 |> println
21.96152422706632
21.961524227066317
function draw(R, more)
pyplot(size=(700, 700), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, x2, y2, x0) = res[1]
y0 = sqrt(R^2 - x0^2)
θ = atand(y0, x0)
plot([-x0, 0, x0], [y0, 0, y0], color=:red, lw=0.5)
circle(0, R - r1, r1, :green)
circle(0, 0, R, beginangle=θ, endangle=180 - θ)
circle2(x2, y2, r2, :blue)
circle(0, R - 2r1 - r2, r2, :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 - r1, "大円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
point(0, R - 2r1 - r2, "等円:r2\n(0,R-2r1-r2)", :black, :center, :bottom, delta=delta/2)
point(x2, y2, "等円:r2\n(x2,y2)", :black, :center, :bottom, delta=delta/2)
point(0, R, "扇長:R", :red, :center, :bottom, delta=delta/2)
point(x0, sqrt(R^2 - x0^2), " (x0,√(R^2-x0^2))", :red, :left, :vcenter)
xlims!(-x0 - 2delta, x0 + 30delta)
end
end;
draw(30, true)
