https://sangaku0418.hatenablog.com/entry/2024/12/10/002108

藤田貞資(1789)：神壁算法巻上

http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
所懸于東都不忍辯才天堂者一事天明8年(1788)
末術關流神谷幸吉定令門人東都本郷二丁目佐久間直次郎善久
キーワード：三角形，矢
#Julia #SymPy #算額 #和算 #数学

三角形内を 3 本の線分で区画して，元の三角形と相似な小さい三角形を作る。対応する頂点を結ぶ線分を，甲矢，乙矢，丙矢と呼ぶ。

元の三角形の 3 辺（大斜，中斜，小斜）が 6 寸，4 寸，3 寸で，甲矢，乙矢，丙矢の和が 5 寸のとき，甲矢の長さは何寸か。

図を構成する座標点を，(0, 0), (x1, 0), (x2, y2), (x3, y3), (x4, y4), (x5, y5) とおき，以下の連立方程式を解く。

解析解を求めるのは Julia の SymPy.solve では能力的に無理なので，　nlsolve により数値解を求める。

甲矢は 2 寸である。ちなみに，乙矢は 2.25 寸，丙矢は 0.75 寸である。
全てのパラメータは以下のとおりである。
x1 = 6; x2 = 2.41667; y2 = 1.77756; x3 = 2.06342; y3 = 0.897116; x4 = 4.05; y4 = 0.44439; x5 = 2.99769; y5 = 1.30331

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

using SymPy

@syms x1, x2, y2, x3, y3, x4, y4, x5, y5,
大, 中, 小, 大2, 中2, 小2, 甲矢, 乙矢, 丙矢, 計
(大, 中, 小, 計) = (6, 4, 3, 5)
x1 = 大
大2 = sqrt( (x4 - x3)^2 + (y3 - y4)^2)
中2 = sqrt( (x4 - x5)^2 + (y5 - y4)^2)
小2 = sqrt( (x5 - x3)^2 + (y5 - y3)^2)
甲矢 = sqrt( (x1 - x4)^2 + y4^2)
乙矢 = sqrt(x3^2 + y3^2)
丙矢 = sqrt( (x5 - x2)^2 + (y2 - y5)^2)
eq1 = 中^2 - ( (x1 - x2)^2 + y2^2)
eq2 = 小^2 - (x2^2 + y2^2)
eq3 = 小2/小 - 中2/中
eq4 = 小2/小 - 大2/大
eq5 = 甲矢 + 乙矢 + 丙矢 - 計
eq6 = (甲矢 + 大2)^2 - ( (x1 - x3)^2 + y3^2)
eq7 = (丙矢 + 中2)^2 - ( (x4 - x2)^2 + (y2 - y4)^2)
eq8 = (乙矢 + 小2)^2 - (x5^2 + y5^2);

# 連立方程式の解は以下で求められるべきであるが，SymPy の能力的に不可能である
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x2, y2, x3, y3, x4, y4, x5, y5))

# x2, y2 については，以下で求めることができる
res = solve([eq1, eq2], (x2, y2))[2]

(29/12, sqrt(455)/12)

function H(u)
(x2, y2, x3, y3, x4, y4, x5, y5) = u
return [
-y2^2 - (6 - x2)^2 + 16, # eq1
-x2^2 - y2^2 + 9, # eq2
sqrt( (-x3 + x5)^2 + (-y3 + y5)^2)/3 - sqrt( (x4 - x5)^2 + (-y4 + y5)^2)/4, # eq3
-sqrt( (-x3 + x4)^2 + (y3 - y4)^2)/6 + sqrt( (-x3 + x5)^2 + (-y3 + y5)^2)/3, # eq4
sqrt(x3^2 + y3^2) + sqrt(y4^2 + (6 - x4)^2) + sqrt( (-x2 + x5)^2 + (y2 - y5)^2) - 5, # eq5
-y3^2 - (6 - x3)^2 + (sqrt(y4^2 + (6 - x4)^2) + sqrt( (-x3 + x4)^2 + (y3 - y4)^2))^2, # eq6
-(-x2 + x4)^2 - (y2 - y4)^2 + (sqrt( (-x2 + x5)^2 + (y2 - y5)^2) + sqrt( (x4 - x5)^2 + (-y4 + y5)^2))^2, # eq7
-x5^2 - y5^2 + (sqrt(x3^2 + y3^2) + sqrt( (-x3 + x5)^2 + (-y3 + y5)^2))^2, # eq8
]
end;
x1 = 6
iniv = BigFloat[30, 23, 21, 10, 55, 3, 37, 17] .* 6/74.6
res = nls(H, ini=iniv)

([2.4166666666666665, 1.7775607506417952, 2.0634154925612784, 0.8971156586851534, 4.04999554843813, 0.4443901876604488, 2.9976902975298994, 1.3033123554115043], true)

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

function draw(j, k, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
x1 = 6
(x2, y2, x3, y3, x4, y4, x5, y5) = [30, 23, 21, 10, 55, 3, 37, 17] .* 6/74.6
(x2, y2, x3, y3, x4, y4, x5, y5) = res[1]
甲矢 = sqrt( (x1 - x4)^2 + y4^2)
乙矢 = sqrt(x3^2 + y3^2)
丙矢 = sqrt( (x5 - x2)^2 + (y2 - y5)^2)
@printf("甲矢 = %g; 乙矢 = %g; 丙矢 = %g\n", 甲矢, 乙矢, 丙矢)
@printf("x1 = %g; x2 = %g; y2 = %g; x3 = %g; y3 = %g; x4 = %g; y4 = %g; x5 = %g; y5 = %g\n", x1, x2, y2, x3, y3, x4, y4, x5, y5)
plot([0, x1, x2, 0], [0, 0, y2, 0], color=:blue, lw=0.5)
segment(0, 0, x5, y5)
segment(x3, y3, x1, 0)
segment(x2, y2, x4, y4)
segment(x1, 0, x4, y4, :green, lw=1.5)
segment(0, 0, x3, y3, :green, lw=1.5)
segment(x2, y2, x5, y5, :green, lw=1.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)
point(0, 0, "0", :red, :center, delta=-3delta)
point(x1, 0, "x1", :red, :center, delta=-3delta)
point(x2, y2, "(x2,y2)", :red, :center, :bottom, delta=3delta)
point(x3, y3, "(x3,y3)", :red, :left, delta=-4delta, deltax=-8delta)
point(x4, y4, "(x4,y4)", :red, :left, delta=-3delta, deltax=-15delta)
point(x5, y5, "(x5,y5)", :red, :right, :vcenter, deltax=-5delta)
point(x2/2, y2/2, "小斜 ", :blue, :right, :vcenter, mark=false)
point( (x1 + x2)/2, y2/2, " 中斜", :blue, :left, :vcenter, mark=false)
point(x1/2, 0, "大斜", :blue, :center, delta=-3delta, mark=false)
point( (x4+x1)/2, y4/2, "甲矢", :green, :right, mark=false)
point(x3/2, y3/2, "乙矢", :green, :left, mark=false)
point( (x5 + x2)/2, (y2 + y5)/2, "丙矢 ", :green, :right, :vcenter, mark=false)
ylims!(-13delta, y2 + 3delta)
end
end;