https://sangaku0418.hatenablog.com/entry/2025/03/08/153229

45 岩手県一関市舞川字龍ヶ沢観福寺観音堂明治34年(1901)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円2個，楕円，正三角形
#Julia #SymPy #算額 #和算 #数学

楕円の中に曲率円を 2 個設け，それらの円と外接するように正三角形を容れる。楕円の長径と短径が与えられたとき，正三角形の一辺の長さを求める術を述べよ。

楕円の長半径，短半径，中心座標を \(a, b, (0, 0)\)
曲率円の半径と中心座標を \(r, (a - r, 0), (r - a, 0)\)
正三角形の一辺の長さと頂点の座標を \(c, (2r - a, y_1), (2r - a, y_1 + c), (x_0, y_0)\)
とおき，以下の連立方程式を解く。

正三角形の頂点 \( (x_0, y_0)\) は楕円の周上にある。

SymPy の能力では解析解を求めることができないので，特定の \(a, b\) における数値解を求める。

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

using SymPy

function getc(a, b)
@syms r::positive,
x0::positive, y0::positive, y1::negative,
c::positive;
eq1 = r - b^2/a
eq2 = x0^2/a^2 + y0^2/b^2 - 1
eq3 = dist(2r - a, y1, x0, y0, a - r, 0) - r^2
eq4 = √Sym(3)*(y0 - y1) -(x0 - 2r + a)
eq5 = √Sym(3)*(y0 - y1 - c) + (x0 - 2r + a)
res = solve([eq1, eq2, eq3, eq4, eq5], (r, x0, y0, y1, c))[1]
@printf("楕円の長径，短径が %g，%g のとき，正三角形の一辺の長さは %.15g である。\n", 2a, 2b, res[5])
return res
end;

楕円の長径，短径が 18，10 のとき，正三角形の一辺の長さは 11.5725008614986 である。

「術」の解説は以下のようなっており，11.447328388628943 という，微妙に違う解を与える。

function jutu(長径, 短径)
天 = 長径^2 - 短径^2
地 = sqrt(天 - 短径^2)
人 = sqrt(天)
A = √3長径
B = 短径^2/A
C = 地 + 人 - B
println(C/2)
end;
jutu(18, 10)

11.447328388628943

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

function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, x0, y0, y1, c) = float.(getc(a, b))
plot([x0, 2r - a, 2r - a, x0], [y0, y1 + c, y1, y0], color=:green, lw = 0.5)
circle2(a - r, 0, r)
ellipse(0, 0, a, b, color=:blue)
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(0, b, "b", :blue, :center, :bottom, delta=delta/2, deltax=delta/2)
point(a, 0, "a", :blue, :left, :bottom, delta=delta/2, deltax=delta/2)
point(a - r, 0, "a-r", :red, :center, :bottom, delta=delta/2)
point(2r - a, 0, "2r-a", :blue, :right, :bottom, delta=delta/2, deltax=-delta/2)
point(2r - a, y1, "(2r-a,y1)", :green, :center, delta=-delta/2, deltax=-delta/2)
point(2r - a, y1 + c, "(2r-a,y1+c)", :green, :right, delta=-delta/2, deltax=-delta/2)
point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
end
end;