算額（その575）
福岡県久留米高良山者
藤田貞資(1789): 神壁算法
http://hyonemitsu.web.fc2.com/Korataisha.pdf
外円内に斜線を隔てて甲円 1 個，乙円 2 個，丙円 1 個がある。外円，甲円，丙円の直径が 60 寸，20 寸，28 寸のとき，乙円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x, y)
丙円の半径と中心座標を r3, (0, r3 - R)
斜線と外円の交点座標を (x1, sqrt(R^2 - x1^2)), (x2, -sqrt(R^2 - x2^2))
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, r2::positive, x::positive, y::positive,
r3::positive, x1::positive, x2::positive
(R, r1, r3) = (60, 20, 28) .// 2
eq1 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, R - r1) - r1^2
eq2 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), x, y) - r2^2
eq3 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, r3 - R) - r3^2
eq4 = distance(-x2, -sqrt(R^2 - x2^2), x1, sqrt(R^2 - x1^2), x, y) - r2^2
eq5 = x^2 + y^2 - (R - r2)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x1, x2))
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)
(r2, x, y, x1, x2) = u
return [
(20 - (-x1^2*sqrt(900 - x2^2) - 20*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) - 20*x2^2 + 40*sqrt(900 - x1^2)*sqrt(900 - x2^2) + 36000)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) - 20*sqrt(900 - x1^2) - 20*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 100, # eq1
-r2^2 + (x - (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (2*sqrt(900 - x1^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) - (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2, # eq2
(-16 - (-x1^2*sqrt(900 - x2^2) + 16*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) + 16*x2^2 - 32*sqrt(900 - x1^2)*sqrt(900 - x2^2) - 28800)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 16*sqrt(900 - x1^2) + 16*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 196, # eq3
-r2^2 + (x - (-x2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (-2*sqrt(900 - x2^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2, # eq4
x^2 + y^2 - (30 - r2)^2, # eq5
]
end;
iniv = BigFloat[29, 40, 12, 21, 28] .* (30/70)
res = nls(H, ini=iniv)
(BigFloat[12.50000000000000000000000000000000000000000000000000000001913542531009364124636, 16.77050983124842272306880251548457176580463769708644293204781524730502299339952, 4.999999999999999999999999999999999999999999999999999999963738320809945557470043, 17.39163982499836430540468409013214849787147613031186674435362589343567963493901, 22.3606797749978969640917366873127623544061835961152572427143001079596068757832], true)
乙円の直径は 25 寸である。
その他のパラメータは以下の通り。
乙円の直径 = 25; r2 = 12.5; x = 16.7705; y = 5; x1 = 17.3916; x2 = 22.3607
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r3) = (60, 20, 28) .// 2
(r2, x, y, x1, x2) = res[1]
@printf("乙円の直径 = %g; r2 = %g; x = %g; y = %g; x1 = %g; x2 = %g\n", 2r2, r2, x, y, x1, x2)
plot()
circle(0, 0, R)
circle(0, R - r1, r1, :blue)
circle(x, y, r2, :magenta)
circle(-x, y, r2, :magenta)
circle(0, r3 - R, r3, :green)
segment(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), :gray)
segment(x1, sqrt(R^2 - x1^2), -x2, -sqrt(R^2 - x2^2), :gray)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(0, R - r1, " 甲円:r1\n (0,R-r1)", :blue, :left, :vcenter)
point(0, r3 - R, " 丙円:r1\n (0,r3-R)", :green, :left, :vcenter)
point(x, y, "乙円:r2,(x,y)", :magenta, :center, :top, delta=-delta/2)
point(x1, √(R^2 - x1^2), "(x1,√(R^2 - x1^2))", :black, :center, :bottom, delta=delta/2)
point(x2, -√(R^2 - x2^2), "(x2,-√(R^2 - x2^2))", :black, :center, :bottom, delta=delta/2)
end
end;
