算額（その1162）
八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，外円，二等辺三角形
#Julia, #SymPy, #算額, #和算
外円の中に二等辺三角形を容れ，甲円 3 個，乙円 4 個を容れる。乙円の直径が 2 寸のとき，甲円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
二等辺三角形の底辺の頂点の座標を (x, y); x = sqrt(R^2 - y^2)
甲円の半径と中心座標を r1, (0, y + r1), (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
@syms R::positive, x::positive, y::negative, r1::positive, x1::positive,
y1::positive, r2::positive, x2::positive,
y2::negative
@syms h, r1, r2, r3, x3, y31, y32, y33
eq1 = dist2(0, R, x, y, 0, y + r1, r1)
eq2 = dist2(0, R, x, y, x1, y1, r1)
eq3 = dist2(0, R, x, y, x2, y2, r2)
eq4 = x1^2 + y1^2 - (R - r1)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = (x2 - x1)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq7 = y1/x1 - x/(R - y)
eq8 = x^2 + y^2 - R^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)
(R, x, y, r1, x1, y1, x2, y2) = u
return [
-R^2*r1^2 + R^2*x^2 + 2*R*r1^2*y - 2*R*r1*x^2 - 2*R*x^2*y - r1^2*y^2 + 2*r1*x^2*y + x^2*y^2, # eq1
-R^2*r1^2 + R^2*x^2 - 2*R^2*x*x1 + R^2*x1^2 + 2*R*r1^2*y - 2*R*x^2*y1 + 2*R*x*x1*y + 2*R*x*x1*y1 - 2*R*x1^2*y - r1^2*x^2 - r1^2*y^2 + x^2*y1^2 - 2*x*x1*y*y1 + x1^2*y^2, # eq2
-R^2*r2^2 + R^2*x^2 - 2*R^2*x*x2 + R^2*x2^2 + 2*R*r2^2*y - 2*R*x^2*y2 + 2*R*x*x2*y + 2*R*x*x2*y2 - 2*R*x2^2*y - r2^2*x^2 - r2^2*y^2 + x^2*y2^2 - 2*x*x2*y*y2 + x2^2*y^2, # eq3
x1^2 + y1^2 - (R - r1)^2, # eq4
x2^2 + y2^2 - (R - r2)^2, # eq5
-(r1 + r2)^2 + (-x1 + x2)^2 + (y1 - y2)^2, # eq6
-x/(R - y) + y1/x1, # eq7
-R^2 + x^2 + y^2, # eq8
]
end;
r2 = 2/2
iniv = BigFloat[19, 10, -16.6, 7, 11.5, 3, 12, -8]
res = nls(H, ini=iniv)
([4.266666666666667, 2.0655911179772892, -3.7333333333333334, 1.6, 2.581988897471611, 0.6666666666666666, 2.6334969275741744, -1.9328230761165115], true)
甲円の直径は乙円の直径の 1.6 倍である。
乙円の直径が 2 のとき，甲円の直径は 3.2 である。
その他のパラメータは以下の通りである。
R = 4.26667; x = 2.06559; y = -3.73333; r1 = 1.6; x1 = 2.58199; y1 = 0.666667; r2 = 1; x2 = 2.6335; y2 = -1.93282
二等辺三角形の斜辺の長さ = 4.26667
function draw(r2, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, x, y, r1, x1, y1, x2, y2) = [19, 9.4, -16.6, 7, 11.5, 3, 12, -8]
(R, x, y, r1, x1, y1, x2, y2) = res[1]
@printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r2, 2r1)
@printf("R = %g; x = %g; y = %g; r1 = %g; x1 = %g; y1 = %g; r2 = %g; x2 = %g; y2 = %g\n", R, x, y, r1, x1, y1, r2, x2, y2)
@printf("二等辺三角形の斜辺の長さ = %g\n", sqrt(x^2 + y^2))
plot()
circle(0, 0, R)
circle2(x1, y1, r1, :blue)
circle(0, y + r1, r1, :blue)
circle2(x2, y2, r2, :green)
plot!([x, 0, -x, x], [y, R, y, y], color=:magenta, lw=0.5)
θ = atand(y1, x1) + atand(-y2, x2)
(x22, y22) = transform(x2, y2, deg = 2θ)
circle2(x22, y22, r2, :green)
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, y + r1, "甲円:r1,(0,y+r1)", :blue, :center, delta=-delta/2)
point(x1, y1, "甲円:r1,(x1,y1)", :blue, :center, delta=-delta/2)
point(x2, y2, "乙円:r2\n(x2,y2)", :green, :center, delta=-delta/2)
point(0, R, " R", :red, :left, :bottom, delta=delta/2)
point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
end
end;
---
甲円の半径 だけを求めるには以下のようにすれば，解析解を求めることができる。
a を二等辺三角形の斜辺の長さとする。
using SymPy
@syms a, r1, r2, R
eq1 = a^2 - 16r1*(R - r1)
eq2 = -4R*r1*(R - r1)/(R - 2r1) + a^2
eq3 = r2 - (R*r1 - r1^2)/R
res = solve([eq1, eq2, eq3], (r1, a, R))[2]
(8*r2/5, 32*sqrt(15)*r2/15, 64*r2/15)
res[1](r2 => 1).evalf() |> println
res[2](r2 => 1).evalf() |> println
res[3](r2 => 1).evalf() |> println
1.60000000000000
8.26236447190916
4.26666666666667
甲円の半径 r1 は，乙円の半径 r2 の 8/5 = 1.6 倍である。
乙円の直径が 2 寸のとき，甲円の直径は 3.2 寸である。
