算額（その756）
山形県鶴岡市茨新田 鶴岡神明宮 文久2年(1862)
http://www.wasan.jp/yamagata/sinmei.html
底辺を共有し交わる二等辺三角形と直角三角形がある。二等辺三角形に内接する全円と，区分された領域に甲円乙円，丙円が入っている。
全円，甲円，乙円の直径がそれぞれ 5 寸，2 寸，4 寸のとき，丙円の直径はいかほどか。
底辺の長さを 2a，二等辺三角形と直角三角形の高さを h, h2 とする
全円の半径と中心座標を r1, (a, r1)
甲円の半径と中心座標を r2, (a, y2)
乙円の半径と中心座標を r3, (x3, r3)
丙円の半径と中心座標を r4, (2a - r4, y4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, h::positive, h2::positive,
r1::positive, r2::positive, y2::positive,
r3::positive, x3::negative,
r4::positive, y4::positive
(r1, r2, r3) = (5, 2, 4) .// 2
eq1 = dist(a, h, 2a, 0, a, r1) - r1^2
eq2 = dist(a, h, 2a, 0, a, y2) - r2^2
eq3 = dist(a, h, 2a, 0, x3, r3) - r3^2
eq4 = dist(a, h, 2a, 0, 2a - r4, y4) - r4^2
eq5 = dist(0, 0, 2a, h2, a, y2) - r2^2
eq6 = dist(0, 0, 2a, h2, x3, r3) - r3^2
eq7 = dist(0, 0, 2a, h2, 2a - r4, y4) - r4^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, h, h2, y2, x3, r4, y4) = u
return [
a^2*h^2*(5/2 - h)^2/(a^2 + h^2)^2 + (-h^2*(5/2 - h)/(a^2 + h^2) - h + 5/2)^2 - 25/4, # eq1
a^2*h^2*(-h + y2)^2/(a^2 + h^2)^2 + (-h^2*(-h + y2)/(a^2 + h^2) - h + y2)^2 - 1, # eq2
(-a - a*(a*(-a + x3) - h*(2 - h))/(a^2 + h^2) + x3)^2 + (-h + h*(a*(-a + x3) - h*(2 - h))/(a^2 + h^2) + 2)^2 - 4, # eq3
-r4^2 + (a - a*(a*(a - r4) - h*(-h + y4))/(a^2 + h^2) - r4)^2 + (-h + h*(a*(a - r4) - h*(-h + y4))/(a^2 + h^2) + y4)^2, # eq4
(-2*a*(2*a^2 + h2*y2)/(4*a^2 + h2^2) + a)^2 + (-h2*(2*a^2 + h2*y2)/(4*a^2 + h2^2) + y2)^2 - 1, # eq5
(-2*a*(2*a*x3 + 2*h2)/(4*a^2 + h2^2) + x3)^2 + (-h2*(2*a*x3 + 2*h2)/(4*a^2 + h2^2) + 2)^2 - 4, # eq6
-r4^2 + (-h2*(2*a*(2*a - r4) + h2*y4)/(4*a^2 + h2^2) + y4)^2 + (2*a - 2*a*(2*a*(2*a - r4) + h2*y4)/(4*a^2 + h2^2) - r4)^2, # eq7
]
end;
(r1, r2, r3) = (5, 2, 4) ./ 2
iniv = BigFloat[4.0, 7, 8, 5, 5, 1, 5]
res = nls(H, ini=iniv)
([5.0, 6.666666666666667, 7.5, 5.0, 6.0, 1.5, 4.5], true)
丙円の直径は 3 寸である。
その他のパラメータは以下のとおり。
r1 = 2.5; r2 = 1; r3 = 2; a = 5; h = 6.66667; h2 = 7.5; y2 = 5; x3 = 6; r4 = 1.5; y4 = 4.5
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
# (r1, r2, r3) = (5, 2, 4) ./ 2
(a, h, h2, y2, x3, r4 , y4) = res[1]
@printf("丙円の直径 = %g; 全円の直径 = %g; 甲円の直径 = %g; 乙円の直径 = %g\n", 2r4, 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; r3 = %g; a = %g; h = %g; h2 = %g; y2 = %g; x3 = %g; r4 = %g; y4 = %g\n", r1, r2, r3, a, h, h2, y2, x3, r4, y4)
plot([0, 2a, a, 0], [0, 0, h, 0], color=:black, lw=0.5)
plot!([0, 2a, 2a, 0], [0, 0, h2, 0], color=:black, lw=0.5)
circle(a, r1, r1)
circle(a, y2, r2, :blue)
circle(x3, r3, r3, :green)
circle(2a - r4, y4, r4, :magenta)
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(a, r1, "全円:r1,(a,r1)", :red, :center, :bottom, delta=delta)
point(a, y2, "甲円:r2\n(a,y2)", :blue, :center, :bottom, delta=delta)
point(x3, r3, "乙円:r3,(x3,r3)", :green, :center, delta=-delta)
point(2a - r4, y4, "丙円:r4,(2a-r4,y4)", :magenta, :center, delta=-delta)
point(2a, 0, " 2a", :black, :left, :bottom, delta=delta)
point(a, h, "(a,h)", :black, :center, :bottom, delta=delta)
point(2a, h2, "(2a,h2)", :black, :center, :bottom, delta=delta)
plot!(xlims=(-0.2, 10.5))
end
end;
