算額(その1231)
京都府京都市上京区北野馬喰町 北野天満宮 明治12年(1879)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円4個,二等辺三角形
二等辺三角形の中に,甲円 2 個,乙円と丙円を 1 個ずつ容れる。上斜(斜辺)が 1014 寸,下斜(底辺)が 1428 寸のとき,甲円の直径はいかほどか。
斜辺,底辺を c, 2a とおく。高さは b = sqrt(c^2 - a^2) となる。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (0, 2r3 + r2)
丙円の半径と中心座標を r3, (0, r3)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, c::positive,
r1::positive, x1::positive,
r2::positive, r3::positive
@syms a, b, c, r1, x1, r2, r3
# b = sqrt(c^2 - a^2) 式には代入しない
eq1 = dist2(a, 0, 0, b, 0, 2r3 + r2, r2)
eq2 = dist2(a, 0, 0, b, x1, r1, r1)
eq3 = x1^2 + (2r3 + r2 - r1)^2 - (r1 + r2)^2
eq4 = x1^2 + (r1 - r3)^2 - (r1 + r3)^2;
# res = solve([eq1, eq2, eq3], (r1, r2, r3))
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(r1, x1, r2, r3) = u
return [
a^2*b^2 - 2*a^2*b*r2 - 4*a^2*b*r3 + 4*a^2*r2*r3 + 4*a^2*r3^2 - b^2*r2^2, # eq1
b*(a^2*b - 2*a^2*r1 - 2*a*b*x1 + 2*a*r1*x1 - b*r1^2 + b*x1^2), # eq2
x1^2 - (r1 + r2)^2 + (-r1 + r2 + 2*r3)^2, # eq3
x1^2 + (r1 - r3)^2 - (r1 + r3)^2, # eq4
]
end;
a = 1428//2
c = 1014
b = sqrt(c^2 - a^2)
iniv = BigFloat[180, 290, 203, 115]
res = nls(H, ini=iniv)
([178.5, 285.6, 203.09333333333333, 114.24], true)
甲円,乙円,丙円の直径はそれぞれ 357 寸,406.18666666666667寸,228.48 寸である。
function draw(a, c, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
#b = sqrt(c^2 - a^2)
(r1, x1, r2, r3) = res[1]
@printf("上斜が %g,下斜が %g のとき,甲円の直径は %g である。\n", 2a, c, 2r1)
plot([a, 0, -a, a], [0, b, 0, 0], color=:green, lw=0.5)
circle2(x1, r1, r1)
circle(0, 2r3 + r2, r2, :blue)
circle(0, r3, r3, :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(0, b, " b", :green, :left, :bottom, delta=delta)
point(a, 0, " a", :green, :left, :bottom, delta=delta)
point(x1, r1, "甲円:r1\n(x1,r1)", :red, :center, delta=-delta)
point(0, 2r3 + r2, "乙円:r2\n(0,2r3+r2)", :blue, :center, delta=-delta)
point(0, r3, "丙円:r3\n(0,r3)", :magenta, :center, delta=-delta)
end
end;
draw(1428/2, 1014, true)