算額（その830）
宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf
直線の上に甲乙丙丁戊己庚の 7 個の円が載って隣同士外接している。庚円の直径が 3 寸のとき，庚円の頂点までの高さはいかほどか。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, r5)
己円の半径と中心座標を r6, (0, r6)
庚円の半径と中心座標を r7, (x7, y7)
庚円の頂点までの高さを h
として，以下の連立方程式を解く。
なお，eq13, eq14 は算法助術の公式47 を援用した。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive,
r3::positive, x3::positive, r4::positive, x4::positive,
r5::positive, x5::negative, r6::positive,
r7::positive, x7::negative, y7::positive, h::positive
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = x1^2 + (r1 - r6)^2 - (r1 + r6)^2 |> expand
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand
eq5 = (x4 - x5)^2 + (r4 - r5)^2 - (r4 + r5)^2 |> expand
eq6 = x5^2 + (r5 - r6)^2 - (r5 + r6)^2 |> expand
eq7 = x7^2 + (r6 - y7)^2 - (r6 + r7)^2 |> expand
eq8 = (x1 - x7)^2 + (r1 - y7)^2 - (r1 + r7)^2 |> expand
eq9 = (x2 - x7)^2 + (r2 - y7)^2 - (r2 + r7)^2 |> expand
eq10 = (x3 - x7)^2 + (r3 - y7)^2 - (r3 + r7)^2 |> expand
eq11 = (x4 - x7)^2 + (r4 - y7)^2 - (r4 + r7)^2 |> expand
eq12 = (x5 - x7)^2 + (r5 - y7)^2 - (r5 + r7)^2 |> expand
eq13 = r2*h + r3*h - 2r2*r3 - 2sqrt(2r2*r3*r7*h) |> simplify
eq14 = r4*h + r3*h - 2r4*r3 - 2sqrt(2r4*r3*r7*h) |> simplify;
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, x1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x7, y7, h) = u
return [
-4*r1*r2 + x1^2 - 2*x1*x2 + x2^2, # eq1
-4*r1*r6 + x1^2, # eq2
-4*r2*r3 + x2^2 - 2*x2*x3 + x3^2, # eq3
-4*r3*r4 + x3^2 - 2*x3*x4 + x4^2, # eq4
-4*r4*r5 + x4^2 - 2*x4*x5 + x5^2, # eq5
-4*r5*r6 + x5^2, # eq6
-2*r6*r7 - 2*r6*y7 - r7^2 + x7^2 + y7^2, # eq7
-2*r1*r7 - 2*r1*y7 - r7^2 + x1^2 - 2*x1*x7 + x7^2 + y7^2, # eq8
-2*r2*r7 - 2*r2*y7 - r7^2 + x2^2 - 2*x2*x7 + x7^2 + y7^2, # eq9
-2*r3*r7 - 2*r3*y7 - r7^2 + x3^2 - 2*x3*x7 + x7^2 + y7^2, # eq10
-2*r4*r7 - 2*r4*y7 - r7^2 + x4^2 - 2*x4*x7 + x7^2 + y7^2, # eq11
-2*r5*r7 - 2*r5*y7 - r7^2 + x5^2 - 2*x5*x7 + x7^2 + y7^2, # eq12
-2*sqrt(2h*r2*r3*r7) + h*r2 + h*r3 - 2*r2*r3, # eq13
-2*sqrt(2h*r3*r4*r7) + h*r3 + h*r4 - 2*r3*r4, # eq14
]
end;
r7 = 1.5
iniv = BigFloat[6.9429, 14.95, 0.9806, 9.7315, 0.5331, 8.2854, 0.5388, 7.2135, 1.0202, 5.7307, 8.0479, 7.7707, 2.5, 4.0]
res = nls(H, ini=iniv)
([6.942751374957064, 14.95003944293182, 0.9806016370341291, 9.731582443471599, 0.5331209415506033, 8.28551291814206, 0.538803297760604, 7.213603740249264, 1.0201813325955165, 5.730799550457136, 8.048094598036768, 7.770762947374868, 2.5, 4.0], true)
庚円の直径が 3 寸のとき，庚円の頂点までの高さは 4 寸である。
その他のパラメータは以下のとおりである。
r1 = 6.94275; x1 = 14.95; r2 = 0.980602; x2 = 9.73158; r3 = 0.533121; x3 = 8.28551; r4 = 0.538803; x4 = 7.2136; r5 = 1.02018; x5 = 5.7308; r6 = 8.04809; x7 = 7.77076; y7 = 2.5; h = 4
function draw(more)
pyplot(size=(700, 700), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
r7 = 1.5
(r1, x1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x7, y7, h) = res[1]
@printf("庚円の頂点の高さ = %g\n", h)
@printf("r1 = %g; x1 = %g; r2 = %g; x2 = %g; r3 = %g; x3 = %g; r4 = %g; x4 = %g; r5 = %g; x5 = %g; r6 = %g; x7 = %g; y7 = %g; h = %g\n", r1, x1, r2, x2, r3, x3, r4, x4, r5, x5, r6, x7, y7, h)
plot()
circle(x1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :green)
circle(x4, r4, r4, :magenta)
circle(x5, r5, r5, :orange)
circle(0, r6, r6, :tomato)
circle(x7, y7, r7, :brown)
segment(-r6, 0, x1 + r1, 0, :black,lw=1)
if more == true
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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta/2)
point(x2, r2, "乙円:r2,(x2,r2)", :blue, :left, :bottom, delta=delta/2)
point(x3, r3, " 丙円:r3,(x3,r3)", :green, :left, :vcenter)
point(x4, r4, "丁円:r4,(x4,r4) ", :magenta, :right, :vcenter)
point(x5, r5, "戊円:r5,(x5,r5) ", :black, :right, :bottom, delta=delta/2)
point(0, r6, "己円:r6,(0,r6) ", :brown, :center, delta=-delta/2)
point(x7, y7, "庚円:r7,(x7,y7) ", :brown, :center, :bottom, delta=delta/2)
point(x7, h, " (x7,h)", :black, :left, :bottom, delta=delta/2)
end
end;
