算額（その446）
寛政7年乙卯10月(1795) 久留米 鈴木半平正晟
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
外円内に甲，乙，丙，丁の 4 種類の円が入っている。それぞれ隣り合う円に外接し，外円に内接している。甲円の直径が 3寸0分5厘のとき，外円の直径はいかほどか。
外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (0, r0 - 2r1 - r2, 0)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (r4, 0)
とおき，連立方程式の解を求める。
include("julia-source.txt")
using SymPy
@syms r0::positive, r1::positive,
r2::positive, r3::positive, x3::positive, y3::positive,
r4::positive;
r1 = 305//200
r4 = r0//2
eq1 = x3^2 + (r0 - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + (y3 - r0 +2r1 + r2)^2 - (r2 + r3)^2
eq3 = r4^2 + (r0 - 2r1 - r2)^2 - (r2 + r4)^2
eq4 = (x3 - r4)^2 + y3^2 - (r3 + r4)^2
eq5 = x3^2 + y3^2 - (r0 - r3)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r2, r3, x3, y3));
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)
(r0, r2, r3, x3, y3) = u
return [
x3^2 - (r3 + 61/40)^2 + (r0 - y3 - 61/40)^2, # eq1
x3^2 - (r2 + r3)^2 + (-r0 + r2 + y3 + 61/20)^2, # eq2
r0^2/4 - (r0/2 + r2)^2 + (r0 - r2 - 61/20)^2, # eq3
y3^2 + (-r0/2 + x3)^2 - (r0/2 + r3)^2, # eq4
x3^2 + y3^2 - (r0 - r3)^2, # eq5
]
end;
r1 = 305//200
iniv = [big"10.0", 2, 2, 2, 7]
res = nls(H, ini=iniv);
names = ("r0", "r2", "r3", "x3", "y3")
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;
r0 = 9.51; r2 = 1.86053; r3 = 2.05932; x3 = 3.33205; y3 = 6.66409;
外円の直径は 2r0 = 19.02 つまり，19寸0分2厘。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 305//200
(r0, r2, r3, x3, y3) = [24.0, 5, 5, 9, 17]
(r0, r2, r3, x3, y3) = res[1]
r4 = r0/2
@printf("r0 = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g\n", r0, r1, r2, r3, x3, y3)
@printf("外円の直径 = %g\n", 2r0)
plot()
circle(0, 0, r0, :red)
circle(0, r0 - r1, r1, :blue)
circle(0, r1 - r0, r1, :blue)
circle(0, r0 - 2r1 - r2, r2, :brown)
circle(0, 2r1 + r2 - r0, r2, :brown)
circle4(x3, y3, r3, :green)
circle(r4, 0, r4, :orange)
circle(-r4, 0, r4, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(0, r0 - r1, " 甲円:r1\n r0-r1", :black, :left, :vcenter)
point(0, r0 - 2r1 - r2, " 乙円:r2\n r0-2r1-r2", :black, :left, :vcenter)
point(x3, y3, " 丙円:r3\n (x3,y3)", :black, :left, :vcenter)
point(r4, 0, "丁円:r4,(r4,0)", :orange, :center, :bottom, delta=delta/2)
point(0, r0, " r0", :red, :left, :bottom, delta=delta/2)
hline!([0], color=:gray, lw=0.5)
vline!([0], color=:gray, lw=0.5)
else
plot!(showaxis=false)
end
end;
