静岡県磐田市鎌田 医王寺（鎌田山金剛院醫王寺） 安永8年(1779)
https://www.iouji.net/
https://www.iouji.net/visit/treasure/t02/
キーワード：円3個，直角三角形，正方形
#Julia #SymPy #算額 #和算 #数学
直角三角形の中に正方形と大円，中円，小円を容れる。大中小 3 円の直径の和が 94 寸，斜辺と正方形の一辺の差が 125 寸のとき，正方形の一辺の長さはいかほどか。
大円，中円，小円の半径を \(r_1,\ r_2, r_3\)
正方形の一辺の長さを \(s\)
直角三角形の 3 辺の長さを \(c,\ d,\ g\)
大中小 3 円の直径の和を \(K_1\)
斜辺と正方形の一辺の差を \(K_2\)
そのほか図のように記号を定め，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, c::positive, d::positive,
e::positive, f::positive, g::positive, s::positive,
r1::positive, r2::positive, r3::positive
K1 = 94
K2 = 125
eq1 = 2(r1 + r2 + r3) - K1
eq2 = g - s - K2
eq3 = s + f - (c - a) - 2r1
eq4 = s + e - (d - b) - 2r2
eq5 = a + b - s - 2r3
eq6 = (c - a)/s - r1/r3
eq7 = (d - b)/(c - a) - r2/r1
eq8 = (s + e + f) - g
eq9 = c^2 + d^2 - (s + e + f)^2
eq10 = a^2 + b^2 - s^2
eq11 = (s*f + s*e + a*b)/2 + s^2 - c*d/2;
\(K_1 = 94,\ K_2 = 125\) と実値を与えれば，連立方程式を解くことができ，すべてのパラメータは整数値になる。
このときの図が上に示したものである。
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11],
(a, b, c, d, e, f, g, s, r1, r2, r3))[4] # 4 of 4
(48, 36, 148, 111, 45, 80, 185, 60, 20, 15, 12)
\(K_1,\ K_2\) に任意の値を与えた場合，方程式は解けないこともある（有限の時間内に解が得られない）。
初期値が適切であれば数値解を求めることができる（多くの場合，関数中で定義した初期値で十分のようである）。
\(K_1 = 94,\ K_2 = 125\) の場合の数値解は，解析解と同じものが得られる。
\(正方形の一辺の長さ = 60\)
\(直径の和 = 94,\ 正方形の一辺と弦の差 = 125\)
\(a = 48,\ b = 36,\ c = 148,\ d = 111,\ e = 45\)
\(f = 80,\ g = 185,\ s = 60\)
\(r_1 = 20,\ r_2 = 15,\ r_3 = 12\)
function driver(K1, K2)
function H(u)
(a, b, c, d, e, f, g, s, r1, r2, r3) = u
return [
2(r1 + r2 + r3) - K1, # eq1
g - s - K2, # eq2
s + f - (c - a) - 2r1, # eq3
s + e - (d - b) - 2r2, # eq4
a + b - s - 2r3, # eq5
(c - a)/s - r1/r3, # eq6
(d - b)/(c - a) - r2/r1, # eq7
(s + e + f) - g, # eq8
c^2 + d^2 - (s + e + f)^2, # eq9
a^2 + b^2 - s^2, # eq10
(s*f + s*e + a*b)/2 + s^2 - c*d/2, # eq11
]
end;
iniv = 37 *BigFloat[1.29836, 0.97287, 4, 3, 1.2178, 2.16163, 5,
1.62057, 0.53916, 0.4032, 0.32628]
res = nls(H, ini=iniv)[1]
end;
術は，以下の \(x\) についての 5 次式を解けばよいとしている。
\( ( ( (2x^2 + 2K_1*K_2 - 8K_2^2)*x +6K_1*K_2^2 - 8K_2^3)*x + 6K_1*K_2^3 - 2K_2^4)*x + 2K_1*K_2^4 - (x+ K_2)^2*K_1^2*K_2 = 0\)
@syms K1, K2, x
f = ( ( (2x^2 + 2K1*K2 - 8K2^2)*x +6K1*K2^2 - 8K2^3)*x + 6K1*K2^3 - 2K2^4)*x + 2K1*K2^4 - (x+ K2)^2*K1^2*K2
\(- K_{1}^{2} K_{2} \left(K_{2} + x\right)^{2} + 2 K_{1} K_{2}^{4} + x \left(6 K_{1} K_{2}^{3} - 2 K_{2}^{4} + x \left(6 K_{1} K_{2}^{2} - 8 K_{2}^{3} + x \left(2 K_{1} K_{2} - 8 K_{2}^{2} + 2 x^{2}\right)\right)\right)\)
\(f\) は \( (K_2 + x)\) という因子を持つので，実際には \(x\) の 3 次式を解けばよい。
f |> factor
\(\left(K_{2} + x\right)^{2} \left(- K_{1}^{2} K_{2} + 2 K_{1} K_{2}^{2} + 2 K_{1} K_{2} x - 2 K_{2}^{2} x - 4 K_{2} x^{2} + 2 x^{3}\right)\)
g = -K1^2*K2 + 2*K1*K2^2 + 2*K1*K2*x - 2*K2^2*x - 4*K2*x^2 + 2*x^3
\(- K_{1}^{2} K_{2} + 2 K_{1} K_{2}^{2} + 2 K_{1} K_{2} x - 2 K_{2}^{2} x - 4 K_{2} x^{2} + 2 x^{3}\)
\(K_1 = 94,\ K_2 = 125\) の場合には確かに，正方形の一辺の長さは 60 という正解を与える。
solve(g(K1 => 94, K2 => 125))[1]
\(60\)
また，たとえば \(K_1 = 80,\ K_2 = 160\) の場合にも，正方形の一辺の長さは 55.0399356501993 という正解を与える。
res = solve(g(K1 => 80, K2 => 150))[1].evalf()
\(55.0399356501993 - 7.0 \cdot 10^{-21} i\)
上の数値解を求めるプログラムでの結果と一致する。
正方形の一辺の長さ = 55.0399356501993
直径の和 = 80, 正方形の一辺と弦の差 = 150
a = 50.4343, b = 22.0402, c = 187.883, d = 82.1063, e = 24.0529
f = 125.947, g = 205.04, s = 55.0399
r1 = 21.7693, r2 = 9.51337, r3 = 8.71732
描画関数プログラムのソースを見る
function draw(K1, K2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, c, d, e, f, g, s, r1, r2, r3) = driver(K1, K2)
@printf("正方形の一辺の長さ = %.15g\n", s)
@printf("直径の和 = %.15g, 正方形の一辺と弦の差 = %g\n", 2(r1 + r2 + r3), g - s)
@printf("a = %g, b = %g, c = %g, d = %g, e = %g\nf = %g, g = %g, s = %g\nr1 = %g, r2 = %g, r3 = %g\n",
a, b, c, d, e, f, g, s, r1, r2, r3)
plot([0, c, 0, 0], [0, 0, d, 0], color=:green, lw=0.5)
plot!([b, 0, a, a + b], [b + a, b, 0, a], color=:green, lw=0.5)
circle(r3, r3, r3)
x1 = c - r1/tan(atan(d, c)/2)
circle(x1, r1, r1, :blue)
y2 = d - r2/tan(atan(c, d)/2)
circle(r2, y2, r2, :magenta)
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(x1, r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta)
point(r2, y2, "中円:r2\n(r2,y2)", :magenta, :center, delta=-delta)
point(r3, r3, "小円:r3\n(r3,r3)", :red, :center, delta=-delta)
point(a, 0, "a", :green, :center, delta=-delta)
point(c, 0, "c", :green, :center, delta=-delta)
point(0, b, "b ", :green, :right, :vcenter)
point(0, d, "d ", :green, :right, :vcenter)
dimension_line(0, d, b, a + b, "e", :tomato, delta=2delta, deltax=2delta, length=s/15)
dimension_line(b, a + b, a + b, a, "s", :darkorange, delta=2delta, deltax=2delta, length=s/15)
dimension_line(a + b, a, c, 0, "f", :purple, delta=2delta, deltax=2delta, length=s/15)
plot!(xlims=(-5delta, c + delta), ylims=(-5delta, d + delta))
end
end;
draw(94, 125, true)
正方形の一辺の長さ = 60
直径の和 = 94, 正方形の一辺と弦の差 = 125
a = 48, b = 36, c = 148, d = 111, e = 45
f = 80, g = 185, s = 60
r1 = 20, r2 = 15, r3 = 12
以下のアイコンをクリックして応援してください
