算額（その772）
宮城県白石市小原 小原温泉薬師堂 大正5年(1916)
徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf
外円の中に大円，中円，小円，甲円，乙円が入っている。中円，小円，甲円の直径がそれぞれ 9 寸，6 寸，3 寸のとき，乙円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (x3, y3)
甲円の半径と中心座標を r4, (x4, y4)
乙円の半径と中心座標を r5, (x5, y5)
とおき，以下の連立方程式を解く。
算額にある図とは大幅に異なる結果になる。そもそも「大円」は「中円」はおろか，「小円」よりも小さい。
include("julia-source.txt");
using SymPy
@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x4^2 + (R - r1 - y4)^2 - (r1 + r4)^2
eq4 = x5^2 + (R - r1 - y5)^2 - (r1 + r5)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x3^2 + y3^2 - (R - r3)^2
eq7 = x4^2 + y4^2 - (R - r4)^2
eq8 = x5^2 + y5^2 - (R - r5)^2
eq9 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq10 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq11 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^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)
(R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = u
return [
x2^2 - (r1 + r2)^2 + (R - r1 - y2)^2, # eq1
x3^2 - (r1 + r3)^2 + (R - r1 - y3)^2, # eq2
x4^2 - (r1 + r4)^2 + (R - r1 - y4)^2, # eq3
x5^2 - (r1 + r5)^2 + (R - r1 - y5)^2, # eq4
x2^2 + y2^2 - (R - r2)^2, # eq5
x3^2 + y3^2 - (R - r3)^2, # eq6
x4^2 + y4^2 - (R - r4)^2, # eq7
x5^2 + y5^2 - (R - r5)^2, # eq8
-(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2, # eq9
-(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2, # eq10
-(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2, # eq11
]
end;
(r2, r3, r4) = (9, 6, 3) .// 2
iniv = BigFloat[7.5, 2.43, 2.7, -1.31, -4.2, 1.62, 3.9, 4.56, 1.0, -3.4, 5.54]
res = nls(H, ini=iniv)
([7.5031904642363045, 2.4327679290250375, 2.7, -1.3149726097831362, -4.2, 1.6244150815566774, 3.9, 4.563802772896491, 1.0, -3.4, 5.543598670009763], true)
乙円の半径は 1 寸である（直径は 2 寸）。
その他のパラメータは以下のとおりである。
R = 7.50319; r1 = 2.43277; x2 = 2.7; y2 = -1.31497; x3 = -4.2; y3 = 1.62442; x4 = 3.9; y4 = 4.5638; r5 = 1; x5 = -3.4; y5 =5.5436
function draw(more=false)
pyplot(size=(400, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3, r4) = (9, 6, 3) .// 2
(R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = res[1]
@printf("r2 = %g; r3 = %g; r4 = %g のとき，乙円の直径は %g\n", r2, r3, r4, 2r5)
@printf("R = %g; r1 = %g; x2 = %g; y2 = %g; x3 = %g; y3 = %g; x4 = %g; y4 = %g; r5 = %g; x5 = %g; y5 =%g\n", R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5)
plot()
circle(0, 0, R)
circle(0, R - r1, r1, :blue)
circle(x2, y2, r2, :green)
circle(x3, y3, r3, :magenta)
circle(x4, y4, r4, :purple)
circle(x5, y5, r5, :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(0, R - r1, "大円:r1\n(0,R-r1)", :blue, :center, :bottom, delta=delta)
point(x2, y2, "中円:r2\n(x2,y2)", :green, :center, :bottom, delta=delta)
point(x3, y3, "小円:r3\n(x3,y3)", :magenta, :center, :bottom, delta=delta)
point(x4, y4, "甲円:r4,(x4,y4)", :purple, :left, delta=-delta)
point(x5, y5, "乙円:r5,(x5,y5)", :black, :right, delta=-delta)
point(R, 0, " R", :red, :left, :bottom, delta=delta)
end
end;
