算額（その1331）
四十 岩手県一関市牧沢 牧沢八幡神社 明治5年(1872)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円5個，外円，斜線2本
#Julia, #SymPy, #算額, #和算
大円（外円）の中に小円 1 個と中円 3 個を容れる。小円と大円の直径が与えられたとき，中円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (0, R - r1), (x1, y1)
小円の半径と中心座標を r2, (0, r2 - R)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive,
x1::positive, y1::positive,
r2::positive, x0::positive, y0::positive;
eq1 = dist2(0, R, x0, -sqrt(R^2 - x0^2), x1, y1, r1)
eq2 = dist2(0, R, x0, -sqrt(R^2 - x0^2), 0, r2 - R, r2)
eq3 = x1^2 + y1^2 - (R - r1)^2
eq4 = x1^2 + (R - r1 - y1)^2 - 4r1^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)
(r1, x1, y1, x0) = u
t = sqrt(R^2 - x0^2)
return [
(-4*R^3*r1^2 + R^3*x0^2 - 4*R^3*x0*x1 + 4*R^3*x1^2 - 4*R^2*r1^2*sqrt(R^2 - x0^2) - 2*R^2*x0^2*y1 + R^2*x0^2*sqrt(R^2 - x0^2) + 4*R^2*x0*x1*y1 - 4*R^2*x0*x1*sqrt(R^2 - x0^2) + 4*R^2*x1^2*sqrt(R^2 - x0^2) + 2*R*r1^2*x0^2 + 2*R*x0^3*x1 - 3*R*x0^2*x1^2 + R*x0^2*y1^2 - 2*R*x0^2*y1*sqrt(R^2 - x0^2) + 4*R*x0*x1*y1*sqrt(R^2 - x0^2) - 2*x0^3*x1*y1 - x0^2*x1^2*sqrt(R^2 - x0^2) + x0^2*y1^2*sqrt(R^2 - x0^2))/(2*R*(2*R^2 + 2*R*sqrt(R^2 - x0^2) - x0^2)), # eq1
(-4*R^3*r2^2 + 4*R^3*x0^2 - 4*R^2*r2^2*sqrt(R^2 - x0^2) - 4*R^2*r2*x0^2 + 4*R^2*x0^2*sqrt(R^2 - x0^2) + 3*R*r2^2*x0^2 - 4*R*r2*x0^2*sqrt(R^2 - x0^2) + r2^2*x0^2*sqrt(R^2 - x0^2))/(2*R*(2*R^2 + 2*R*sqrt(R^2 - x0^2) - x0^2)), # eq2
x1^2 + y1^2 - (R - r1)^2, # eq3
-4*r1^2 + x1^2 + (R - r1 - y1)^2, # eq4
]
end;
(R, r2) = (1/2, 0.2/2)
iniv = BigFloat[0.2, 0.3, 0.06, 0.2]
res = nls(H, ini=iniv)
res |> println
2res[1][1]|> println
([0.21665006633518888, 0.2792628680296955, -0.04795242900710796, 0.11042310999998962], true)
0.43330013267037776
大円と小円の直径が 1 寸, 0.2 寸のとき，中円の直径は 0.43330013267037776 寸である。
function draw(R, r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, y1, x0) = res[1]
y0 = -sqrt(R^2 - x0^2)
@printf("大円と小円の直径が %g, %g のとき，中円の直径は %g である。\n", 2R, 2r2, 2r1)
plot([-x0, 0, x0], [y0, R, y0], color=:green, lw=0.5)
circle(0, 0, R, :magenta)
circle(0, R - r1, r1)
circle2(x1, y1, r1)
circle(0, r2 - R, r2, :blue)
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(0, R, "R", :magenta, :center, :bottom, delta=delta/2)
point(0, R - r1, "中円:r1\n(0,R-r1)", :red, :left, :vcenter, deltax=3delta)
point(x1, y1, "中円:r1,(x1,y1)", :red, :center, delta=-delta/2)
point(0, r2 - R, "小円:r2\n(0,r2-R)", :blue, :center, delta=-delta/2)
point(x0, -√(R^2 - x0^2), " -(x0,√(R^2-x0^2))", :green, :left, :vcenter)
end
end;
draw(1/2, 0.2/2, true)
