算額（その813）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
不等辺三角形内に全円，大円，中円，小円を入れる。大円，中円，小円の直径がそれぞれ 256 寸，225 寸，144 寸のとき，全円の直径はいかほどか。
底辺の長さ a，頂点の座標を (b, h)
全円の半径と中心座標を r0, (x0, r0)
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, h::positive,
r0::positive, x0::positive,
r1::positive, x1::positive,
r2::positive, x2::positive,
r3::positive, x3::positive, y3::positive,
d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x3 - x2)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = (a + sqrt(b^2 + h^2) + sqrt((a - b)^2 + h^2))r0 - a*h
eq5 = r0/(a - x0) - r1/(a - x1)
eq6 = r0/x0 - r2/x2
eq7 = dist(0, 0, b, h, x2, r2) - r2^2
eq8 = dist(0, 0, b, h, x3, y3) - r3^2
eq9 = dist(a, 0, b, h, x3, y3) - 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)
(a, b, h, x0, r0, x1, x2, x3, y3) = u
return [
(r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2, # eq1
(-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2, # eq2
(-r2 + y3)^2 - (r2 + r3)^2 + (-x2 + x3)^2, # eq3
-a*h + r0*(a + sqrt(b^2 + h^2) + sqrt(h^2 + (a - b)^2)), # eq4
r0/(a - x0) - r1/(a - x1), # eq5
r0/x0 - r2/x2, # eq6
-r2^2 + (-b*(b*x2 + h*r2)/(b^2 + h^2) + x2)^2 + (-h*(b*x2 + h*r2)/(b^2 + h^2) + r2)^2, # eq7
-r3^2 + (-b*(b*x3 + h*y3)/(b^2 + h^2) + x3)^2 + (-h*(b*x3 + h*y3)/(b^2 + h^2) + y3)^2, # eq8
-r3^2 + (-h*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2) + y3)^2 + (-a + x3 - (-a + b)*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2))^2, # eq9
]
end;
(r1, r2, r3) = (256, 225, 144) .// 2
iniv = BigFloat[1000, 355, 358, 385, 162, 513, 269, 368, 269]
res = nls(H, ini=iniv)
res |> println
2res[1][5] |> println
([1014.0, 354.88757396449705, 357.8698224852071, 384.0, 160.0, 510.0, 270.0, 367.9881656804734, 268.8284023668639], true)
320.0
大円，中円，小円の直径がそれぞれ 256 寸，225 寸，144 寸のとき，全円の直径は 320 寸である。
その他のパラメータは以下のとおりである。
a = 1014; b = 354.888; h = 357.87; x0 = 384; r0 = 160; x1 = 510; x2 = 270; x3 = 367.988; y3 = 268.828
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (256, 225, 144) .// 2
(a, b, h, x0, r0, x1, x2, x3, y3) = res[1]
@printf("全円の直径 = %g\n", 2r0)
@printf("a = %g; b = %g; h = %g; x0 = %g; r0 = %g; x1 = %g; x2 = %g; x3 = %g; y3 = %g\n", a, b, h, x0, r0, x1, x2, x3, y3)
plot([0, a, b, 0], [0, 0, h, 0], color=:black, lw=0.5)
circle(x0, r0, r0, :orange)
circle(x1, r1, r1, :blue)
circle(x2, r2, r2, :green)
circle(x3, y3, r3)
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(x0, r0, " 全円:r0,(x0,r0)", :black, :left, :vcenter)
point(x1, r1, "大円:r1\n(x1,r1)", :blue, :center, delta=-delta)
point(x2, r2, "中円:r2\n(x2,r2)", :green, :center, delta=-delta)
point(x3, y3, " 小円:r3,(x3,y3)", :red, :left, :vcenter)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(b, h, "(b,h)", :black, :left, :bottom, delta=delta)
end
end;
