https://sangaku0418.hatenablog.com/entry/2024/02/03/174536

『不朽算法』上巻第六問 [日下貞八郎誠嗣辺・安島万蔵直円遺稿]

米山忠興: 側円術(IV) ---『安子遺稿側円解二条第一』の現代語訳 ---，東洋大学紀要. 自然科学篇，Vol. 52, pp. 97-116, 2008-03
http://id.nii.ac.jp/1060/00002535/
キーワード：円1個，楕円，直角三角形
#Julia #SymPy #算額 #和算 #数学

直角三角形（鈎股弦）内に円と楕円が入っている。円は直角三角形の三辺に接しており，楕円と外接している。楕円の長径は底辺（股）に平行で，直角三角形の斜辺（弦）に接している。
直角三角形の辺の長さと楕円の短径を与えて，長径を求めよ。

直角三角形の直角を挟む二辺のうち，短い方を「鈎」，長い方を「股」とする（斜辺「弦」は \(\sqrt{鈎^2 + 股^2}\) である）。
円の半径 \(r\) は「\( (鈎 + 股 - 弦)/2\)」，中心座標は \( (r,\ r)\)
楕円の長径と短径を \(2a,\ 2b\)，中心座標を \( (x_0,\ y_0)\) とする。
楕円と斜辺，楕円と円の接点座標を \( (x_1,\ y_1),\ (x_2,\ y_2)\)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms 鈎, 股, 弦, a, b, r, x, x1, y1, x2, y2, x0, y0, slope0

(鈎, 股, b) = (2240, 2808, 455//2)
弦 = 3592 # sqrt(鈎^2 + 股^2)
r = (鈎 + 股 - 弦)//2
y0 = b
長径，短径が a b，中心座標が (x0, y0) の楕円
(x1, y1) が楕円上にある
eq1 = (x1 - x0)^2/a^2 + (y1 - y0)^2/b^2 - 1
(x1, y1) における楕円の接線の傾き
slope0 = -鈎//股 # 目的とする接線。符号に注意
eq2 = b^2*(x1 - x0)/(a^2*(y0 - y1)) - slope0 # slope
もう一個の条件
eq3 = y1/(股 - x1) + slope0
楕円と円が外接する
eq4 = (x2 - x0)^2/a^2 + (y2 - y0)^2/b^2 - 1
slope2 = (x2 - r)/(r - y2)
eq5 = b^2*(x2 - x0)/(a^2*(y0 - y2)) - slope2
eq6 = (x2 - r)^2 + (y2 - r)^2 - r^2;

function H(u)
(x1, y1, x0, x2, y2, a) = u
return [
4*(y1 - 455/2)^2/207025 - 1 + (-x0 + x1)^2/a^2, # eq1
280/351 + (-207025*x0/4 + 207025*x1/4)/(a^2*(455/2 - y1)), # eq2
y1/(2808 - x1) - 280/351, # eq3
4*(y2 - 455/2)^2/207025 - 1 + (-x0 + x2)^2/a^2, # eq4
-(x2 - 728)/(728 - y2) + (-207025*x0/4 + 207025*x2/4)/(a^2*(455/2 - y2)), # eq5
(x2 - 728)^2 + (y2 - 728)^2 - 529984, # eq6
]
end;

iniv = BigFloat[2398, 327, 1870, 1300, 230, 580]
res = nls(H, ini=iniv)

([2397.8426966292136, 327.19101123595505, 1872.0, 1310.4, 291.2, 585.0], true)

楕円の長径(a) = 585（差渡し径(2a) = 1170）である。

その他のパラメータは以下の通り。

\(鈎 = 2240; 股 = 2808; b = 227.5\)
\(x_1 = 2397.84; y_1 = 327.191; x_0 = 1872; y_0 = 227.5\)
\(x_2 = 1310.4; y_2 = 291.2; r = 728\)

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

function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(鈎, 股, b) = (2240, 2808, 455//2)
弦 = sqrt(鈎^2 + 股^2)
r = (鈎 + 股 - 弦)/2
y0 = b
(x1, y1, x0, x2, y2, a) = res[1]
@printf("鈎 = %g; 股 = %g; b = %g\n", 鈎, 股, b)
@printf("楕円の長径(a) = %g　（差渡し径(2a) = %g\n", a, 2a)
@printf("x1 = %g; y1 = %g; x0 = %g; y0 = %g\n", x1, y1, x0, y0)
@printf("x2 = %g; y2 = %g; r = %g\n", x2, y2, r)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
ellipse(x0, b, a, b, color=:red)
circle(r, r, r, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(x0, y0, "(x0,y0)")
point(r, r, "(r,r)")
point(0, 鈎, " 鈎", :black, :left, :bottom)
point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
point(x1, y1, "(x1,y1) ", :red, :right)
point(x2, y2, "(x2,y2) ", :red, :right)
end
end;