算額（その296）
中村信弥「算額への招待 追補2」
http://www.wasan.jp/syotai/syotai.html
長野県千曲市八幡 武水別神社（八幡神社） 天保8年(1837)
天円，地円，人円が互いに外接し，三角形に内接している。天円と地円の直径がそれぞれ 8 寸，4.5 寸のとき，三角形の三辺の長さを求めよ。
図に示すように記号を定め，以下の連立方程式を nlsolve() により解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, r2::positive, r3::positive,
x1::positive, x2::positive, x3::positive, x32::positive,
a::positive, b::positive, c::positive;
(r1, r2) = (8//2, 9//4)
eq1 = r3/x3 - r2/x2
eq2 = r3/x3 - r1/x1
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x32 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq6 = (a + sqrt((a - b)^2 + c^2) + sqrt(b^2 + c^2)) * r1 - a * c
eq7 = r3 / (a - x32) - r1 / (a - x1)
eq8 = distance(b, c, a, 0, x1, r1) - r1^2
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])
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)
(r3, x1, x2, x3, x32, a, b, c) = u
return [
r3/x3 - 9/(4*x2), # eq1
r3/x3 - 4/x1, # eq2
(9/4 - r3)^2 - (r3 + 9/4)^2 + (x2 - x3)^2, # eq3
(x1 - x2)^2 - 36, # eq4
(4 - r3)^2 - (r3 + 4)^2 + (-x1 + x32)^2, # eq5
-a*c + 4*a + 4*sqrt(b^2 + c^2) + 4*sqrt(c^2 + (a - b)^2), # eq6
r3/(a - x32) - 4/(a - x1), # eq7
(x1 - (a^2*x1 - 2*a*b*x1 + a*c^2 - 4*a*c + b^2*x1 + 4*b*c)/(a^2 - 2*a*b + b^2 + c^2))^2 + (-c*(a^2 - a*b - a*x1 + b*x1 + 4*c)/(a^2 - 2*a*b + b^2 + c^2) + 4)^2 - 16, # eq8
]
end;
iniv = [big"1.3", 14, 8, 4, 18, 20, 15, 10]
res = nls(H, ini=iniv);
using Printf
names = ("r3", "x1", "x2", "x3", "x32", "a", "b", "c")
for (i, name) in enumerate(names)
@printf("%s = %.6f\n", name, res[1][i])
end
println(res[2])
r3 = 1.265625
x1 = 13.714286
x2 = 7.714286
x3 = 4.339286
x32 = 18.214286
a = 20.297143
b = 15.250421
c = 9.723228
true
大斜 = a = 20.29714285714286
中斜 = sqrt(b^2 + c^2) = 18.08636238036625
小斜 = sqrt((a - b)^2 + c^2)) = 10.954933808937678
算額の答えは，大斜 = 41.427，中斜 = 29.533，小斜 = 15.534 であるが，これはありえない。図を描くとわかる。
中村信弥「算額への招待 追補2」では，大斜の式として，天地人の直径をそれぞれ a, x, b として
l = a/(a - b)*(sqrt(a*sqrt(a*b)) + sqrt(b*sqrt(a*b)) + sqrt(a*b))
を導いている。この式は正しいが，
「a = 8, b = 4.5 とすると l = 16*(6 + 7√3)/7 = 41.427097...」としているが，4.5 は地円の直径であり，人円の直径はこの時点では未知である。
実は，途中で x = sqrt(a*b) とあるので，既知の a, x から b を求めればよい。すなわち，b = x^2/a である。b = 4.5 を代入したために間違った解を求めてしまったのである。
上の式で，a = 8, b = 4.5^2/8 = 2.53125 を代入すれば，大斜 = 20.297142857142855 になる。元の算額の術も同じ間違いをおかしている。
a = 8
b = 4.5^2/8 # = 2.53125
l = a/(a - b)*(sqrt(a*sqrt(a*b)) + sqrt(b*sqrt(a*b)) + sqrt(a*b))
20.297142857142855
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r3, x1, x2, x3, x32, a, b, c) = res[1]
println("r1 = $(Float64(r1)); r2 = $(Float64(r2))")
println("大斜 = $(Float64(a)); 中斜 = $(Float64(sqrt(b^2 + c^2))); 小斜 = $(Float64(sqrt((a - b)^2 + c^2)))")
plot([0, a, b, 0], [0, 0, c, 0], color=:black, lw=0.5)
circle(x1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :green)
circle(x32, r3, r3, :green)
if more
point(x1, r1, "天:r1,(x1,r1)", :red, :center)
point(x2, r2, "地:r2,(x2,r2)", :blue, :center)
point(x3, r3, "人:r3,(x3,r3)", :green, :center)
point(x32, r3, "人:r3,(x32,r3)", :green, :center)
point(b, c, " (b,c)", :green, :left, :bottom)
point(a, 0, " a", :green, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
