算額（その782）
福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html
直角三角形の各辺を直径とする 3 個の半円と，それらの隙間に乾円 1 個，坤円 2 個を入れる。坤円の直径が与えられたとき，乾円の直径を求めよ。
直角三角形の直角を挟む二辺のうち短い方の長さを 2r1，長い方の長さを 2r2，対辺の長さを 2R とする。
乾円の半径と中心座標を r3, (r3, y3)
坤円の半径と中心座標を r4, (x4, y4), (r2, -r4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive, y3::positive,
r4::positive, x4::positive, y4::positive, x0::positive, y0::positive, d
R = sqrt(4r1^2 + 4r2^2)/2
eq1 = r1 + 2r4 - R
eq2 = x0^2 + (y0 - r1)^2 - r1^2
eq3 = y0/(2r2 - x0) - 2r1/2r2
eq4 = x4^2 + (r1 - y4)^2 - (r1 - r4)^2
eq5 = (r2 - r3)^2 + y3^2 - (r2 + r3)^2
eq6 = (r2 - x4)^2 + y4^2 - (r2 - r4)^2
eq7 = dist(0, 2r1, 2r2, 0, r3, y3) - r3^2
eq7 = div(numerator(apart(eq7, d)), r2);
eq8 = 2R*sqrt(x0^2 + y0^2) - 4r1*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, r3, y3, x4, y4, x0, y0) = u
return [
r1 + 2*r4 - sqrt(4*r1^2 + 4*r2^2)/2, # eq1
-r1^2 + x0^2 + (-r1 + y0)^2, # eq2
-r1/r2 + y0/(2*r2 - x0), # eq3
x4^2 - (r1 - r4)^2 + (r1 - y4)^2, # eq4
y3^2 + (r2 - r3)^2 - (r2 + r3)^2, # eq5
y4^2 - (r2 - r4)^2 + (r2 - x4)^2, # eq6
4*r1^2*r2 - 4*r1^2*r3 - 4*r1*r2*y3 + 2*r1*r3*y3 - r2*r3^2 + r2*y3^2, # eq7
-4*r1*r2 + sqrt(4*r1^2 + 4*r2^2)*sqrt(x0^2 + y0^2), # eq8
]
end;
r4 = 1/2
iniv = BigFloat[1.5, 2, 0.5, 2, 0.8, 0.9, 36/25, 48/25]
res = nls(H, ini=iniv);
println(res);
([1.5, 2.0, 0.5, 2.0, 0.8000000000009287, 0.9000000000012384, 1.44, 1.92], true)
坤円の半径が 0.5 のとき，乾円の半径は 0.5 である（坤円の半径に等しい）。
その他のパラメータは以下のとおりである。
r1 = 1.5; r2 = 2; r3 = 0.5; y3 = 2; x4 = 0.8; y4 = 0.9; x0 = 1.44; y0 = 1.92; R = 2.5
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r4 = 0.5
(r1, r2, r3, y3, x4, y4, x0, y0) = (1.5, 2, 0.5, 2, 0.8, 0.9, 36/25, 48/25)
(r1, r2, r3, y3, x4, y4, x0, y0) = res[1]
R = sqrt(4r1^2 + 4r2^2)/2
@printf("r1 = %g; r2 = %g; r3 = %g; y3 = %g; x4 = %g; y4 = %g; x0 = %g; y0 = %g; R = %g\n", r1, r2, r3, y3, x4, y4, x0, y0, R)
@printf("乾円の直径 = %g; 坤円の直径 = %g\n", r4, r3)
plot([0, 2r2, 0, 0], [0, 0, 2r1, 0], color=:black, lw=0.5)
circle(0, r1, r1, beginangle=270, endangle=450)
circle(r2, 0, r2, :blue, beginangle=0, endangle=180)
circle(r3, y3, r3, :green)
circle(x4, y4, r4, :magenta)
circle(r2, -r4, r4, :magenta)
circle(r2, r1, R, :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(x0, y0, " (x0,y0)", :black, :left, :vcenter)
point(r2, r1, " (r2,r1)", :black, :left, :vcenter)
point(r3, y3, "乾円:r3\n(r3,y3)", :green, :center, delta=-delta/2)
point(x4, y4, "坤円:r4\n(x4,y4)", :magenta, :center, delta=-delta/2)
point(r2, -r4, "坤円:r4\n(r2,-r4)", :magenta, :center, delta=-delta/2)
point(0, 2r1, "2r1 ", :red, :right, :vcenter)
point(0, r1, "r1 ", :red, :right, :vcenter)
point(r2, 0, " r2", :blue, :left, :bottom, delta=delta/2)
point(2r2, 0, " 2r2", :blue, :left, :bottom, delta=delta/2)
else
plot!(showaxis=false)
end
end;
