算額（その381）
静岡県伊豆市 江川邸 享和2年(1802)
http://www.wasan.jp/sizuoka/egawa.html
長野県長野市若穂 清水寺観音堂 寛政6年(1794)
中村信弥「改訂増補 長野県の算額」(45)
http://www.wasan.jp/zoho/zoho.html
正三角形内に 3 本の斜線を引き，分割された領域に 4 個の等円を内接させる。正三角形の一辺の長さがが与えられるとき，等円の直径はいかほどか。
正三角形の一辺の長さを a，等円の半径を r とする。
以下の連立方程式を解けばよいが，solve() では無理のようなので，nlsove() で数値解を求める。解析解は後半に記載。
include("julia-source.txt");
using SymPy
@syms a::positive, r::positive, x::positive, y::positive,
x1::positive,y1::positive;
eq1 = distance(0, √Sym(3)a/2, a/2, 0, x, y) - r^2
eq2 = distance(0, √Sym(3)a/2, x1, y1, x, y) - r^2
eq3 = distance(0, √Sym(3)a/2, x1, y1, 0, a/2√Sym(3)) - r^2
eq4 = distance(x1, y1, a/2, 0, x, y) - r^2
eq5 = distance(x1, y1, a/2, 0, 0, a/2√Sym(3)) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5])
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, x1, y1) = u
return [
-r^2 + (-3*a/8 + 3*x/4 + sqrt(3)*y/4)^2 + (-sqrt(3)*a/8 + sqrt(3)*x/4 + y/4)^2, # eq1
-r^2 + (x - x1*(3*a^2 - 2*sqrt(3)*a*y - 2*sqrt(3)*a*y1 + 4*x*x1 + 4*y*y1)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2 + (y - (3*a^2*y - 2*sqrt(3)*a*x*x1 + 2*sqrt(3)*a*x1^2 - 4*sqrt(3)*a*y*y1 + 4*x*x1*y1 + 4*y*y1^2)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2, # eq2
4*a^2*x1^2*(3*a - 2*sqrt(3)*y1)^2/(9*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (sqrt(3)*a/6 - sqrt(3)*a*(3*a^2 - 4*sqrt(3)*a*y1 + 12*x1^2 + 4*y1^2)/(6*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)))^2, # eq3
-r^2 + (x - (a^2*x - 4*a*x*x1 - 2*a*y*y1 + 2*a*y1^2 + 4*x*x1^2 + 4*x1*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2 + (y - y1*(a^2 - 2*a*x - 2*a*x1 + 4*x*x1 + 4*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2, # eq4
a^2*y1^2*(-sqrt(3)*a + 2*sqrt(3)*x1 + 6*y1)^2/(9*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (-a*y1*(3*a - 6*x1 + 2*sqrt(3)*y1)/(3*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)) + sqrt(3)*a/6)^2, # eq5
]
end;
a = 10
iniv = [big"0.09", 0.27, 0.19, 0.17, 0.09] .* a
res = nls(H, ini=iniv);
println(res);
(BigFloat[1.142125629369641587973968514127294180492236375359075288930728366920198239037438, 2.50000000000000007079306505041260465903195757502042722173815452658909699459355, 2.045875760182909850411741235498302045060062420025658180702662505806048892327303, 1.510890190652600068372820337416072773295486747795635718444765203143264858414334, 1.173562901893666266337957465962928304869915646101974796747607212272768807266053], true)
r = 1.14213; a = 10; x = 2.5; y = 2.04588; x1 = 1.51089; y1 = 1.17356
小円の直径 = 2r = 2.28425
解析解
正三角形の面積を 3 個の相似な三角形と中央の小さな正三角形の面積に分割する。
中央の正三角形の面積は 3√3r^2 である（注）。
下側の三角形に注目する。鈍角は 120°であるので，3 辺の和は 2a + 2*r/√3 で，面積は (2a + 2*r/√3)r/2である。
よって，
√3a^2/4 = 3(2a + 2*r/√3)r/2 + 3√3r^2
を r について解けばよい。
注: 正三角形の内接円の半径が r のとき，正三角形の面積は 3√3r^2 である。
@syms a::positive, r::positive
eq = 3(2a + 2*r/√Sym(3))r/2 + 3√Sym(3)r^2 - √Sym(3)a^2/4 |> simplify
eq |> println
-sqrt(3)*a^2/4 + 3*a*r + 4*sqrt(3)*r^2
solve(eq, r)[1] |> println
a*(-sqrt(3) + sqrt(7))/8
直径は (√7 - √3)/4 * a である。
a = 10 のとき 2.2842512587392836 である。
(√7 - √3)/4 * 10
2.2842512587392836
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 10
(r, x, y, x1, y1) = res[1]
@printf("r = %g; a = %g; x = %g; y = %g; x1 = %g; y1 = %g\n", r, a, x, y, x1, y1)
@printf("小円の直径 = 2r = %g\n", 2r)
plot([a/2, 0, -a/2, a/2], [0, √3a/2, 0, 0], color=:black, linewidth=0.25)
circle(0, a/2√3, r)
for deg = 0:120:240
A = [cosd(deg) -sind(deg); sind(deg) cosd(deg)]
(x2, y2, p1, p2, p3, p4) = A * [x 0 x1; y-a/2√3 √3a/2-a/2√3 y1-a/2√3]
segment(p1, p2 + a/2√3, p3, p4 + a/2√3, :blue)
circle(x2, y2 + a/2√3, r)
end
if more == true
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(0, a/2√3, " a/2√3", :red)
point(x, y, "(x,y)", :red)
point(x1, y1, "(x1,y1) ", :black, :right)
point(0, √3a/2, " √3a/2", :black)
point(a/2, 0, " a/2", :black, :left, :bottom)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
