算額（その1322）
一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，斜線2本
#Julia, #SymPy, #算額, #和算
直線上に，2 本の接線を共有する大円 3 個，中円 1 個，小円 1 個を描く。大円と小円の直径の差が 2 寸のとき，中円の直径はいかほどか。
上の図は，大円と小円の直径の差が 5 寸のときのものである。問題の通り大円と小円の直径の差が 2 寸のときのものは下の方に示すが，似ても似つかない図になる。
大円と小円の直径の差を K
大円の半径と中心座標を r1, (x1, r1), (0, y1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (0, y1 - r1 - r3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, x1::positive, y1::positive,
r2::positive, r3::positive, a::positive,
K::positive;
eq1 = dist2(-a, 0, 0, r1, x1, r1, r1)
eq2 = dist2(-a, 0, 0, r1, 0, y1, r1)
eq3 = dist2(-a, 0, 0, r1, 0, r2, r2)
eq4 = dist2(-a, 0, 0, r1, 0, y1 - r1 - r3, r3)
eq5 = r2*(a + x1) - r1*a # r2/a - r1/(a + x1)
eq6 = (2r1 - 2r3) - K;
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, r2, r3, a) = u
return [
r1^2*(-a^2 - r1^2 + x1^2), # eq1
-2*a^2*r1*y1 + a^2*y1^2 - r1^4, # eq2
r1*(a^2*r1 - 2*a^2*r2 - r1*r2^2), # eq3
4*a^2*r1^2 + 4*a^2*r1*r3 - 4*a^2*r1*y1 - 2*a^2*r3*y1 + a^2*y1^2 - r1^2*r3^2, # eq4
-a*r1 + r2*(a + x1), # eq5
-K + 2*r1 - 2*r3, # eq6
]
end;
K = 2.0
iniv = BigFloat[1.5, 6.3, 3.1, 0.75, 0.023, 6.1]
#iniv = BigFloat[3, 4.3, 7.2, 1.25, 0.5, 3.1] # K = 5 のときの初期値
res = nls(H, ini=iniv)
([1.0034409020547763, 8.597223980953418, 2.0137872878330065, 0.5, 0.003440902054776372, 8.538463944689585], true)
大円と小円円の直径の差が 2 寸のとき，中円の直径は 1 寸である。
全てのパラメータは以下のとおりである。
K = 2; r1 = 1.00296; x1 = 9.25092; y1 = 2.01188; r2 = 0.5; r3 = 0.00296477; a = 9.19639
function draw(K, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, x1, y1, r2, r3, a) = res[1]
@printf("大円と小円円の直径の差が %g のとき，中円の直径は %g である。\n", K, 2r2)
@printf("K = %g; r1 = %g; x1 = %g; y1 = %g; r2 = %g; r3 = %g; a = %g\n", K, r1, x1, y1, r2, r3, a)
plot()
circle2(x1, r1, r1)
circle(0, y1, r1)
circle(0, r2, r2, :blue)
circle(0, y1 - r1 - r3, r3, :green)
abline(-a, 0, r1/a, -a, 1.2x1, :magenta)
abline(a, 0, -r1/a, a, -1.2x1, :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(0, y1, "大円:r1,(0,y1)", :red, :center, delta=-delta/2)
point(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
point(0, r2, " 中円:r2\n(0,r2)", :blue, :center, delta=-delta/2)
point(0, y1 - r1 - r3, " 小円:r3,(0,y1-r1-r3)", :black, :left, :vcenter)
end
end;
draw(2, false)
