算額(その785)
福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html
半円内に 2 本の斜線を描き,天円,地円,人円を入れる。人円の直径が与えられたとき,天円の直径を求めよ。
半円の半径と中心座標を r1, (0, 0)
天円の半径と中心座標を r2, (0, r2)
地円の半径と中心座標を r3, (x3, r3)
人円の半径と中心座標を r4, (x4, r4)
斜線と半円の交点座標を (x01, y01), (x02, y02)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, r2::positive,r3::positive, x3::positive,
r4::positive, x4::positive, y4::positive,
x01::positive, y01::positive, x02::positive, y02::positive
eq1 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq2 = x3^2 + r3^2 - (r1 - r3)^2
eq3 = x4^2 + y4^2 - (r1 - r4)^2
eq4 = 2sqrt(r1^2 - (r1 - 2r4)^2) - sqrt((x01 - x02)^2 + (y01 - y02)^2)
eq5 = x01^2 + y01^2 - r1^2
eq6 = x02^2 + y02^2 - r1^2
eq7 = (x01 - x02)/(y02 - y01) - y4/x4
eq8 = r1 - 2r2
eq9 = dist(x01, y01, x02, y02, 0, r2) - r2^2
eq10 = dist(x01, y01, x02, y02, x3, r3) - r3^2;
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, x3, x4, y4, x01, y01, x02, y02) = u
return [
x3^2 + (r2 - r3)^2 - (r2 + r3)^2, # eq1
r3^2 + x3^2 - (r1 - r3)^2, # eq2
x4^2 + y4^2 - (r1 - r4)^2, # eq3
2*sqrt(r1^2 - (r1 - 2*r4)^2) - sqrt((x01 - x02)^2 + (y01 - y02)^2), # eq4
-r1^2 + x01^2 + y01^2, # eq5
-r1^2 + x02^2 + y02^2, # eq6
(x01 - x02)/(-y01 + y02) - y4/x4, # eq7
r1 - 2*r2, # eq8
-r2^2 + (-x01 - (-x01 + x02)*(-x01*(-x01 + x02) + (r2 - y01)*(-y01 + y02))/((-x01 + x02)^2 + (-y01 + y02)^2))^2 + (r2 - y01 - (-y01 + y02)*(-x01*(-x01 + x02) + (r2 - y01)*(-y01 + y02))/((-x01 + x02)^2 + (-y01 + y02)^2))^2, # eq9
-r3^2 + (r3 - y01 - (-y01 + y02)*((r3 - y01)*(-y01 + y02) + (-x01 + x02)*(-x01 + x3))/((-x01 + x02)^2 + (-y01 + y02)^2))^2 + (-x01 + x3 - (-x01 + x02)*((r3 - y01)*(-y01 + y02) + (-x01 + x02)*(-x01 + x3))/((-x01 + x02)^2 + (-y01 + y02)^2))^2, # eq10
]
end;
r4 = 1/2
iniv = BigFloat[9.1, 4.6, 2.2, 6.3, 5.3, 6.6, 8.2, 3.6, 1.8, 8.8]
res = nls(H, ini=iniv)
([9.0, 4.5, 2.25, 6.363961030678928, 5.342584568965026, 6.611111111111111, 8.23517481947363, 3.630688046735422, 1.8214549574017131, 8.813756397709023], true)
人円の直径が 1 のとき,天円の直径は 9 である。
その他のパラメータは以下のとおりである。
r1 = 9; r2 = 4.5; r3 = 2.25; x3 = 6.36396; x4 = 5.34258; y4 = 6.61111; x01 = 8.23517; y01 = 3.63069; x02 = 1.82145; y02 = 8.81376
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r4 = 1/2
(r1, r2, r3, x3, x4, y4, x01, y01, x02, y02) = res[1]
@printf("人円の直径が %g のとき,天円の直径は %g\n", 2r4, 2r2)
@printf("r1 = %g; r2 = %g; r3 = %g; x3 = %g; x4 = %g; y4 = %g; x01 = %g; y01 = %g; x02 = %g; y02 = %g\n",
r1, r2, r3, x3, x4, y4, x01, y01, x02, y02)
plot()
circle(0, 0, r1, beginangle=0, endangle=180)
circle(0, r2, r2, :blue)
circle2(x3, r3, r3, :green)
circle2(x4, y4, r4, :magenta)
segment(x01, y01, x02, y02, :orange)
segment(-x01, y01, -x02, y02, :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(r1, 0, "r1 ", :red, :right, :bottom, delta=delta/2)
point(0, 0, "半円:r1,(0,0)", :red, :center, :bottom, delta=delta/2)
point(r1, 0, "r1 ", :red, :right, :bottom, delta=delta/2)
point(0, r2, "天円:r2,(0,r2)", :blue, :center, :bottom, delta=delta)
point(x3, r3, "地円:r3,(x3,r3)", :green, :center, :bottom, delta=delta)
point(x4, y4, "人円:r4,(x4,r4) ", :black, :right, :vcenter)
point(x01, y01, "(x01,y01)", :black, :right, :vcenter)
point(x02, y02, "(x02,y02)", :black, :left, :bottom, delta=delta)
end
end;