算額(その1103)
九十九 江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード:円8個,外円,弦,斜線
全円の中に水平な弦と斜線を引き,甲円 3 個,乙円 2 個,丙円 2 個を容れる。丙円の直径が 8 寸のとき,乙円の直径はいかほどか。
全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (0, R - 3r1), (0, r1 - R)
乙円の半径と中心座標を r2,(x2, y2); y2 = -r1
丙円の半径と中心座標を r3, (x3, R - 2r1 + r3)
斜線と全円の交点座標を (x0, -sqrt(R^2 - x0^2))
とおき,以下の連立方程式を解く。
include("julia-source.txt")
using SymPy
@syms R::positive, r1::positive, r2::positive,
x2::positive, y2::positive, r3::positive,
x3::positive, x0::positive
y2 = -r1 # (R - 2r1) - (2R - 2r1)/2
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x3^2 + (R - 2r1 + r3)^2 - (R - r3)^2
eq4 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), 0, R - 3r1, r1) |> numerator
eq5 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), 0, r1 - R, r1) |> numerator
eq6 = dist2(sqrt(R^2 - (R - 2r1)^2), R - 2r1, -x0, -sqrt(R^2 - x0^2), x2, y2, r2) |> numerator;
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, r1, r2, x2, x3, x0) = u
s = sqrt(R^2 - x0^2)
t = sqrt(R - r1)
u = sqrt(r1)
v = r1^(3/2)
w = r1^(5/2)
return [
x3^2 + (r1 - r3)^2 - (r1 + r3)^2, # eq1
r1^2 + x2^2 - (R - r2)^2, # eq2
x3^2 - (R - r3)^2 + (R - 2*r1 + r3)^2, # eq3
8R^5*r1 - 42R^4*r1^2 + 8R^4*r1*s + 4R^3*v*x0*t + 78R^3*r1^3 - 42R^3*r1^2*s - 6R^3*r1*x0^2 - 26R^2*w*x0*t + 4R^2*v*x0*t*s - 46R^2*r1^4 + 78R^2*r1^3*s + 47R^2*r1^2*x0^2/2 - 2R^2*r1*x0^2*s + 60R*r1^(7/2)*x0*t - 26R*w*x0*t*s - 2R*v*x0^3*t - 78R*r1^4*s - 36R*r1^3*x0^2 + 5R*r1^2*x0^2*s/2 - 36*r1^(9/2)*x0*t + 20*r1^(7/2)*x0*t*s + w*x0^3*t + 36*r1^5*s + 20*r1^4*x0^2 - r1^3*x0^2*s, # eq4
12R^4*r1^2 + 4R^4*x0^2 + 24R^3*v*x0*t - 20R^3*r1^3 - 20R^3*r1^2*s + 8R^3*r1*x0^2 + 4R^3*x0^2*s - 52R^2*w*x0*t - 40R^2*v*x0*t*s + 16R^2*u*x0^3*t + 4R^2*r1^4 + 44R^2*r1^3*s - 73R^2*r1^2*x0^2 - 40R^2*r1*x0^2*s + 24R*r1^(7/2)*x0*t + 76R*w*x0*t*s - 60R*v*x0^3*t - 28R*r1^4*s + 88R*r1^3*x0^2 + 85R*r1^2*x0^2*s - 8*r1^(9/2)*x0*t - 24*r1^(7/2)*x0*t*s + 50*w*x0^3*t + 8*r1^5*s - 24*r1^4*x0^2 - 50*r1^3*x0^2*s, # eq5
4R^5*r1 - 4R^4*u*x0*t - 8R^4*u*x2*t - 12R^4*r1^2 + 4R^4*r1*s - 4R^4*r2^2 + R^4*x0^2 + 4R^4*x0*x2 + 4R^4*x2^2 + 24R^3*v*x0*t + 24R^3*v*x2*t - 4R^3*u*x0*t*s - 8R^3*u*x2*t*s + 28R^3*r1^3 - 20R^3*r1^2*s + 8R^3*r1*r2^2 - 6R^3*r1*x0^2 - 20R^3*r1*x0*x2 - 12R^3*r1*x2^2 - 4R^3*r2^2*s + R^3*x0^2*s + 4R^3*x0*x2*s + 4R^3*x2^2*s - 20R^2*w*x0*t - 32R^2*w*x2*t + 8R^2*v*x0*t*s + 24R^2*v*x2*t*s - 8R^2*u*r2^2*x0*t + 6R^2*u*x0^3*t + 12R^2*u*x0^2*x2*t + 4R^2*u*x0*x2^2*t - 20R^2*r1^4 + 20R^2*r1^3*s - 8R^2*r1^2*r2^2 + 21R^2*r1^2*x0^2 + 24R^2*r1^2*x0*x2 + 12R^2*r1^2*x2^2 + 8R^2*r1*r2^2*s - 16R^2*r1*x0^2*s - 20R^2*r1*x0*x2*s - 12R^2*r1*x2^2*s + 2R^2*r2^2*x0^2 - 2R^2*x0^3*x2 - 3R^2*x0^2*x2^2 + 8R*r1^(7/2)*x0*t - 28R*w*x0*t*s - 16R*w*x2*t*s - 24R*v*x0^3*t - 32R*v*x0^2*x2*t - 8R*v*x0*x2^2*t - 8R*u*r2^2*x0*t*s + 8R*u*x0^2*x2*t*s + 4R*u*x0*x2^2*t*s - 12R*r1^4*s - 40R*r1^3*x0^2 - 24R*r1^3*x0*x2 + 33R*r1^2*x0^2*s + 48R*r1^2*x0*x2*s + 12R*r1^2*x2^2*s - 16R*r1*r2^2*x0^2 + 14R*r1*x0^3*x2 + 8R*r1*x0^2*x2^2 - R*x0^2*x2^2*s - 8*r1^(9/2)*x0*t + 24*r1^(7/2)*x0*t*s + 16*r1^(7/2)*x2*t*s + 18*w*x0^3*t + 32*w*x0^2*x2*t + 8*w*x0*x2^2*t + 16*v*r2^2*x0*t*s - 12*v*x0^2*x2*t*s - 8*v*x0*x2^2*t*s - 2*u*x0^3*x2^2*t + 8*r1^5*s + 24*r1^4*x0^2 + 16*r1^4*x0*x2 - 18*r1^3*x0^2*s - 32*r1^3*x0*x2*s - 8*r1^3*x2^2*s + 16*r1^2*r2^2*x0^2 - 12*r1^2*x0^3*x2 - 8*r1^2*x0^2*x2^2 + 2*r1*x0^2*x2^2*s, # eq6
]
end;
r3 = 8/2
iniv = BigFloat[18.8, 5.8, 6.4, 10.8, 9.6, 11.8]
res = nls(H, ini=iniv)
([18.77777777777778, 5.777777777777778, 6.5, 10.833333333333334, 9.614803401237305, 11.786666666666667], true)
甲円の直径は乙円の直径の 13/8 倍である。
乙円の直径が 8 寸のとき,甲円の直径は 13 寸である。
その他のパラメータは以下のとおりである。
R = 169/9, r1 = 52/9, r2 = 13/2, x2 = 65/6, x3 = 9.614803401237305, x0 = 884/75
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = 8/2
(R, r1, r2, x2, x3, x0) = [169/9, 52/9, 13/2, 65/6, 9.614803401237305, 884/75]
(R, r1, r2, x2, x3, x0) = [18.77777777777778, 5.777777777777778, 6.5, 10.833333333333334, 9.614803401237305, 11.786666666666667]
y2 = (R - 2r1) - (2R - 2r1)/2
plot()
circle(0, 0, R, :blue)
circle(0, R - r1, r1, :magenta)
circle(0, R - 3r1, r1, :magenta)
circle(0, r1 -R, r1, :magenta)
circle2(x2, y2, r2)
circle2(x3, R - 2r1 + r3, r3, :orange)
segment(-sqrt(R^2 - (R - 2r1)^2), R - 2r1, sqrt(R^2 - (R - 2r1)^2), R - 2r1)
segment(x0, -sqrt(R^2 - x0^2), -sqrt(R^2 - (R - 2r1)^2), R - 2r1)
segment(-x0, -sqrt(R^2 - x0^2), sqrt(R^2 - (R - 2r1)^2), R - 2r1)
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", :blue, :left, :bottom, delta=delta/2)
point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
point(sqrt(R^2 - (R - 2r1)^2), R - 2r1, "(√(R^2-(R-2r1)^2),R-2r1)", :black, :right, delta=-delta/2)
point(-x0, -sqrt(R^2 - x0^2), "(-x0,-√(R^2-x0^2))", :black, :center, delta=-delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :magenta, :center, delta=-delta/2)
point(0, R - 3r1, "甲円:r1,(0,R-3r1)", :magenta, :center, delta=-delta/2)
point(0, r1 - R, "甲円:r1,(0,r1-R)", :magenta, :center, delta=-delta/2)
point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, delta=-delta/2)
point(x3, R - 2r1 + r3, "丙円:r3\n(x3,R-2r1+r3)", :black, :center, delta=-delta/2)
end
end;