https://sangaku0418.hatenablog.com/entry/2023/07/25/175608

岐阜県大垣市赤坂町金生山明星輪寺元治2年(1865)

http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF
キーワード：円10個，半円，長方形
#Julia #SymPy #算額 #和算 #数学

第四問　長方形内に青の半円 2 個，赤円 1 個，緑円 2 個，浅青円 1 個，白円 2 個，および薄紅円 2 個，無名円 2 個を容れる。緑円は青，赤，浅青のいずれにも外接し，緑，浅青円は長方形の辺にも接している。

注：算額（少なくとも水野の図）では，緑円と赤円，浅青円と赤円が接していないような，誤解を与えかねない図になっている。

長方形の左下隅を原点とし，以下のように記号を定める。
青円の半径と中心座標 \(r_1,\ (r_1,\ 0)\)
赤円の半径と中心座標 \(r_2,\ (2r_1 - r_2,\ y_2)\)
緑円の半径と中心座標 \(r_3,\ (r_3,\ y_3)\)
浅青円の半径と中心座標 \(r_5,\ (r_5,\ y_2)\)
白円の半径と中心座標 \(r_4,\ (2r_1 - r_4,\ y_4)\)
薄紅円の半径と中心座標 \(r_6,\ (x_6,\ y_6)\)
無名円の半径と中心座標 \(r_7,\ (x_7,\ y_7)\)
長方形の幅は \(2r_1，高さは 2y_2 である。\)

\(r_1\) を既知とすれば，残りの未知数は 13 変数なので，13 元連立方程式になるが，solve() では一度には解けないので，薄紅円と無名円を除いてまず解を求め，しかる後に得られた解を既知として，薄紅円と無名円のパラメータを求める。

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

using SymPy

@syms r1::positive, r2::positive, y2::positive, r3::positive, y3::positive,
r4::positive, y4::positive, r5::positive, r6::positive, x6::positive,
y6::positive, r7::positive, x7::positive, y7::positive;

eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = (r1 - r4)^2 + y4^2 - (r1 + r4)^2
eq5 = (2r1 - r2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (r2 - r4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq9 = (r3 - r5)^2 + (y2 - y3)^2 - (r3 + r5)^2
eq13 = 2r5 + 2r2 - 2r1
res1 = solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13], (r2, y2, r3, y3, r4, y4, r5))
solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13, eq4, eq7, eq8, eq10, eq11, eq12], (r2, y2, r3, y3, r4, y4, r5, r6, x6, y6, r7, x7, y7))[1]

(3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)

浅青円の半径 \(r_5\) は \(r_1/4\) すなわち，青円の半径 \(r_1\) の \(1/4\) である。

以下は図を描くためだけに必要なものである。

得られた \(r_2,\ y_2,\ r_3,\ y_3,\ r_4,\ y_4,\ r_5\) を既知として，次に，\(r_6,\ x_6,\ y_6,\ r_7,\ x_7,\ y_7\) に関する連立方程式を解く。

using SymPy

(r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(Sym(3))*r1, r1/3, 2*sqrt(Sym(3))*r1/3, 3*r1*(7 - 4*sqrt(Sym(3))), 2*r1*(-3 + 2*sqrt(Sym(3))), r1/4)
eq4 = (x6 - r1)^2 + y6^2 - (r1 + r6)^2
eq7 = (2r1 - r2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq8 = (2r1 - r2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2
eq10 = (x6 - r3)^2 + (y6 - y3)^2 - (r3 + r6)^2
eq11 = (x7 - r3)^2 + (y7 - y3)^2 - (r3 + r7)^2
eq12 = (x7 - r5)^2 + (y2 - y7)^2 - (r5 + r7)^2
res2 = solve([eq4, eq7, eq8, eq10, eq11, eq12], (r6, x6, y6, r7, x7, y7))[1]

(3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)

二段階に分けて求めたパラメータで図を描き，解が正しいことを確認する。

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

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1
(r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)
(r6, x6, y6, r7, x7, y7)= (3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)
plot([0, 2r1, 2r1, 0, 0], [0, 0, 2y2, 2y2, 0], color=:black, lw=0.5)
circlef(r1, 0, r1, :dodgerblue, beginangle=0, endangle=180)
circlef(r1, 2y2, r1, :dodgerblue, beginangle=180, endangle=360)
circlef(2r1 - r2, y2, r2, :red)
circlef(r3, y3, r3, :green)
circlef(r3, 2y2 - y3, r3, :green)
circle(2r1 - r4, y4, r4, :orange)
circle(2r1 - r4, 2y2 - y4, r4, :orange)
circlef(r5, y2, r5, :skyblue1)
circlef(x6, y6, r6, :pink)
circlef(x6, 2y2 - y6, r6, :pink)
circlef(x7, y7, r7, :orange)
circlef(x7, 2y2 - y7, r7, :orange)
if more
point(r1, 0, "r1", :black, :center, :bottom)
point(2r1 - r2, y2, "(2r1-r2,y2)", :black, :center, :bottom)
point(r3, y3, "(r3,y3)", :black, :center, :bottom)
point(2r1 - r4, y4, "(2r1-r4,y4) ", :black, :center, :bottom)
point(r5, y2, "(r5,y2)", :black, :center, :bottom)
point(x6, y6, " (x6,y6)", :black, :left, :bottom)
point(x7, y7, " (x7,y7)", :black, :left, :vcenter)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;