算額（その608）
高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/11
半円内に2本の斜線を引き，区画された領域に甲円，乙円，丙円を入れる。乙円，丙円の直径が与えられたとき，甲円の直径を求めよ。
算額（その511）より一段階複雑になっている。
半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
2 本の斜線が半円の円周上で交差する座標を (x, y)
とおき以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, x::negative, y::positive, r1::positive, x1::negative, r2::positive, x2::positive, y2::positive, r3::positive, x3::negative, y3::positive
y = sqrt(R^2 - x^2)
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = distance(R, 0, x, y, x1, r1) - r1^2
eq4 = y/(x + R) * y3/x3 + 1
eq5 = y/(x - R) * y2/x2 + 1
eq6 = (sqrt((x + R)^2 + y^2) + sqrt((R - x)^2 + y^2) + 2R)*r1 - 2R*y
eq7 = distance(R, 0, x, y, x2, y2) - r2^2
eq8 = distance(-R, 0, x, y, 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 v, r.f_converged
end;
function H(u)
(R, x, r1, x1, x2, y2, x3, y3) = u
return [
x2^2 + y2^2 - (R - r2)^2, # eq1
x3^2 + y3^2 - (R - r3)^2, # eq2
-r1^2 + (r1 - (R*r1 + R*sqrt(R^2 - x^2) + r1*x - x1*sqrt(R^2 - x^2))/(2*R))^2 + (x1 - (R^2 + R*x + R*x1 - r1*sqrt(R^2 - x^2) - x*x1)/(2*R))^2, # eq3
1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)), # eq4
1 + y2*sqrt(R^2 - x^2)/(x2*(-R + x)), # eq5
-2*R*sqrt(R^2 - x^2) + r1*(2*R + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2)), # eq6
-r2^2 + (x2 - (R^2 + R*x + R*x2 - x*x2 - y2*sqrt(R^2 - x^2))/(2*R))^2 + (y2 - (R*y2 + R*sqrt(R^2 - x^2) + x*y2 - x2*sqrt(R^2 - x^2))/(2*R))^2, # eq7
-r3^2 + (x3 - (-R^2 + R*x + R*x3 + x*x3 + y3*sqrt(R^2 - x^2))/(2*R))^2 + (y3 - (R*y3 + R*sqrt(R^2 - x^2) - x*y3 + x3*sqrt(R^2 - x^2))/(2*R))^2, # eq8
]
end;
(r2, r3) = (10, 5)
iniv = BigFloat[50, -10, 20, -7, 23, 33, -35, 27]
res = nls(H, ini=iniv)
(BigFloat[50.00000000000000000000000000000000000000000000000000000000000000000000000630095, -14.00000000000000000000000000000000000000000000000000000000000000000000000571507, 20.00000000000000000000000000000000000000000000000000000000000000000000000238607, -10.00000000000000000000000000000000000000000000000000000000000000000000000404504, 24.00000000000000000000000000000000000000000000000000000000000000000000000337377, 32.00000000000000000000000000000000000000000000000000000000000000000000000631642, -36.00000000000000000000000000000000000000000000000000000000000000000000000619483, 27.00000000000000000000000000000000000000000000000000000000000000000000000251126], true)
乙円，丙円の直径が 10, 5 のとき，甲円の直径は 40 である。
その他のパラメータは以下の通り。
R = 50; x = -14; y = 48; r1 = 20; x1 = -10; r2 = 10; x2 = 24; y2 = 32; r3 = 5; x3 = -36; y3 = 27
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (10, 5)
(R, x, r1, x1, x2, y2, x3, y3) = res[1]
y = sqrt(R^2 - x^2)
@printf("甲円の直径 = %g\n", 2r1)
@printf("R = %g; x = %g; y = %g; r1 = %g; x1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; x3 = %g; y3 = %g\n",
R, x, y, r1, x1, r2, x2, y2, r3, x3, y3)
plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
circle(0, 0, R, :black, beginangle=0, endangle=180)
circle(x1, r1, r1)
circle(x2, y2, r2, :blue)
circle(x3, y3, r3, :green)
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(x1, r1, "甲円:r1 \n(x1,r1) ", :red, :right, :vcenter)
point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, :top, delta=-delta)
point(x3, y3, " 丙円:r3,(x3,y3)", :green, :left, :bottom, delta=delta)
point(0, R, " R", :black, :left, :bottom, delta=delta)
end
end;
