算額（その465）
愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)
平田浩一: 伊佐爾波神社の算額にみる江戸末期の和算，愛媛大学教育学部紀要，第60巻，195-206, 2013.
https://www.ed.ehime-u.ac.jp/~kiyou/2013/pdf/20.pdf
正三角形の中に2本の斜線で区切られた領域それぞれに直径の等しい円が内接している。下側の斜線と底辺を左に伸ばし，その2直線と正三角形の左側の斜辺に内接する外接円の直径を求めよ。
正三角形の左隅を原点におき，一辺の長さを a とする。
3個の等円の半径と中心座標を r, (x1, r), (x2, y2), (x3, y3)
外円の半径と中心座標を r0, (x0, r0)
斜線と正三角形の斜辺の交点座標を (x4, √3x4), (x5, (a - x5)√3)
とおき以下の連立方程式の数値解を求める。
include("julia-source.txt")
using SymPy
@syms a::positive, r::positive, x1::positive, x2::negative, y2::positive,
x3::positive, y3::negative, x4::positive, x5::positive,
r0::positive, x0::negative;
eq1 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x1, r) - r^2
eq2 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x3, y3) - r^2
eq3 = distance(a, 0, a/2, sqrt(Sym(3))a/2, x2, y2) - r^2
eq4 = distance(a, 0, a/2, sqrt(Sym(3))a/2, x3, y3) - r^2
eq5 = distance(x4, sqrt(Sym(3))x4, x5, (a - x5)sqrt(Sym(3)), x2, y2) - r^2
eq6 = distance(x4, sqrt(Sym(3))x4, x5, (a - x5)sqrt(Sym(3)), x3, y3) - r^2
eq7 = distance(x4, sqrt(Sym(3))x4, a, 0, x1, r) - r^2
eq8 = distance(x4, sqrt(Sym(3))x4, a, 0, x2, y2) - r^2
eq9 = distance(x4, sqrt(Sym(3))x4, a, 0, x0, r0) - r0^2 # 外円
eq10 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x0, r0) - r0^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=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 v, r.f_converged
end;
function H(u)
(r0, x0, r, x1, x2, y2, x3, y3, x4, x5) = u
return [
-r^2 + (r/4 - sqrt(3)*x1/4)^2 + (-sqrt(3)*r/4 + 3*x1/4)^2, # eq1
-r^2 + (3*x3/4 - sqrt(3)*y3/4)^2 + (-sqrt(3)*x3/4 + y3/4)^2, # eq2
-r^2 + (-3*a/4 + 3*x2/4 + sqrt(3)*y2/4)^2 + (-sqrt(3)*a/4 + sqrt(3)*x2/4 + y2/4)^2, # eq3
-r^2 + (-3*a/4 + 3*x3/4 + sqrt(3)*y3/4)^2 + (-sqrt(3)*a/4 + sqrt(3)*x3/4 + y3/4)^2, # eq4
-r^2 + (x2 - (3*a^2*x4 - 3*a*x4^2 - 9*a*x4*x5 - sqrt(3)*a*x4*y2 + sqrt(3)*a*x5*y2 + x2*x4^2 - 2*x2*x4*x5 + x2*x5^2 + 6*x4^2*x5 + sqrt(3)*x4^2*y2 + 6*x4*x5^2 - sqrt(3)*x5^2*y2)/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2 + (y2 - (y2*(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2) + (x4 - x5)*(-sqrt(3)*a*x2 + sqrt(3)*a*x4 + sqrt(3)*x2*x4 + sqrt(3)*x2*x5 - 2*sqrt(3)*x4*x5 - x4*y2 + x5*y2))/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2, # eq5
-r^2 + (x3 - (3*a^2*x4 - 3*a*x4^2 - 9*a*x4*x5 - sqrt(3)*a*x4*y3 + sqrt(3)*a*x5*y3 + x3*x4^2 - 2*x3*x4*x5 + x3*x5^2 + 6*x4^2*x5 + sqrt(3)*x4^2*y3 + 6*x4*x5^2 - sqrt(3)*x5^2*y3)/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2 + (y3 - (y3*(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2) + (x4 - x5)*(-sqrt(3)*a*x3 + sqrt(3)*a*x4 + sqrt(3)*x3*x4 + sqrt(3)*x3*x5 - 2*sqrt(3)*x4*x5 - x4*y3 + x5*y3))/(3*a^2 - 6*a*x4 - 6*a*x5 + 4*x4^2 + 4*x4*x5 + 4*x5^2))^2, # eq6
-r^2 + (r - sqrt(3)*x4*(a^2 - a*x1 - a*x4 + sqrt(3)*r*x4 + x1*x4)/(a^2 - 2*a*x4 + 4*x4^2))^2 + (x1 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x1 + a*x4 + sqrt(3)*r*x4 + x1*x4 - 4*x4^2))/(a^2 - 2*a*x4 + 4*x4^2))^2, # eq7
-r^2 + (x2 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x2 + a*x4 + x2*x4 - 4*x4^2 + sqrt(3)*x4*y2))/(a^2 - 2*a*x4 + 4*x4^2))^2 + (-sqrt(3)*x4*(a^2 - a*x2 - a*x4 + x2*x4 + sqrt(3)*x4*y2)/(a^2 - 2*a*x4 + 4*x4^2) + y2)^2, # eq8
-r0^2 + (r0 - sqrt(3)*x4*(a^2 - a*x0 - a*x4 + sqrt(3)*r0*x4 + x0*x4)/(a^2 - 2*a*x4 + 4*x4^2))^2 + (x0 - (x4*(a^2 - 2*a*x4 + 4*x4^2) - (a - x4)*(-a*x0 + a*x4 + sqrt(3)*r0*x4 + x0*x4 - 4*x4^2))/(a^2 - 2*a*x4 + 4*x4^2))^2, # eq9
-r0^2 + (r0/4 - sqrt(3)*x0/4)^2 + (-sqrt(3)*r0/4 + 3*x0/4)^2, # eq10
]
end;
a = 20
iniv = [big"16.0", -10, 10, 19, 47, 21, 35, 42, 14, 52] .* (a/71)
iniv = a .* [big"0.225", -0.141, 0.141, 0.268, 0.662, 0.296, 0.493, 0.592, 0.197, 0.732]
res = nls(H, ini=iniv);
names = ("r0", "x0", "r", "x1", "x2", "y2", "x3", "y3", "x4", "x5")
if res[2]
for (name, x) in zip(names, res[1])
@printf("%s = %g; ", name, round(Float64(x), digits=6))
end
println()
else
println("収束していない")
end;
r0 = 4.50197; x0 = -2.59921; r = 2.96213; x1 = 5.13055; x2 = 13.0676; y2 = 6.08309; x3 = 10; y3 = 11.3963; x4 = 3.86488; x5 = 14.4977;
三角形の一辺の長さが 20 寸のとき，外円の直径は 9.00393 寸である。ちなみに，等円の直径は 5.92425 寸である。
using Plots
function abline(x0, y0, slope, minx, maxx, color=:black; lw=0.5)
f(x) = slope * (x - x0) + y0
segment(minx, f(minx), maxx, f(maxx), color; lw=lw)
end;
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 20
(r0, x0, r, x1, x2, y2, x3, y3, x4, x5) = res[1]
@printf("r0 = %g; x0 = %g; r = %g; x1 = %g; x2 = %g; y2 = %g; x3 = %g; y3 = %g; x4 = %g; x5 = %g\n", r0, x0, r, x1, x2, y2, x3, y3, x4, x5)
@printf("外円の直径 = %g; 等円の直径 = %g\n", 2r0, 2r)
plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:green, lw=0.5)
circle(x1, r, r, :blue)
circle(x2, y2, r, :blue)
circle(x3, y3, r, :blue)
circle(x0, r0, r0)
segment(x4, √3x4, x5, (a - x5)√3)
abline(x4, √3x4, √3x4/(x4 - a), -0.4a, a)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(x1, r, "(x1,r)", :blue, :center, :top, delta=-delta)
point(x2, y2, "(x2,y2)", :blue, :center, :top, delta=-delta)
point(x3, y3, "(x3,y3)", :blue, :center, :top, delta=-delta)
point(x0, r0, "(x0,r0)", :red, :center, :top, delta=-delta)
point(x4, √3x4, " (x4,√3x4)", :black, :left, :vcenter)
point(x5, √3(a-x5), " (x5,√3(a-x5))", :black, :left, :vcenter)
point(a/2, √3a/2, " (a/2,√3a/2)", :black, :left, :vcenter)
hline!([0], color=:gray, lw=0.5)
vline!([0], color=:gray, lw=0.5)
else
plot!(showaxis=false)
end
end;
