算額(その914)
一〇八 加須市騎西町 玉敷神社 大正4年(1915)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
外円の中に甲円,乙円を 2 個ずつ入れ,甲円,乙円に外接する丙円を 2 個描く。甲円,乙円,外円の直径がそれぞれ 9 寸,3 寸,26 寸のとき,丙円の直径はいかほどか。
後述するが,上の図は外円,甲円,乙円の直径が 20, 6, 3 のときのもので,丙円の直径は 12.3435 である。
その他のパラメータは,R = 10; r1 = 3; r2 = 1.5; y1 = -6.32456; y2 = 8.3666; r3 = 6.17176; y3 = 2.28132
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (r3, y3)
とおき,以下の連立方程式を解く。
まとめて解くには SymPy の能力が足りない。
include("julia-source.txt");
using SymPy
@syms R::positive,
r1::positive, y1::negative,
r2::positive, y2::positive,
r3::positive, y3::positive
eq1 = r1^2 + y1^2 - (R - r1)^2 |> expand
eq2 = r2^2 + y2^2 - (R - r2)^2 |> expand
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;
外円,甲円,乙円の直径が与えられるので,eq1, eq2 はそれぞれだけで未知数 y1, y2 を求めることができる。
(ans_y1, ans_y2) = solve([eq1, eq2], (y1, y2))[1]
(-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))
y1, y2 が既知として,eq3, eq4 の連立方程式を解く。
(y1, y2) = (-sqrt(R)*sqrt(R - 2*r1), sqrt(R)*sqrt(R - 2*r2))
eq3 = (r3 - r1)^2 + (y3 - y1)^2 - (r1 + r3)^2 |> expand
eq4 = (r3 - r2)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand;
(ans_r3, ans_y3) = solve([eq3, eq4], (r3, y3))[1] # 1 of 2
((-R*r1 + R*r2 + (sqrt(R)*sqrt(R - 2*r1) + sqrt(R)*sqrt(R - 2*r2))*(-sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2)))/(2*(r1 - r2)), -sqrt(2)*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2 + sqrt(R - 2*r1)*sqrt(R - 2*r2))/(r1 - r2) + sqrt(R)*(r1*sqrt(R - 2*r2) + r2*sqrt(R - 2*r1))/(r1 - r2))
R, r1, r2 を与えれば,それぞれのパラメータの数値解が得られる。
甲円,乙円,外円の直径がそれぞれ 9 寸,3 寸,26 寸のとき,以下のようになる。
ans_y1(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y2(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_r3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
ans_y3(R => 13, r1 =>4.5, r2 => 1.5).evalf() |> println
-7.21110255092798
11.4017542509914
7.73565831477412
4.58897582448691
「術」では,「1385.4 を 89.57 で割れば丙円の直径が得られる」とあるが,「甲円,乙円,外円の直径がそれぞれ 9 寸,3 寸,26 寸のとき」の解であるし 1385.4/89.57 = 15.467232332254104 で,「答」の「丙円径十五寸四分二厘」とも小数点以下第2位で不一致である。
丙円の直径の正確な値は上に示すように 2*7.73565831477412 = 15.47131662954824 である。
任意の 甲円,乙円,外円の直径に対しては,上述の数式解に定数を代入すればよい。
y1, y2, r3, y3 の式は長く,SymPy では簡約化できないが,以下のようにすればプログラム的には短くなる。
s = √(R - 2r1)
t = √(R - 2r2)
u = -√(2R*r1*r2*(R - r1 - r2 + s*t))
v = √R*(r1*t + r2*s)
w = r1 - r2
(y1, y2, r3, y3) = (
-s√R,
t√R,
(s + t)*(u + v)√R/(2w^2) - R/2,
(u + v)/w
)
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
#(R, r1, r2) = (26, 9, 3) ./ 2
#(y1, y2, r3, y3) = (-7.21110255092798,
# 11.4017542509914,
# 7.73565831477412,
# 4.58897582448691)
(R, r1, r2) = (20, 6, 3) ./ 2
s = √(R - 2r1)
t = √(R - 2r2)
u = -√(2R*r1*r2*(R - r1 - r2 + s*t))
v = √R*(r1*t + r2*s)
w = r1 - r2
(y1, y2, r3, y3) = (
-s√R,
t√R,
(s + t)*(u + v)√R/(2w^2) - R/2,
(u + v)/w
)
@printf("外円,甲円,乙円の直径が %g, %g, %g のとき,丙円の直径は %g\n", 2R, 2r1, 2r2, 2r3)
@printf("R = %g; r1 = %g; r2 = %g; y1 = %g; y2 = %g; r3 = %g; y3 = %g\n", R, r1, r2, y1, y2, r3, y3)
plot()
circle(0, 0, R, :brown)
circle2(r1, y1, r1, :green)
circle2(r2, y2, r2)
circle2(r3, y3, r3, :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", :green, :left, :bottom, delta=delta/2)
point(r1, y1, "甲円:r1,(r1,y1)", :green, :center, delta=-delta/2)
point(r2, y2, "乙円:r2,(r2,y2)", :black, :center, :bottom, delta=delta/2)
point(r3, y3, "丙円:r3,(r3,y3)", :blue, :center, delta=-delta/2)
end
end;