算額（その466）
愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)
平田浩一: 伊佐爾波神社の算額にみる江戸末期の和算，愛媛大学教育学部紀要，第60巻，195-206, 2013.
https://www.ed.ehime-u.ac.jp/~kiyou/2013/pdf/20.pdf
円弧内に青，黃，赤，白，黒の 5 個の円がある。青，赤，黒の 3 個の円の直径が与えられたとき，白円の直径を求めよ。
円弧の半径と中心座標を r0, (x0, y0)
青円の半径と中心座標を r1, (x1, y1)
黃円の半径と中心座標を r2, (x2, y2)
赤円の半径と中心座標を r3, (x3, y3); y3 = r3
白円の半径と中心座標を r4, (x4, y4); y4 = r4
黒円の半径と中心座標を r5, (x5, y5); x5 = 0, y5 = r5
とおき，以下の連立方程式を解き，数値解を求める。
include("julia-source.txt")
using SymPy
@syms r1::positive, x1::positive, y1::positive,
r2::positive, x2::positive, y2::positive,
r3::positive, x3::positive,
r4::positive, x4::positive,
r5::positive,
r0::positive, x0::positive, y0::positive;
x5 = 0
y5 = r5
y3 = r3
y4 = r4
eq1 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq3 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq6 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq7 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = (x0 - x1)^2 + (y0 - y1)^2 - (r0 - r1)^2
eq10 = (x0 - x3)^2 + (y0 - y3)^2 - (r0 - r3)^2
eq11 = (x0 - x5)^2 + (y0 - y5)^2 - (r0 - r5)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0))
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 v, r.f_converged
end;
function H(u)
(x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = u
return [
-(r1 + r2)^2 + (x1 - x2)^2 + (y1 - y2)^2, # eq1
-(r1 + r3)^2 + (-r3 + y1)^2 + (x1 - x3)^2, # eq2
x1^2 - (r1 + r5)^2 + (-r5 + y1)^2, # eq3
-(r2 + r3)^2 + (-r3 + y2)^2 + (x2 - x3)^2, # eq4
(r3 - r4)^2 - (r3 + r4)^2 + (x3 - x4)^2, # eq5
-(r2 + r4)^2 + (-r4 + y2)^2 + (x2 - x4)^2, # eq6
x2^2 - (r2 + r5)^2 + (-r5 + y2)^2, # eq7
x4^2 + (r4 - r5)^2 - (r4 + r5)^2, # eq8
-(r0 - r1)^2 + (x0 - x1)^2 + (y0 - y1)^2, # eq9
-(r0 - r3)^2 + (-r3 + y0)^2 + (x0 - x3)^2, # eq10
x0^2 - (r0 - r5)^2 + (-r5 + y0)^2, # eq11
]
end;
青，赤，黒の 3 個の円の直径がそれぞれ 15, 30, 18 のときの解を求める。
(r1, r3, r5) = (15, 30, 18)
iniv = (r1/31) .* [big"15.0", 45, 6, 21, 25, 62, 10, 25, 57, 37, -15]
iniv = r1 .* [big"7.258", 21.774, 2.903, 10.161, 12.097, 30.0, 4.839, 12.097, 27.581, 17.903, -7.258]
res = nls(H, ini=iniv)
names = ("x1", "y1", "r2", "x2", "y2", "x3", "r4", "x4", "r0", "x0", "y0")
if res[2]
for (name, x) in zip(names, res[1])
@printf("%s = %g; ", name, round(Float64(x), digits=6))
end
println()
else
println("収束していない")
end;
@printf("白円の直径は %g である。\n", 2res[1][7])
x1 = 17.5989; y1 = 45.9156; r2 = 6.54237; x2 = 23.4663; y2 = 25.1876; x3 = 59.6904; r4 = 9.42818; x4 = 26.0544; r0 = 59.7271; x0 = 40.362; y0 = 7.41423;
白円の直径は 18.8564 である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0) = res[1]
@printf("x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g; x3 = %g; r4 = %g; x4 = %g; r0 = %g; x0 = %g; y0 = %g\n", x1, y1, r2, x2, y2, x3, r4, x4, r0, x0, y0)
@printf("白円の直径は %g である。\n", 2r4)
plot()
circle(x1, y1, r1)
circle(x2, y2, r2)
circle(x3, r3, r3)
circle(x4, r4, r4)
circle(0, r5, r5)
a = sqrt(r0^2 - y0^2)
segment(x0 - a, 0, x0 + a, 0, :magenta, lw=0.5)
θ = atand(y0/a)
circle(x0, y0, r0, :magenta, beginangle=-θ, endangle=181 + θ)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(x1, y1, " 青:r1\n (x1,y1)", :blue, :center, :bottom, delta=delta)
point(x2, y2, " 黃:r2,(x2,y2)", :black, :left, :vcenter)
point(x3, r3, " 赤:r3,(x3,r3)", :red, :center, :bottom, delta=delta)
point(x4, r4, "白:r4,(x4,r4)", :black, :center, :top, delta=-delta)
point(x5, r5, "黒:r5,(x5,r5)", :black, :center, :top, delta=-delta)
point(x0, y0, " 円弧:r0,(x0,y0)", :magenta, :left, :bottom, delta=delta)
hline!([0], color=:gray, lw=0.5)
vline!([0], color=:gray, lw=0.5)
else
plot!(showaxis=false)
end
end;
