裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

和算の心(その006)

2024年07月27日 | Julia

和算の心(その006)

直線上の 3 個の円

直線上に,甲円,乙円,丙円が互いに接しており,この順序で載っている。甲円,乙円と直線の接点間の距離(以下では,水平距離と呼ぶ)を求めよ。

甲円,乙円,丙円の半径を r1,r2,r3 とする。

1. 一般的な場合

次節で述べる特別な場合,すなわち 3 円が互いに接している(甲円と乙円も接している)場合よりも甲円と丙円が離れている場合である。

甲円と乙円,乙円と丙円の水平距離は である。(和算の心(その002)

結果として,甲円と丙円の水平距離は になる。

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

include("julia-source.txt");

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

SymPy では,これ以上は自動的に簡約化されないが,2√r2(√r1 + √r3) である。

res[2][2] |> simplify

   

res[2][2](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

using LaTeXStrings
function draw(r1=3, r2=2, r3=1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   x2 = 2√(r1*r2)
   x3 = 2√r2*(√r1 + √r3)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3+r3, 0, :black)
   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(0, r1, L"r_1,(0,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)
       point(x3 - r3, 1.3*2r3, L"x_3 = 2\sqrt{r_2}(\sqrt{r_1} + \sqrt{r_3})", mark=false)
       point(x3 - r3, 1.15*2r3, L"x_2 = 2\sqrt{r_1 r_2})", mark=false)
   end
   display(p1)
end;

draw(3, 1, 2, true)

   r1 = 3;  r2 = 1;  x2 = 3.4641;  r3 = 2;  x3 = 6.29253


2.  3 個の円が互いに接する場合

この状況では,乙円の半径には制約が付く(どのような値でも取れるわけではない)。そのため,甲円と乙円,乙円と丙円の水平距離は簡単な数式では表せない。しかし,甲円と丙円の水平距離は と簡単なままである。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおき以下の連立方程式を解く。

include("julia-source.txt");

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

SymPy でこの連立方程式を解くと,r2 と x2 は複雑な式になる。

res[1][1]  # r2

   

√r1 => r1, √r3 => r3 と変換したうえで因数分解し,r1 => √r1, r3 => √r3 でもとに戻すと簡約化された式になる。

res[1][1](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

res[1][2]  # x2

   

res[1][2](√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

x3 は簡潔な式のままである。

res[1][3]  #x3

   

別の関係式を使って連立方程式を解くと,特別なことをしないでも簡潔な式が得られる。

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive
eq1 = x2 - 2√(r1*r2)
eq2 = x3 - 2√(r1*r3)
eq3 = (x3 - x2) - 2√(r2*r3)
res = solve([eq1, eq2, eq3], (r2, x2, x3));

res[1][1] |> display
res[1][2] |> display
res[1][3] |> display


   

   

   

デカルトの円定理による場合

デカルトの円定理は,半径の関係のみである。半径を求めた後で,位置を求めればよい。

@syms r1::positive, r2::positive, r3::positive
r2 = 1/(1/r1 + 1/r3 + 2√(1/(r1*r3)))

   

これも,√r1 => r1, √r3 => r3 と変換したうえで因数分解し,r1 => √r1, r3 => √r3 でもとに戻すと簡約化された式になる。

r2(√r1 => r1, √r3 => r3) |> factor |> (x -> x(r1 => √r1, r3 => √r3))

   

using LaTeXStrings
function draw(r1=3, r3=1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   (r2, x2, x3) = (r1*r3/(sqrt(r1) + sqrt(r3))^2, 2*r1*sqrt(r3)/(sqrt(r1) + sqrt(r3)), 2*sqrt(r1)*sqrt(r3))
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3+r3, 0, :black)
   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(0, r1, L"r_1,(0,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)
       point(x3 - r3, 1.45*2r3, L"r_2 = \frac{r_1 r_3}{(\sqrt{r1} + \sqrt{r3})^2}", mark=false)
       point(x3 - r3, 1.15*2r3, L"x_3 = 2\sqrt{r_1 r_3}", mark=false)
   end
   display(p1)
end;

draw(3, 2, true)

   r1 = 3;  r2 = 0.606123;  x2 = 2.69694;  r3 = 2;  x3 = 4.89898

 

3. 3 個の円が互いに接し,両端の円の半径が等しい場合

乙円,丙円の中心座標を (x2, r2), (x3, r3) とするが,2√(r1*r2) + 2√(r2*r3) = 2√(r1*r3) となる。
x2 = r1, x3 = 2r1 は自明である。
さらに,r1 = r3 である。
以下の方程式を解くと r2 = r1/4 であることがわかる。

include("julia-source.txt");

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive
x2 = r1
eq1 = r1^2 + (x2 - r2)^2 - (r1 + r2)^2
solve(eq1, r2)[1]  # r2

   

または,前節で得られた r2 の式を r3 = r1 として簡約化すればよい。

r3 = r1
r1*r3/(sqrt(r1) + sqrt(r3))^2  # r2

   

または,もう少し前の方程式をr3 = r1 として解いてもよい。

r3 = r1
eq = 2√(r1*r2) + 2√(r2*r3) - 2√(r1*r3) 
solve(eq, r2)[1]  # r2

   

using LaTeXStrings
function draw(r1, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   r3 = r1
   x3 = 2r1
   r2 = r1/4
   x2 = r1
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   p1 = plot(ylims=(-0.05r1, 1.05*2r1), showaxis=false)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   segment(-r1, 0, x3 + r3, 0, :black)
   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(0, r1, L"r_1,(0,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)
       point(x3 - r3, 1.15*2r3, L"$r_3 = r_1,\ x_3 = 2r_1$", mark=false)
       point(x3 - r3, 1.3*2r3, L"$r_2 = \frac{r_1}{4},\ x_2 = r_1$", mark=false)
   end
   display(p1)
end;

draw(3, true)

   r1 = 3;  r2 = 0.75;  x2 = 3;  r3 = 3;  x3 = 6

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村