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

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

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

問1. 長方形内に 2 本の斜線を容れ，それぞれに接する甲円と乙円を容れる。長方形の短辺の長さが 10 寸，甲円の直径が 8 寸のとき，長方形の長辺の長さはいかほどか。

長方形の長辺の長さを \(x\)，乙円の半径を \(r\)，斜線が \(x\) 軸，\(y\) 軸と交わる座標を \(x_1,\ y_1\) とする。各円の中心から斜線までの距離が半径に等しいことを条件として，方程式を解く。

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

using SymPy

@syms x::positive, x1::positive, y1::positive, r::positive;

eq1 = distance(0, y1, x, 10, r, r) - r^2 |> simplify
eq2 = distance(0, 10, x1, 0, r, r) - r^2 |> simplify
eq3 = distance(0, y1, x, 10, x - 4, 4) - 4^2 |> simplify
eq4 = distance(0, 10, x1, 0, x - 4, 4) - 4^2 |> simplify

res = solve([eq1, eq2, eq3, eq4], (x, x1, y1, r))[1]

(8, 24, 20/3, 4)

今までに経験のないことであるが，solve で得られる解が不適切解である。得られる乙円の半径が 4 という解は甲円と乙円が重なるという解である。

そこで，nlsolve() で解いてみる。

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")

x*(2*r^2*y1 - 20*r^2 - 2*r*x*y1 - 2*r*y1^2 + 20*r*y1 + x*y1^2)/(x^2 + y1^2 - 20*y1 + 100),
20*x1*(r^2 - r*x1 - 10*r + 5*x1)/(x1^2 + 100),
4*x*(5*x + 12*y1 - 120)/(x^2 + y1^2 - 20*y1 + 100),
20*(5*x^2 - 6*x*x1 - 40*x + x1^2 + 24*x1)/(x1^2 + 100),

function H(u)
(x, x1, y1, r) = u
return [
x*(2*r^2*y1 - 20*r^2 - 2*r*x*y1 - 2*r*y1^2 + 20*r*y1 + x*y1^2)/(x^2 + y1^2 - 20*y1 + 100),
20*x1*(r^2 - r*x1 - 10*r + 5*x1)/(x1^2 + 100),
4*x*(5*x + 12*y1 - 120)/(x^2 + y1^2 - 20*y1 + 100),
20*(5*x^2 - 6*x*x1 - 40*x + x1^2 + 24*x1)/(x1^2 + 100),
]
end;

iniv = [big"13.8", 7.8, 4.3, 2.6]
res = nls(H, ini=iniv)
println(res)
([13.759000729485331, 7.807200583588266, 4.267083029381112, 2.560249817628667], true)

\(x = 13.75900,\ x_1 = 7.80720,\ y_1 = 4.26708,\ r = 2.56025\)
であり，\(x = 13.75900\) は元の単位でいえば 13寸7分5厘9毛である。

算額の答えでは 13 寸 7 分 5 厘 6 毛あまりありとなっている。

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

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x, x1, y1, r) = res[1]
@printf("x = %.5f, x1 = %.5f, y1 = %.5f, r = %.5f\n", x, x1, y1, r)
plot([0, x, x, 0, 0], [0, 0, 10, 10, 0], color=:blue, lw=0.5)
circle(r, r, r)
circle(x - 4, 4, 4, :blue)
segment(0, 10, x1, 0, :green)
segment(0, y1, x, 10, :green)
if more
point(x - 4, 4, " 甲: (x-4,4)")
point(r, r, " 乙: (r,r)")
point(x1, 0, "x1")
point(0, y1, "y1 ", :green, :right)
point(x, 0, " x")
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5, xlims=(-1, 14.5),ylims=(-0.5, 10.5))
end
end;