算額(その139)
大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html
鉤,股,弦がそれぞれ 3寸1分,3寸6分,4寸8分である三角形が2本の斜線で区切られ,径の同じ円が 2 個入っている。円の径を求めよ。
鉤股弦とはいっているが,直角三角形ではない。鉤と股をそのままの数値として受け入れれば,弦は 4.750789408087881 寸でなければならない。
鉤 = 31, 股 = 36 円の半径を r,それぞれの中心座標を (r, y), (x, r) とおき,直線までの距離に関する方程式を立て,解く。
solve() では解けないので,nlsolve() で数値解を求める。
using SymPy
function distance(x1, y1, x2, y2, x0, y0)
l = sympy.Line(sympy.Point(x1, y1), sympy.Point(x2, y2))
l.distance(sympy.Point(x0, y0))^2 # 二乗距離を返す!!
end;
@syms r::positive, x::positive, y::positive, x0::positive, y0::positive;
eq1 = distance(0, 31, x0, y0, r, y) - r^2 # 円1の中心から BD への距離
eq2 = distance(0, 31, x0, y0, x, r) - r^2 # 円2の中心から BD への距離
eq3 = distance(0, 31, 36, 0, x, r) - r^2 # 円2の中心から AB への距離
eq4 = distance(0, 0, x0, y0, r, y) - r^2 # 円1の中心から OC への距離
eq5 = distance(0, 0, x0, y0, x, r) - r^2; # 円2の中心から OC への距離
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x, y, x0, y0))
println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
-r^2 + (r - x0*(r*x0 + y*y0 - 31*y - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (y - (r*x0*y0 - 31*r*x0 + 31*x0^2 + y*y0^2 - 62*y*y0 + 961*y)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (r - (r*y0^2 - 62*r*y0 + 961*r + x*x0*y0 - 31*x*x0 + 31*x0^2)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (x - x0*(r*y0 - 31*r + x*x0 - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (1116*r/2257 + 961*x/2257 - 34596/2257)^2 + (1296*r/2257 + 1116*x/2257 - 40176/2257)^2,
-r^2 + (r - x0*(r*x0 + y*y0)/(x0^2 + y0^2))^2 + (y - y0*(r*x0 + y*y0)/(x0^2 + y0^2))^2,
-r^2 + (r - y0*(r*y0 + x*x0)/(x0^2 + y0^2))^2 + (x - x0*(r*y0 + x*x0)/(x0^2 + y0^2))^2,
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=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini)#, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(r, x, y, x0, y0) = u
return [
-r^2 + (r - x0*(r*x0 + y*y0 - 31*y - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (y - (r*x0*y0 - 31*r*x0 + 31*x0^2 + y*y0^2 - 62*y*y0 + 961*y)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (r - (r*y0^2 - 62*r*y0 + 961*r + x*x0*y0 - 31*x*x0 + 31*x0^2)/(x0^2 + y0^2 - 62*y0 + 961))^2 + (x - x0*(r*y0 - 31*r + x*x0 - 31*y0 + 961)/(x0^2 + y0^2 - 62*y0 + 961))^2,
-r^2 + (1116*r/2257 + 961*x/2257 - 34596/2257)^2 + (1296*r/2257 + 1116*x/2257 - 40176/2257)^2,
-r^2 + (r - x0*(r*x0 + y*y0)/(x0^2 + y0^2))^2 + (y - y0*(r*x0 + y*y0)/(x0^2 + y0^2))^2,
-r^2 + (r - y0*(r*y0 + x*x0)/(x0^2 + y0^2))^2 + (x - x0*(r*y0 + x*x0)/(x0^2 + y0^2))^2,
]
end;
iniv = [5.0, 20, 10, 12, 10]
res = nls(H, ini=iniv)
println(res)
([5.6146321148429355, 20.875286969374045, 9.74605295956063, 13.244959542108463, 7.680342537201841], true)
直径は 2 × 5.6146321148429355 ≒ 11寸2分2厘9毛である。算額では 9分4厘5毛になっている。
~井於神社奉納「算額」について~ http://io-jinja.org/sangakusetumei.html に引用されている,「算術稽古大全」宅間流五世松岡能一,文化5年(1808) に以下が示されているとあるが,これはまた以上の 2 通りの解とも違う。
(鉤, 股, 弦) = (3.1, 3.6, 4.8)
天 = 鉤 + 股 + 弦
等円径 = (sqrt(2天 * 鉤) + 天) / 2(鉤*股)
0.8935453733438091
弦に正しい(?)数値を代入しても異なる解になる。
(鉤, 股, 弦) = (3.1, 3.6, 4.750789408087881)
天 = 鉤 + 股 + 弦
等円径 = (sqrt(2天 * 鉤) + 天) / 2(鉤*股)
0.8905302961492191
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, lw=0.5, linestyle=:solid)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=lw, linestyle=linestyle)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function intersection(x1, y1, x2, y2, x3, y3, x4, y4)
denominator = (x1*y3 - x1*y4 - x2*y3 + x2*y4 - x3*y1 + x3*y2 + x4*y1 - x4*y2)
X = (x1*x3*y2 - x1*x3*y4 - x1*x4*y2 + x1*x4*y3 - x2*x3*y1 + x2*x3*y4 + x2*x4*y1 - x2*x4*y3) / denominator
Y = (x1*y2*y3 - x1*y2*y4 - x2*y1*y3 + x2*y1*y4 - x3*y1*y4 + x3*y2*y4 + x4*y1*y3 - x4*y2*y3) / denominator
(X, Y)
end
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, x, y, x0, y0) = res[1]
@printf("円の直径 = %.3f\n", 2r)
plot([0, 36, 0, 0], [0, 0, 31, 0])
circle(r, y, r)
circle(x, r, r)
(X1, Y1) = intersection(0, 0, x0, y0, 0, 31, 36, 0)
segment(0, 0, X1, Y1)
(X2, Y2) = intersection(0, 31, x0, y0, 0, 0, 36, 0)
segment(0, 31, X2, Y2)
if more == true
point(r, y, "円1(r,y;r)", :red, :center, :top)
point(x, r, "円2(x,r;r)", :red, :center, :top)
point(x0, y0, "(x0,y0)", :red)
point(0, 0, "O ", :blue, :right, :bottom)
point(36, 0, " A", :blue, :left, :bottom)
point(0, 31, " B", :blue, :left, :bottom)
point(X1, Y1, " C", :blue, :left, :bottom)
point(X2, Y2, "D ", :blue, :right, :bottom)
hline!([0], color=:black, lw=0.5, xlims=(-2, 38))
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
draw(false)
円の直径 = 11.229