算額(その1001)
番外六 広野村川嶋(現嵐山町) 鬼神社
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
外円(半円)の中に斜線 2 本と,甲円(半円),乙円,丙円を入れる。乙円の直径が与えられたときに,丙円の直径を求めよ。
外円の半径と中心座標を R, (0, 0); R = 2r1
甲円の半径と中心座標を r1, (r1, 0)
乙円の半径と中心座標を r2, (R - r2)
丙円の半径と中心座標を r3, (x3, y3)
とおき,以下の連立法て式を解く。
SymPy の能力的に,一度には解けないので,丙円は後回しにして,「r1 を未知数として」外円,甲円,乙円,および斜線のパラメータを決定する。
include("julia-source.txt");
using SymPy
@syms R::poitive, r1::poitive, r2::poitive, r3::poitive,
x3::poitive, y3::poitive, x0::poitive, y0::poitive
R = 2r1
#y0 = (x0 + R)/3
eq3 = x0^2 + y0^2 - R^2
eq4 = dist2(-R, 0, x0, y0, 0, R - r2, r2)
eq6 = dist2(-R, 0, x0, y0, r1, 0, r1);
#= 以下の 3 本の式は丙円に関するもの
eq1 = (x3 - r1)^2 + y3^2 - (r1 - r3)^2
eq2 = y3/(x3 - r1) - (x0 + R)/y0
eq5 = dist2(-R, 0, x0, y0, -x3, y3, r3);
=#
# solve([eq1, eq2, eq3, eq4, eq5], (r1, r3, x3, y3, x0))
res1 = solve([eq3, eq4, eq6], (r2, x0, y0))[5] # 5 of 5
(-65*r1/8 + 9*sqrt(2)*r1/4 + 3*sqrt(2)*sqrt(r1*(65*r1 + 63*sqrt(r1^2)))*sqrt(9 - 4*sqrt(2))/8 - 63*sqrt(r1^2)/8 + 7*sqrt(2)*sqrt(r1^2)/4, 14*sqrt(r1^2)/9, 8*sqrt(2)*r1/9)
5 組の解が得られるが,5 番目のものが適解である。更に得られる r2 の式は非常に短い形に簡約化できる。
次に,r2, x0, y0 が既知というスタート地点に立って,eq1, eq2, eq5 を解いて丙円の半径 r3 と中心座標 9x3, y3) を決定する。
@syms R::poitive, r1::poitive, r2::poitive, r3::poitive,
x3::poitive, y3::poitive, x0::poitive, y0::poitive
R = 2r1
r2 = 2r1*(8sqrt(Sym(2)) - 11) # r2 は超簡約化できる
x0 = 14r1/9
y0 = 8sqrt(Sym(2))*r1/9
eq1 = (x3 - r1)^2 + y3^2 - (r1 - r3)^2
eq2 = y3/(x3 - r1) - (x0 + R)/y0
eq5 = dist2(-R, 0, x0, y0, -x3, y3, r3);
res2 = solve([eq1, eq2, eq5], (r3, x3, y3))[1] # 1 of 2
(r1/3, 11*r1/9, 4*sqrt(2)*r1/9)
2 組の解が得られるが,最初のものが適解である。
さて,算額の問は,「乙円の直径がわかっているときに,丙円の直径を求めよということである。
連立方程式を解いて,乙円と丙円の半径 r2, r3 が甲円の半径 r1 でどのように表されるかがわかったので,r2から直接 r3 を求める式 r3/r2 を作ればよい。
@syms r1, r2, r3
r2 = 2r1*(8sqrt(Sym(2)) - 11)
r3 = r1/3
r3/r2 |> simplify |> factor |> println
(11 + 8*sqrt(2))/42
丙円の直径は,乙円の直径の (11 + 8√2)/42 倍である。
「術」は「置十八箇開平方ニ以減六箇餘以除乙円径得丙円径」と書いてあるが,実際の計算方法は上で述べたものとは異なるようだ。正しいかどうかは筆者にはわからない。
r2 = 3 のときの数値解を求めると r3 = 1.5938363213560542 となり,3 * (11 + 8√2)/42 = 1.593836321356054 と一致する。
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, r3, x3, y3, x0, y0) = u
return [
y3^2 + (-r1 + x3)^2 - (r1 - r3)^2, # eq1
y3/(-r1 + x3) - (2*r1 + x0)/y0, # eq2
-4*r1^2 + x0^2 + y0^2, # eq3
16*r1^4 - 16*r1^3*r2 + 16*r1^3*x0 - 16*r1^3*y0 - 16*r1^2*r2*x0 + 8*r1^2*r2*y0 + 4*r1^2*x0^2 - 8*r1^2*x0*y0 + 4*r1^2*y0^2 - 4*r1*r2*x0^2 + 4*r1*r2*x0*y0 - r2^2*y0^2, # eq4
-4*r1^2*r3^2 + 4*r1^2*y0^2 - 8*r1^2*y0*y3 + 4*r1^2*y3^2 - 4*r1*r3^2*x0 - 4*r1*x0*y0*y3 + 4*r1*x0*y3^2 - 4*r1*x3*y0^2 + 4*r1*x3*y0*y3 - r3^2*x0^2 - r3^2*y0^2 + x0^2*y3^2 + 2*x0*x3*y0*y3 + x3^2*y0^2, # eq5
r1^2*(-4*r1^2 - 4*r1*x0 - x0^2 + 8*y0^2), # eq6
]
end;
r2 = 3
iniv = BigFloat[25/2, 4, 15, 8, 20, 11]
res = nls(H, ini=iniv)
(r1, r3, x3, y3, x0, y0) = res[1]
res |> println
r3/r2 |> println
([4.781508964068163, 1.5938363213560542, 5.844066511638866, 3.005366589152766, 7.43790283299492, 6.010733178305532], true)
0.5312787737853514
r1 = 1.234 のとき,それぞれのパラメータは以下のとおりである。
r1 = 1.234; R = 2.468; r2 = 0.774233; x0 = 1.91956; y0 = 1.55124; r3 = 0.411333; x3 = 1.50822; y3 = 0.775618
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1.234
R = 2r1
(r2, x0, y0) = r1 .* (16√2 - 22, 14/9, 8√2/9)
(r3, x3, y3) = r1 .* (1/3, 11/9, 4√2/9)
@printf("r1 = %g; R = %g; r2 = %g; x0 = %g; y0 = %g; r3 = %g; x3 = %g; y3 = %g\n", r1, R, r2, x0, y0, r3, x3, y3)
plot()
circle(0, 0, R, :blue, beginangle=0, endangle=180)
circle(r1, 0, r1, :green, beginangle=0, endangle=180)
circle(-r1, 0, r1, :green, beginangle=0, endangle=180)
circle2(x3, y3, r3, :magenta)
circle(0, R - r2, r2, :orange)
segment(-R, 0, x0, y0, :gray)
segment(R, 0, -x0, y0, :gray)
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(r1, 0, "甲円:r1,(r1,0)", :green, :center, :bottom, delta=delta)
point(0, R - r2, "乙円:r2,(0,R-r2)", :orange, :center, delta=-delta/2)
point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, :bottom, delta=delta/2)
point(x0, y0, " (x0,y0)", :gray, :left, :vcenter)
end
end;