算額（その1087）
二十六 一関市萩荘 赤萩観音寺 弘化4年(1847)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，正三角形，二等辺三角形
正三角形と二等辺三角形が重なってできる領域に等円を 3 個容れる。等円の直径が 1 寸のとき正三角形の一辺の長さはいかほどか。
正三角形の一辺の長さを 2a
等円の半径と中心座標を r, (x, r), (0, -r)
二等辺三角形の底辺の頂点座標を (x0, -2r)
とおいて以下の連立方程式を解く。
しかし，SymPy では一度に解けない。さらに，eq1, eq3 から x，x0 を求めて eq2 に代入し，方程式を解くも，a = √3r/3 という不適切回しかでてこない。方程式のグラフを描くと a = √3 という解があるのを確認できるのに不思議だ。
解析解を諦めて数値解を求めると解は求まる。
しかし，ここでめげてはいけない。初心に帰って図に示すような補助線を描くと，答えが見えてくる。
できる正三角形は全て合同で，一辺の長さは内接円の直径の √3 倍である。
したがって，等円の直径が 1 寸のとき，元の正三角形の一辺の長さは 2√3 = 3.4641016151377544 寸である。
include("julia-source.txt")
using SymPy
@syms a::positive, r::positive, x::positive, x0::positive
eq1 = dist2(a, 0, 0,√Sym(3)a, x, r, r)/√Sym(3)
eq2 = dist(x0, -2r, 0, √Sym(3)a, x, r)- r^2 |> expand |> simplify |> numerator |> factor
eq3 = x0*(√Sym(3)a + r) - r*sqrt(x0^2 + (√Sym(3)a + 2r)^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 Float64.(v), r.f_converged
end;
function H(u)
(a, x, x0) = u
return [
sqrt(3)*(3*a^2/4 - sqrt(3)*a*r/2 - 3*a*x/2 - 3*r^2/4 + sqrt(3)*r*x/2 + 3*x^2/4)/3, # eq1
-27*a^6*r^2 + 27*a^6*x^2 - 54*a^6*x*x0 + 27*a^6*x0^2 - 108*sqrt(3)*a^5*r^3 + 108*sqrt(3)*a^5*r*x^2 - 162*sqrt(3)*a^5*r*x*x0 + 54*sqrt(3)*a^5*r*x0^2 - 540*a^4*r^4 + 540*a^4*r^2*x^2 - 540*a^4*r^2*x*x0 + 54*a^4*r^2*x0^2 + 18*a^4*x^2*x0^2 - 36*a^4*x*x0^3 + 18*a^4*x0^4 - 480*sqrt(3)*a^3*r^5 + 480*sqrt(3)*a^3*r^3*x^2 - 240*sqrt(3)*a^3*r^3*x*x0 - 96*sqrt(3)*a^3*r^3*x0^2 + 48*sqrt(3)*a^3*r*x^2*x0^2 - 60*sqrt(3)*a^3*r*x*x0^3 + 12*sqrt(3)*a^3*r*x0^4 - 720*a^2*r^6 + 720*a^2*r^4*x^2 - 288*a^2*r^4*x0^2 + 144*a^2*r^2*x^2*x0^2 - 72*a^2*r^2*x*x0^3 - 27*a^2*r^2*x0^4 + 3*a^2*x^2*x0^4 - 6*a^2*x*x0^5 + 3*a^2*x0^6 - 192*sqrt(3)*a*r^7 + 192*sqrt(3)*a*r^5*x^2 + 96*sqrt(3)*a*r^5*x*x0 - 96*sqrt(3)*a*r^5*x0^2 + 64*sqrt(3)*a*r^3*x^2*x0^2 + 16*sqrt(3)*a*r^3*x*x0^3 - 20*sqrt(3)*a*r^3*x0^4 + 4*sqrt(3)*a*r*x^2*x0^4 - 2*sqrt(3)*a*r*x*x0^5 - 2*sqrt(3)*a*r*x0^6 - 64*r^8 + 64*r^6*x^2 + 64*r^6*x*x0 - 32*r^6*x0^2 + 32*r^4*x^2*x0^2 + 32*r^4*x*x0^3 - 4*r^4*x0^4 + 4*r^2*x^2*x0^4 + 4*r^2*x*x0^5, # eq2
-r*sqrt(x0^2 + (sqrt(3)*a + 2*r)^2) + x0*(sqrt(3)*a + r), # eq3
]
end;
r = 1/2
iniv = BigFloat[1.7, 0.87, 0.58]
res = nls(H, ini=iniv)
([1.7320508075688774, 0.8660254037844387, 0.5773502691896257], true)
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 1/2
(a, x, x0) = [1.7320508075688774, 0.8660254037844387, 0.5773502691896257]
#@printf("小円の直径が %g のとき，中円の直径は %g である。\n", 2r4, 2r3)
#@printf("r4 = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g\n", r4, r1, r2, r3, x3, y3)
plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
plot!([x0, 0, -x0, x0], [-2r, √3a, -2r, -2r], color=:green, lw=0.5)
plot!([0, x, -x, 0], [0, √3x, √3x, 0], color=:blue, lw=1)
circle2(x, r, r)
circle(0, -r, r)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(x, r, "等円:r,(x,r)", :red, :center, delta=-delta/2)
point(0, -r, "等円:r,(0,-r)", :red, :center, delta=-delta/2)
point(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
point(x0, -2r, " (x0,-2r)", :green, :left, :vcenter)
end
end;
