https://sangaku0418.hatenablog.com/entry/2023/11/22/180643

高山忠直編: 算法評論

国立国会図書館デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8
キーワード：円4個，二等辺三角形
#Julia #SymPy #算額 #和算 #数学

「算額（その164）」と基本的に同じであるが，乙円，丙円の半径を未知数として一般解を求める。

二等辺三角形内に甲円，乙円，丙円を容れる。乙円と丙円の直径が与えられたとき，丙円の直径を求めよ。

一般性を失わずに，丙円の直径を 1 とする。乙円と丙円の直径が与えられたとき，乙円の直径は「乙円の直径/丙円の直径」，丙円の直径は 1 として与える。結果は「丙円の直径/乙円の直径」倍になる。

二等辺三角形の底辺の長さを \(2a\)，高さを \(b とする\)
甲円の半径と中心座標を \(r_1,\ (0,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (0,\ 2r_1 + r_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ r_3)\)
として，以下の連立方程式を解く。

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

using SymPy
@syms a, b, r1, r2, r3, x3
r3 = 1
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = r3/(a - x3) - r1/a
eq3 = r2/(b - 2r1 - r2) - r1/(b - r1)
eq4 = r1/(b - r1) - a/sqrt(a^2 + b^2)
res = solve([eq1, eq2, eq3, eq4], (a, b, r1, x3))[4] # 4 of 4

( (2*sqrt(r2) + 1)^(3/2)/sqrt(r2), 2*(2*sqrt(r2) + 1)^2/(2*sqrt(r2) - r2 + 1), 2*sqrt(r2) + 1, 2*sqrt(2*sqrt(r2) + 1))

甲円の直径は， \(4r_3(\sqrt{r_2/r_3} + 1/2)\) である。

(r2, r3) = (1, 2)
t = r2/r3
u = 2sqrt(t) + 1

\(r_2 = 1;\ r_3 = 2;\ a = 10.6098;\ b = 12.1793;\ r_1 = 4.82843;\ x_3 = 6.2151\)
\(甲円の直径 = 9.65685\)

乙円，丙円の半径を 1, 2 としたとき，甲円の直径は 9.65685 である。

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

function draw(r2, r3, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
t = r2/r3
u = 2sqrt(t) + 1

a = r3*(u^(3/2)/sqrt(t))
b = r3*(2u^2/(u - t))
r1 = r3*u
x3 = r3*(2sqrt(u))
@printf("r2 = %g; r3 = %g; a = %g; b = %g; r1 = %g; x3 = %g\n", r2, r3, a, b, r1, x3)
@printf("甲円の直径 = %g\n", 4r3*(sqrt(r2/r3) + 1/2))
plot([a, 0, -a, 0], [0, b, 0, 0,], color=:black, lw=0.5)
circle(0, r1, r1)
circle(0, 2r1 + r2, r2, :blue)
circle(x3, r3, r3, :green)
circle(-x3, 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(0, r1, " 甲円:r1,(0,r1)", :magenta, :left, :vcenter)
point(0, 2r1 + r2, " 乙円:r2,(0,2r1+r2)", :magenta, :left, :vcenter)
point(x3, r3, " 丙円:r3,(x3,r3)", :magenta, :center, :top, delta=-delta)
end
end;