算額（その280）
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(240)
長野県松本市筑摩 筑摩神社 明治6年(1873)
鉤股弦（直角三角形）に 2 本の斜線を入れ，区切られた領域に 2 個の等円を入れる。鉤が 5 寸のとき，等円の直径を求めよ。
図のように記号を定め，連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, r::positive, x1::positive, y1::positive,
x::positive, y::positive;
eq1 = distance(0, 5, a, 0, r, y1) - r^2
eq2 = distance(0, 5, a, 0, x1, r) - r^2
eq3 = distance(0, 0, x, y, r, y1) - r^2
eq4 = distance(0, 0, x, y, x1, r) - r^2
eq5 = distance(0, 5, a + b, 0, x1, r) - r^2
eq6 = (5 - y)*(a + b) - 5x
eq7 = x*y + x*(5-y)/2 + (a + b - x)*y/2 - 5(a + b)/2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])
nlsove() で解く。
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
-r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2, # eq1
-r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2, # eq2
-r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2, # eq3
-r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2, # eq4
-r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2, # eq5
-5*x + (5 - y)*(a + b), # eq6
-5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2, # eq7
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)
(a, b, r, x1, y1, x, y) = u
return [
-r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2, # eq1
-r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2, # eq2
-r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2, # eq3
-r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2, # eq4
-r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2, # eq5
-5*x + (5 - y)*(a + b), # eq6
-5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2, # eq7
]
end;
iniv = [big"2.8", 2.9, 0.9, 3.3, 1.6, 3.4, 2.1]
res = nls(H, ini=iniv);
println(res);
(BigFloat[2.800221924622118970432095313760788507835397409564486177867680366139659180899345, 2.91487027133284102199601750983229645815568299593371469024204485373220221457499, 0.8974814920915777221723012652934296052881492161222889500985443442946459740211824, 3.326236908516905783403296227571694673430896734617166839237861420836889759024065, 1.560761971632353997076535469661088458153271795581832971123527071275699538470573, 3.431977890341429946507935106849530222260823283947106463592744755845496558390283, 1.997443109692506364316550909933324776725666683744229847556829344314537003376505], true)
a = 2.800222; b = 2.914870; r = 0.897481; x1 = 3.326237; y1 = 1.560762; x = 3.431978; y = 1.997443
2r = 1.794963
等円の直径は 1.794963 寸である。
算額の術では 1.830127 寸としているが，図に描く限り上の方が正しそう。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, r, x1, y1, x, y) = res[1]
@printf("a = %.6f; b = %.6f; r = %.6f; x1 = %.6f; y1 = %.6f; x = %.6f; y = %.6f\n", a, b, r, x1, y1, x, y)
@printf("2r = %.6f\n", 2r)
plot([0, a + b, 0, 0, x], [0, 0, 5, 0, y], color=:black, lw=0.5)
segment(0, 5, a, 0)
circle(r, y1, r, :blue)
circle(x1, r, r, :blue)
if more
point(x, y, " (x,y)", :green, :left, :bottom)
point(x1, r, " (x1,r)")
point(r, y1, " (r,y1)")
point(a, 0, "a ", :green, :right, :bottom)
point(a + b, 0, " a+b", :green, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
