和算の心（その008）

和算の心（その008）

#Julia, #SymPy, #算額, #和算

直線 2 本に挟まれ，互いに外接する 3 個の円の位置関係

図のように，直線（接線） 2 本に甲円，乙円，丙円が接している。
甲円と丙円の直径（2r1, 2r3）が分かっているとき，2 直線の交点 O と甲円と直線の接点 A, B を結ぶ線分の長さ OA = OB を求める術を述べよ。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおく。

eq1, eq2 は，3 個の円が相似であることを表している。
eq3, eq4 は，互いに接する直線上の 2 円の間の水平距離に関する式である。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive, r3::positive, x3::positive
eq1 = r3/x3 - r2/x2
eq2 = r2/x2 - r1/x1
eq3 = (x2 - x3) - 2sqrt(r2*r3)
eq4 = (x1 - x2) - 2sqrt(r1*r2);

連立方程式を解き，r2, x1, x2, x3 を一度に求めることもできるが，逐次的に x1 だけを求めてみよう。

相似比を p とすれば，円の半径には以下の関係が成り立つ。
r2/r1 = r3/r2 = p
r2 = p*r1, r3 = p*r2 より r3 = p^2 * r1
甲円と丙円の半径がわかれば，相似比は「丙円と甲円の比の平方根」である。
p = √(r3/r1) 

直線との接点の座標にも
x2/x1 = r2/r1 = p
(x1 - 2√(r1*r2))/x1 = p
の関係があるので，甲円と直線の接点を結ぶ線分の長さは
x1 = 2r1*√p/(1 - p)
である。

甲円と丙円の直径が 20 寸，5寸のとき，2 直線の交点と甲円と直線の接点を結ぶ線分の長さは，56.5685424949238 寸である。

p = sqrt(r3/r1)
ans_x1 = 2r1*√p/(1 - p) |> simplify
ans_x1 |> println

   2*r1^(5/4)*r3^(1/4)/(sqrt(r1) - sqrt(r3))

ans_x1(r1 => 20, r3 => 5).evalf() |>  println

   56.5685424949238

---

連立方程式を解くと，x1 の式はかなり長く，SymPy では直接は簡約化することができないが，以下のように式変換すると上述の式と同じものが得られる。

res = solve([eq1, eq2, eq3, eq4], (r2, x1, x2, x3))[4];  # 4 of 4

res[1](r1 => 20, r3 => 5).evalf() |> println
res[2](r1 => 20, r3 => 5).evalf() |> println
res[3](r1 => 20, r3 => 5).evalf() |> println
res[4](r1 => 20, r3 => 5).evalf() |> println

   10.0000000000000
   56.5685424949238
   28.2842712474619
   14.1421356237310

res[2] |> println

   2*r1*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/r3

ルートの中の式を再定義する（分数式にするため / を // に変換）。

t = (r1^(7//2)*r3^(5//2) - r1^(5//2)*r3^(7//2) - r1^(3//2)*r3^(9//2) + sqrt(r1)*r3^(11//2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)
t |> println

   (r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)

√r1 を u, √r3 を v に置換し，因数分解する。

@syms u, v
t2 = t(√r1 => u, √r3 => v) |> factor
t2 |> println

   u*v^5/(u - v)^2

元の式に戻す。

t3 = (2r1*sqrt(t2)/r3)(u => √r1, v => √r3)
t3 |> println

   2*r1^(5/4)*r3^(1/4)/Abs(sqrt(r1) - sqrt(r3))

t3(r1 => 20, r3 => 5).evalf()　|> println

   56.5685424949238

using LaTeXStrings
function draw(r1, r3, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(700, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   (r2, x1, x2, x3) = (r1*(2*r1^(7/2)*r3^(5/2) - 2*r1^(5/2)*r3^(7/2) - 2*r1^(3/2)*r3^(9/2) + 2*sqrt(r1)*r3^(11/2) + 4*r1^3*r3^3 - 8*r1^2*r3^4 + 4*r1*r3^5 + r3*(r3 - sqrt((4*r1^(7/2)*r3^(5/2) - 4*r1^(5/2)*r3^(7/2) - 4*r1^(3/2)*r3^(9/2) + 4*sqrt(r1)*r3^(11/2) + 8*r1^3*r3^3 - 16*r1^2*r3^4 + 8*r1*r3^5 + r3^2*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(2*(r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)), 2*r1*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/r3, r1*(2*r1^(7/2)*r3^(5/2) - 2*r1^(5/2)*r3^(7/2) - 2*r1^(3/2)*r3^(9/2) + 2*sqrt(r1)*r3^(11/2) + 4*r1^3*r3^3 - 8*r1^2*r3^4 + 4*r1*r3^5 + r3*(r3 - sqrt((4*r1^(7/2)*r3^(5/2) - 4*r1^(5/2)*r3^(7/2) - 4*r1^(3/2)*r3^(9/2) + 4*sqrt(r1)*r3^(11/2) + 8*r1^3*r3^3 - 16*r1^2*r3^4 + 8*r1*r3^5 + r3^2*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r3*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)), 2*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))
   p = r2/r1
   println(2r1*√p/(1 - p))
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, x1, r2, x2, r3, x3)
   p1 = plot()  # showaxis=false)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, r1, L"r_1,(x_1,r_1)", :red, :center, delta=-delta)
       point(x2, r2, L"r_2,(x_2,r_2)", :blue, :center, delta=-delta)
       point(x3, r3, L"r_3,(x_3,r_3)", :green, :center, delta=-delta)
       θ = 2atand(r1/x1)
       abline(0, 0, tand(θ), 0, x1)
       segment(0, 0, 1.2x1, 0)
       point(x1 - r1*sind(θ), r1 + r1*cosd(θ), "B", :black, :right, :bottom, delta=delta)
       point(x2 - r2*sind(θ), r2 + r2*cosd(θ))
       point(x3 - r3*sind(θ), r3 + r3*cosd(θ))
       point(x1, 0, "A", :black, :center, delta=-delta)
       point(x2, 0)
       point(x3, 0)
       point(0, 0, "O", :black, :center, delta=-delta)
       ylims!(-7delta, 2r1 + delta)
   end
end;

draw(20, 5, true)