https://sangaku0418.hatenablog.com/entry/2024/12/05/001739

街角の数学 Street Wasan　～落書き帳「○△□」～ 934. 『算法天生指南』巻之二（その９）

http://streetwasan.web.fc2.com/math20.04.29.3.html
キーワード：円4個，直線上
#Julia #SymPy #算額 #和算 #数学

問題 III. 直線上に 4 個の円が載っている。甲円，乙円，丁円の直径が 72 寸，24 寸，32 寸のとき，丙円の直径はいかほどか。

甲円の半径と中心座標を \(r_1,\ (x_1,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ r_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
丁円の半径と中心座標を \(r_4,\ (0,\ r_4)\)
とおき，以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive,
r3::positive, x3::positive, y3::positive, r4::positive
# x2 = 2sqrt(r2*r4)
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = x2^2 + (r4 - r2)^2 - (r2 + r4)^2
eq5 = x3^2 + (y3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (x1, x2, r3, x3, y3))[3]; # 3 of 4

それぞれは，長い式になっているので，以下のようにして簡約化する。

1. 「\(\texttt{変数^(n/2)}\)」は，\(\texttt{r1^(1//2) => s, r4^(1//2) => t}\) のように，変数変換を行う。（このようにしないと因数分解できない）
2. 最後に \(\texttt{s => r1^(1//2), t => r4^(1//2)}\) と変数変換し，もとに戻す。

1.　\(x_1\)

変数変換，因数分解，簡約化を経て，以下のようになる。

@syms s, t
ans_x1 = res[1](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は，式を二乗して因数分解すると二乗の分数式になるので，平方根をとって \(\frac{r_2 s + r_2 t - 2s t^2}{r_2 - s t}\) である。

ans_x1.args[3]^2 |> factor |> println

(r2*s + r2*t - 2*s*t^2)^2/(r2 - s*t)^2

根号の外と，分母に同じ項を含むので約分して以下のように簡単になる。

ans_x1 = 2sqrt(r2)* (r2*s + r2*t - 2s*t^2)/(r2 - s*t) *(r2 - s*t)*(s + t)/(r2*s + r2*t- 2s*t^2)
ans_x1 |> println

2*sqrt(r2)*(s + t)

必要ならば，変数変換をもとに戻す。

ans_x1 = ans_x1(s => r1^(1//2), t => r4^(1//2))
ans_x1 |> println

2*sqrt(r2)*(sqrt(r1) + sqrt(r4))

簡約化の前後で同値な式であることを確認する。

res[1](r1 => 72/2, r2 => 24/2, r4 => 32/2) |> println
ans_x1(r1 => 72/2, r2 => 24/2, r4 => 32/2) |> println

69.2820323027551
69.2820323027551

2.　\(x_2\)

@syms s, t
ans_x2 = res[2](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は \(x_1\) と同じ式 \(\frac{r_2 s + r_2 t - 2s t^2}{r_2 - s t}\) であるから，以下のように簡約化される。

ans_x2 = 2sqrt(r2)*t* (r2*s + r2*t - 2s*t^2)/(r2 - s*t) *(r2 - s*t)/(r2*s + r2*t- 2s*t^2)
ans_x2 |> println

2*sqrt(r2)*t

ans_x2 = ans_x2(s => r1^(1//2), t => r4^(1//2))
ans_x2 |> println

2*sqrt(r2)*sqrt(r4)

res[2](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_x2(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

27.7128129211020
27.7128129211020

3.　\(r_3\)

@syms s, t
ans_r3 = res[3](r1^(1//2) => s, r4^(1//2) => t) |> factor
ans_r3 |> println

-r2^2*(s + t)^2/(4*s*t*(r2 - s*t))

ans_r3 = ans_r3(s => r1^(1//2), t => r4^(1//2))
ans_r3 |> println

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

res[3](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_r3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

12.5000000000000
12.5000000000000

4.　\(x_3\)

@syms s, t
ans_x3 = res[4](r1^(1//2) => s, r4^(1//2) => t) |> factor

長い式の平方根は \(x_1\) と同じ式であるから，以下のように簡約化される。

ans_x3 = sqrt(r2)* (r2*s + r2*t - 2s*t^2)/(r2 - s*t)
ans_x3 |> println

sqrt(r2)*(r2*s + r2*t - 2*s*t^2)/(r2 - s*t)

ans_x3 = ans_x3(s => r1^(1//2), t => r4^(1//2))
ans_x3 |> println

sqrt(r2)*(sqrt(r1)*r2 - 2*sqrt(r1)*r4 + r2*sqrt(r4))/(-sqrt(r1)*sqrt(r4) + r2)

res[4](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_x3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

20.7846096908265
20.7846096908265

5.　\(y_3\)

@syms s, t
ans_y3 = res[5](r1^(1//2) => s, r4^(1//2) => t) |> factor
ans_y3 |> println

r2*(r2*s^2 + 2*r2*s*t + r2*t^2 - 8*s^2*t^2)/(4*s*t*(r2 - s*t))

ans_y3 = ans_y3(s => r1^(1//2), t => r4^(1//2))
ans_y3 |> println

r2*(2*sqrt(r1)*r2*sqrt(r4) + r1*r2 - 8*r1*r4 + r2*r4)/(4*sqrt(r1)*sqrt(r4)*(-sqrt(r1)*sqrt(r4) + r2))

res[5](r1 => 36, r2 => 12, r4 => 16).evalf() |> println
ans_y3(r1 => 36, r2 => 12, r4 => 16).evalf() |> println

35.5000000000000
35.5000000000000

一時変数を設定すると若干短く定義できる。
s = sqrt(r1) + sqrt(r4)
t = (r2 - sqrt(r1*r4))
x1 = 2sqrt(r2)*s
x2 = 2sqrt(r2*r4)
r3 = -r2^2*s^2/(4sqrt(r1*r4)*t)
x3 = sqrt(r2)*(sqrt(r1)*r2 - 2sqrt(r1)*r4 + r2*sqrt(r4))/t
y3 = r2*(2sqrt(r1*r4)*r2 + r1*r2 - 8r1*r4 + r2*r4)/(4sqrt(r1*r4)*t)

描画関数プログラムのソースを見る

function draw(r1, r2, r4, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
s = sqrt(r1) + sqrt(r4)
t = (r2 - sqrt(r1*r4))
x1 = 2sqrt(r2)*s
x2 = 2sqrt(r2*r4)
r3 = -r2^2*s^2/(4sqrt(r1*r4)*t)
x3 = sqrt(r2)*(sqrt(r1)*r2 - 2sqrt(r1)*r4 + r2*sqrt(r4))/t
y3 = r2*(2sqrt(r1*r4)*r2 + r1*r2 - 8r1*r4 + r2*r4)/(4sqrt(r1*r4)*t)
@printf("r1 = %g; r2 = %g; r4 = %g; r3 = %g\n", r1, r2, r4, r3)
@printf("x1 = %g; x2 = %g; r3 = %g; x3 = %g; y3 = %g\n", x1, x2, r3, x3, y3)
plot()
circle(x1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, y3, r3, :green)
circle(0, r4, r4, :magenta)
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, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
point(x2, r2, "乙円:r2\n(x2,r2)", :blue, :center, delta=-delta)
point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
point(0, r4, "丁円:r4,(0,r4)", :magenta, :center, delta=-delta)
end
end;