算額（その0191） - 算額あれこれ

長野県中野市田上田上観音堂文化6年(1809)

中村信弥「改訂増補長野県の算額」県内の算額(P.88)
http://www.wasan.jp/zoho/zoho.html
キーワード：円2個，外円，四辺形，対角線
#Julia #SymPy #算額 #和算 #数学

大円内に，天斜，地斜，左斜，右斜で構成される四辺形及びその対角線（乾斜，坤斜）を描く。天斜，乾斜，坤斜で囲まれる領域に小円を描く。
小円と大円の直径の（値の）積が 780 寸（本当は寸^2=歩？），また，天斜，左斜，右斜の長さがそれぞれ 39 寸，60 寸，25 寸のとき，地斜の長さはいかほどか。

四辺形の頂点を \( (x_1,\ y_1),\ (x_2,\ y_2),\ (x_3,\ y_3),\ (x_4,\ y_4)\)，小円の中心座標と半径を \( (x_5,\ y_5),\ r_5\) とおく。そのままでは未知数が 12 個，条件が 11 個で解けないが，図は回転しても一般性を失わないので，\(x_1 = 0\) になるように座標軸を定めれば方程式が解ける。また， \(x_1 = 0\) なら， \(y_1\) は大円の円周上になるので，未知数も方程式も 1 個ずつ減る。

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

using SymPy

@syms x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, r5, R;

x1 = 0
y1 = -R
eq1 = R * r5 - 780//4
eq2 = x2^2 + y2^2 - R^2
eq3 = x3^2 + y3^2 - R^2
eq4 = x4^2 + y4^2 - R^2
eq5 = (x4 - x1)^2 + (y4 - y1)^2 - 39^2
eq6 = (x4 - x3)^2 + (y4 - y3)^2 - 60^2
eq7 = (x1 - x2)^2 + (y1 - y2)^2 - 25^2
eq8 = distance(x4, y4, x1, y1, x5, y5) - r5^2 |> simplify
eq9 = distance(x4, y4, x2, y2, x5, y5) - r5^2 |> simplify
eq10 = distance(x3, y3, x1, y1, x5, y5) - r5^2 |> simplify;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

nlsolve() で解く。

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
println(eq6, ",")
println(eq7, ",")
println(eq8, ",")
println(eq9, ",")
println(eq10, ",")

R*r5 - 195,
-R^2 + x2^2 + y2^2,
-R^2 + x3^2 + y3^2,
-R^2 + x4^2 + y4^2,
x4^2 + (R + y4)^2 - 1521,
(-x3 + x4)^2 + (-y3 + y4)^2 - 3600,
x2^2 + (-R - y2)^2 - 625,
(-r5^2*(R^2 + 2*R*y4 + x4^2 + y4^2)^2 + x4^2*(R*x4 - R*x5 + x4*y5 - x5*y4)^2 + (x4*(R^2 + R*y4 + R*y5 + x4*x5 + y4*y5) - x5*(R^2 + 2*R*y4 + x4^2 + y4^2))^2)/(R^2 + 2*R*y4 + x4^2 + y4^2)^2,
(-r5^2*(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2 + (x2 - x4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2 + (y2 - y4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2)/(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2,
(-R^2*r5^2 + R^2*x3^2 - 2*R^2*x3*x5 + R^2*x5^2 - 2*R*r5^2*y3 + 2*R*x3^2*y5 - 2*R*x3*x5*y3 - 2*R*x3*x5*y5 + 2*R*x5^2*y3 - r5^2*x3^2 - r5^2*y3^2 + x3^2*y5^2 - 2*x3*x5*y3*y5 + x5^2*y3^2)/(R^2 + 2*R*y3 + x3^2 + y3^2),

function H(u)
(x2, y2, x3, y3, x4, y4, x5, y5, r5, R) = u
return [
R*r5 - 195,
-R^2 + x2^2 + y2^2,
-R^2 + x3^2 + y3^2,
-R^2 + x4^2 + y4^2,
x4^2 + (R + y4)^2 - 1521,
(-x3 + x4)^2 + (-y3 + y4)^2 - 3600,
x2^2 + (-R - y2)^2 - 625,
(-r5^2*(R^2 + 2*R*y4 + x4^2 + y4^2)^2 + x4^2*(R*x4 - R*x5 + x4*y5 - x5*y4)^2 + (x4*(R^2 + R*y4 + R*y5 + x4*x5 + y4*y5) - x5*(R^2 + 2*R*y4 + x4^2 + y4^2))^2)/(R^2 + 2*R*y4 + x4^2 + y4^2)^2,
(-r5^2*(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2 + (x2 - x4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2 + (y2 - y4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2)/(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2,
(-R^2*r5^2 + R^2*x3^2 - 2*R^2*x3*x5 + R^2*x5^2 - 2*R*r5^2*y3 + 2*R*x3^2*y5 - 2*R*x3*x5*y3 - 2*R*x3*x5*y5 + 2*R*x5^2*y3 - r5^2*x3^2 - r5^2*y3^2 + x3^2*y5^2 - 2*x3*x5*y3*y5 + x5^2*y3^2)/(R^2 + 2*R*y3 + x3^2 + y3^2),
]
end;

iniv = [big"-23.0", -22, -15, 28, 31, -9, 3, -22, 6, 32]
res = nls(H, ini=iniv)
println(res)
([-23.076923076923077, -22.884615384615383, -15.507692307692308, 28.56153846153846, 31.2, -9.1, 3.6, -22.3, 6.0, 32.5], true)

\(x_2 = -23.07692,\ y_2 = -22.88462,\ x_3 = -15.50769, y_3 = 28.56154\)
\(x_4 = 31.20000,\ y_4 = -9.10000,\ x_5 = 3.60000,\ y_5 = -22.30000\)
\(r_5 = 6.00000,\ R = 32.50000\)
となり，求める地斜の長さは三平方の定理より 52 と計算される。

地斜 \(= \sqrt{(x_3 - x_2)^2 + (y_3 - y_2)^2} = 52.00000\)

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

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x2, y2, x3, y3, x4, y4, x5, y5, r5, R) = res[1]
x1 = 0
y1 = -R
@printf("x1 = %.5f, y1 = %.5f, x2 = %.5f, y2 = %.5f, x3 = %.5f, y3 = %.5f\nx4 = %.5f, y4 = %.5f, x5 = %.5f, y5 = %.5f, r5 = %.5f, R = %.5f\n", x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, r5, R)
@printf("地斜 = sqrt( (x3 - x2)^2 + (y3 - y2)^2) = %.5f\n", sqrt( (x3 - x2)^2 + (y3 - y2)^2))
plot([x1, x2, x3, x4, x1], [y1, y2, y3, y4, y1], color=:blue, lw=0.5)
circle(0, 0, R)
circle(x5, y5, r5)
segment(x1, y1, x3, y3, :green)
segment(x2, y2, x4, y4, :green)
if more
point(x1, y1, "(x1,y1)")
point(x2, y2, "(x2,y2) ", :green, :right)
point(x3, y3, " (x3,y3)")
point(x4, y4, " (x4,y4)")
point(x5, y5, "(x5,y5;r5)", :green, :center)
point(x5, y5, "小円", :red, :center, :bottom)
point(R, 0, " R: 大円", :red)
point( (x3 + x4)/2, (y3 + y4)/2, " 左斜", :blue, mark=false)
point( (x2 + x3)/2, (y2 + y3)/2, " 地斜", :blue, mark=false)
point( (x1 + x2)/2, (y1 + y2)/2, "右斜", :blue, :right, mark=false)
point( (x1 + x4)/2, (y1 + y4)/2, "天斜", :blue, mark=false)
point( (x1 + x3)/2, (y1 + y3)/2, " 乾斜", :green, mark=false)
point( (x2 + x4)/2, (y2 + y4)/2, " 坤斜", :green, mark=false)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;