https://sangaku0418.hatenablog.com/entry/2023/01/25/160421

山形県村山市楯岡最上徳内記念館天明4年(1784)

『続神壁算法』(藤田嘉言)より平成5年(1993)復元
和算の館
http://wasan.jp/yamagata/mogami.html

東京都府中市宮町大國魂神社明治18年(1885)

永井信一: 大國魂神社の算額（下段第5問）についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok1.pdf

キーワード：円4個，外円，三角形
#Julia #SymPy #算額 #和算 #数学　

外円内に三角形，甲円，乙円，丙円が入っている。甲円，乙円，丙円の直径をそれぞれ 49 寸，28 寸，17 寸のとき，外円の直径を求めよ。

三角形の頂点 \({\rm A},\ {\rm B},\ {\rm C}\) の座標をそれぞれ \( (-x_0,\ 2r_1 - R),\ (x_0,\ 2r_1 - R),\ (x_4,\ y_4)\)
甲円の半径と中心座標を \(R,\ (0,\ 49 -R)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ 2r_1 - R + r_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
とおき，以下の連立方程式の数値解を求める。

なお，与えられる条件は，甲円，乙円，丙円の直径の 3 個では多すぎる。

eq7 は，「乙円の中心と外円の中心を結ぶ線分が三角形の辺と直交する」という条件である。すなわち，丙円が最大となる条件である。この条件がないと，丙円と三角形の辺が交差したり，接しないということになる。
この図を描くのに必要なパラメータは \(r_1,\ r_2,\ R,\ r_3,\ x_3,\ y_3,\ x_0,\ x_4, y_4\) の 9 個である。未知数の個数は条件式の個数と同じでなければならない。条件式が 7 個であれば，未知数も 7 個なので，既知数の個数は 2 個であるべきである。乙円が既知であるというのは，過剰な条件である。与えられるべきものは甲円，乙円の直径だけで，丙円の円直径は方程式から決定されるのである。
下図は丙円の拡大図で，算額の通り丙円の直径が 17 寸のときを赤で描き，丙円の直径は計算結果によったものを黒で描いている。

丙円の直径は17寸より小さくなければならない。計算結果による丙円の直径は 16.88961421755243 寸である。
和算家，および，和算解を求めるのを良しとする近現代の和算家（？）は，図のすべてのパラメータを求めて，精密に図を描くということをしないので，未知数と条件数の不一致や，計算結果の不一致に気づかなかったのであろう。

include("julia-source.txt"); # julia-source.txt ソース
function driver(r1, r2)
function H(u)
function parameters()
y2 = 2r1 - R + r2
eq1 = x3^2 + y3^2 - (R - r3)^2
eq2 = R^2 - (2r1 - R)^2 - x0^2
eq3 = x4^2 + y4^2 - R^2;
eq4 = dist3(-x0, 2r1 - R, x4, y4, r2, y2, r2)
eq5 = dist3( x0, 2r1 - R, x4, y4, r2, y2, r2)
eq6 = dist3(-x0, 2r1 - R, x4, y4, x3, y3, r3)
eq7 = y4/(x4 + x0) + x3/y3
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7]
end;
(R, r3, x3, y3, x0, x4, y4) = u
return parameters()
end;
iniv = BigFloat[42.5, -16.0, 30.0, 42.0, 18.0, 38.5].*(r1/24.5)
iniv = BigFloat[42, 17, -15.0, 31.0, 41.0, 19.0, 39].*(r1/24.5)
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
res = driver(49/2, 28/2)
res |> println

[42.5, 8.444807108776216, -18.39145385156791, 28.662006002443494, 42.0, 18.0, 38.5]

丙円の真の直径は 16.88961421755243 であるべきだ。

res[2]*2

16.88961421755243

甲円，乙円（，丙円）の直径がそれぞれ 49 寸，28 寸（，17 寸）のとき，外円の直径は 85 寸である。
外円の直径は 85.0
丙円の直径は 16.88961421755243

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

function draw(r1, r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAMincho")
(R, r3, x3, y3, x0, x4, y4) = driver(r1, r2)
println("外円の直径は ", 2R)
println("丙円の直径は ", 2r3)
plot()
circle(0, 0, R)
circle(0, r1-R, r1, :blue)
circle(r2, 2r1-R+r2, r2, :magenta)
circle(x3, y3, r3, :green)
plot!([x0, x4, -x0, x0], [2r1 - R, y4, 2r1 - R, 2r1 - R], color=:brown, lw=0.5)
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)
@printf("x0 = %.3f, x3 = %.3f, y3 = %.3f, x4 = %.3f, y4 = %.3f", x0, x3, y3, x4, y4)
point(0, r1 - R, " 甲円:r1,(0,r1-R)", :blue, :center, delta=-delta)
point(r2, 2r1 - R + r2, "乙円:r2,(r2,2r1-R+r2)", :magenta, :center, delta=-delta)
point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
point(0, 0, "外円:R,(0,0)", :red, :center, delta=-delta)
point(-x0, 2r1 - R, " A:(-x0,2r1-R)", :brown)
point(x0, 2r1 - R, "B:(x0,2r1-R) ", :brown, :right)
point(x4, y4, " C:(x4,y4)", :brown, :left, :bottom)
end
end;
draw(49/2, 28/2, true)