算額（その608）

算額（その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;