算額（その931）
群馬県安中市下磯部 磯部神社 天保7年(1830)
「算額」第三集 全国調査，香川県算額研究会
台形内に大円，小円，等円が入っている。大円の直径が 10 寸のとき，等円の直径はいかほどか。
横倒しになっているが，台形の高さ，下底，上底の長さを a, b, c
大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (a - r2, r2)
等円の半径と中心座標を r3, (r3, y31), (x3, r3)
とおき，以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, c::positive,
r1::positive, r2::positive, r3::positive,
x3::positive, y31::positive, y32::positive
eq1 = (r1 - r3)^2 + (y31 - r1)^2 - (r1 + r3)^2
eq2 = (x3 - r1)^2 + (y32 - r1)^2 - (r1 + r3)^2
eq3 = (a - r2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (a - r2 - x3)^2 + (y32 - r2)^2 - (r2 + r3)^2
eq5 = dist2(a, c, 0, b, r3, y31, r3)
eq6 = dist2(a, c, 0, b, x3, y32, r3)
eq7 = dist2(a, c, 0, b, r1, r1, r1)
eq8 = dist2(a, c, 0, b, a - r2, r2, 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)
(a, b, c, r2, r3, x3, y31, y32) = u
return [
(-r1 + y31)^2 + (r1 - r3)^2 - (r1 + r3)^2, # eq1
(-r1 + x3)^2 + (-r1 + y32)^2 - (r1 + r3)^2, # eq2
(r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - r2)^2, # eq3
(-r2 + y32)^2 - (r2 + r3)^2 + (a - r2 - x3)^2, # eq4
a^2*b^2 - 2*a^2*b*y31 - a^2*r3^2 + a^2*y31^2 - 2*a*b^2*r3 + 2*a*b*c*r3 + 2*a*b*r3*y31 - 2*a*c*r3*y31, # eq5
a^2*b^2 - 2*a^2*b*y32 - a^2*r3^2 + a^2*y32^2 - 2*a*b^2*x3 + 2*a*b*c*x3 + 2*a*b*x3*y32 - 2*a*c*x3*y32 - b^2*r3^2 + b^2*x3^2 + 2*b*c*r3^2 - 2*b*c*x3^2 - c^2*r3^2 + c^2*x3^2, # eq6
a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*r1 + 2*a*b*c*r1 + 2*a*b*r1^2 - 2*a*c*r1^2, # eq7
a^2*c^2 - 2*a^2*c*r2 + 2*a*b*c*r2 - 2*a*b*r2^2 - 2*a*c^2*r2 + 2*a*c*r2^2, # eq8
]
end;
r1 = 10/2
iniv = BigFloat[18, 11, 7, 4, 1, 10, 10, 8]
res = nls(H, ini=iniv)
([18.121162012570196, 11.12823119917143, 7.415656960521699, 4.083745164389245, 1.1267889889060243, 10.427965642765756, 9.747186511832085, 7.841607343658628], true)
大円の直径が 10 寸のとき，等円の直径は 2*1.1267889889060243 = 2.2535779778120486 寸である。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 10/2
(a, b, c, r2, r3, x3, y31, y32) = res[1]
println("2r3 = $(2r3)")
plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:green, lw=0.5)
circle(r1, r1, r1, :blue)
circle(a - r2, r2, r2, :green)
circle(r3, y31, r3, :magenta)
circle(x3, y32, r3, :magenta)
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(a, 0, " a", :green, :left, :bottom, delta=delta/2)
point(a, c, "(a,c)", :green, :left, :bottom, delta=2delta, deltax=-7delta)
point(0, b, " b", :green, :left, :bottom, delta=delta/2)
point(r1, r1, "大円:r1,(r1,r1)", :blue, :center, delta=-2delta)
point(a - r2, r2, "小円:r2,(a-r2,r2)", :green, :center, delta=-2delta)
point(r3, y31, "等円:r3\n(r3,y31)", :black, :left, delta=-2delta, deltax=-5delta)
point(x3, y32, "等円:r3\n(x3,y32)", :black, :left, delta=-2delta, deltax=-5delta)
end
end;
