算額（その363）

算額（その363）

山形県寒河江市大字柴橋字台下 伊豆神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

直線上に大円，中円，小円が互いに接している。大円，中円の直径をそれぞれ 9 寸，4 寸としたとき，小円の直径を求めよ。

大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive;

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x2, r3, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*sqrt(r1)*sqrt(r2), (-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) - 2*r1^2*r2 + 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))
    (2*sqrt(r1)*sqrt(r2), (2*r1^(7/2)*r2^(3/2) - 4*r1^(5/2)*r2^(5/2) + 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) + 2*r1^2*r2 - 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))

最初の組のものが適解である。

r1, r2 が 4.5, 2 のとき，x2 = 6, r3 = 0.72, x3 = 3.6 である。

(r1, r2) = (9, 4) .// 2
(2*sqrt(r1)*sqrt(r2), (-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4), (2*r1^(5/2)*sqrt(r2) - 2*r1^(3/2)*r2^(3/2) - 2*r1^2*r2 + 2*r1*r2^2)/(r1^2 - 2*r1*r2 + r2^2))

   (6.0, 0.7199999999999942, 3.6)

直線上にあり，互いに接している 2 円間の水平距離については
「和算の心（その002）https://blog.goo.ne.jp/r-de-r/e/eeac1b6f4ed0f7c3a7dfdad299f066b3」
にも記したように，x2 = 2sqrt(r1*r2) = 6 で求めることができる。

2*sqrt(r1*r2)

   6.0

直線上にあり，互いに接している 3 円において，この場合の小円の半径については，デカルトの円定理により，1/(1/r1 + 1/r2 + 2sqrt(1/(r1*r2))) = 0.72 で求めることができる。

1/(1/r1 + 1/r2 + 2sqrt(1/(r1*r2)))

   0.72

なお，連立方程式の解で r3 は
(-2*r1^(7/2)*r2^(3/2) + 4*r1^(5/2)*r2^(5/2) - 2*r1^(3/2)*r2^(7/2) + r1^4*r2 - r1^3*r2^2 - r1^2*r2^3 + r1*r2^4)/(r1^4 - 4*r1^3*r2 + 6*r1^2*r2^2 - 4*r1*r2^3 + r2^4)
という長い式になっており，SymPy では簡約化されない。

大円と小円の中心の水平距離は x3 = 2sqrt(4.5 * 0.72) = 3.6 である。

2sqrt(4.5 * 0.72)

   3.6

r1, r2 を数値として与え連立方程式を解くと以下のようになる。

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive;
(r1, r2) = (9, 4) .// 2
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (x2, r3, x3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (6, 18/25, 18/5)
    (6, 18, 18)

(6, 18/25, 18/5)

   (6, 0.72, 3.6)

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (9, 4) .// 2
   (x2, r3, x3) = (6, 18/25, 18/5)
   @printf("x2 = %g;  r3 = %g;  x3 = %g\n", x2, r3, x3)
   plot()
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r1, " 大円:r1,(0,r1)", :red, :left, :vcenter)
       point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(x3, r3, "小円:r3,(x3,r3)", :magenta, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       hline!([0], color=:black, lw=0.5)
       plot!(showaxis=false)
   end
end;