算額（その502）

算額（その502）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8

甲円，乙円，丙円が図のように外接しあっている。甲円，乙円の直径が与えられたとき，丙円の直径を求めよ。

乙円の半径と中心座標を r2, (0, 0 )
甲円の半径と中心座標を r1, (r1, y1)
丙円の半径と中心座標を r3, (0, r2 + r3)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r1, y1, r2, r3

eq1 = r1^2 + (y1 - r2 - r3)^2 - (r1 + r3)^2
eq2 = r1^2 + y1^2 - (r1 + r2)^2;
res = solve([eq1, eq2], (y1, r3));

2 組の解が得られるが，2 番目のものが適解である。
簡約化すると以下のようになる。

# y1
res[1][1] |> simplify |> println

   sqrt(r2^3*(2*r1 + r2))/r2

# r3: 乙円の半径
res[1][2] |> simplify |> println

   (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2)

   r1 = 5;  r2 = 7;  y1 = 10.9087;  r3 = 0.857477
   甲円の直径が 10，乙円の直径が 14 のとき，丙円の直径は 1.71495 である。

using Plots

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y1, r3) = (sqrt(r2^3*(2*r1 + r2))/r2, (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2))
   @printf("r1 = %g;  r2 = %g;  y1 = %g;  r3 = %g\n", r1, r2, y1, r3)
   @printf("甲円の直径が %g，乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
   plot()
   circle(r1, y1, r1)
   circle(-r1, y1, r1)
   circle(0, 0, r2, :blue)
   circle(0, r2 + r3, 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(r1, y1, " 甲円:r1\n (r1,y1)", :red, :left, :vcenter)
       point(0, 0.4r2, " 乙円:r2,(0,0)", :blue, :left, mark=false)
       point(0, r2 + r3, " 丙円:r3,(0,r2+r3)", :green, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;