算額(その755)
宮城県白石市小原字馬頭山 三瀧神社奉納算額(萬蔵稲荷神社所蔵) 明治8年
徳竹亜紀子,谷垣美保:宮城県白石市小原地区の算額調査,仙台高等専門学校名取キャンパス研究紀要,第57号,2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf
外円の中に,甲円,乙円,大円,中円,小円が入っている。外円,甲円,乙円の直径がそれぞれ 64 寸,52 寸 9 分,36 寸 1 分のとき,中円の直径はいかほどか。
中円は甲円に内接し,乙円と外接している。小円は甲円と外接し,乙円に内接している。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
大円の半径と中心座標を r3, (0, r1 - r2)
中円の半径と中心座標を r4, (r4, y4)
小円の半径と中心座標を r5, (r5, y5)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive,
r4::positive, y4::negative, r5::positive, y5::positive
#(R, r1, r2) = (640, 529, 361) .// 20
eq1 = r4^2 + (r1 - R - y4)^2 - (r1 - r4)^2
eq2 = r5^2 + (y5 - R + r2)^2 - (r2 - r5)^2
eq3 = r4^2 + (R - r2 - y4)^2 - (r2 + r4)^2
eq4 = r5^2 + (y5 - r1 + R)^2 - (r1 + r5)^2
eq5 = 2r1 + 2r2 - 2r3 - 2R
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, r4, y4, r5, y5))
4-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
(-R + r1 + r2, (R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
(-R + r1 + r2, (R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
(-R + r1 + r2, (R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
(-R + r1 + r2, (R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
for i = 1:4
print("#$i: ")
for j = 1:5
print((res[i][j](R => 64.0/2, r1 => 52.9/2, r2 => 36.1/2)).evalf(), "; ")
end
println()
end
#1: 12.5000000000000; 12.0000000000000; -13.6000000000000; -12.0000000000000; -13.6000000000000;
#2: 12.5000000000000; 12.0000000000000; -13.6000000000000; 5.21297815932332; 25.6808988764045;
#3: 12.5000000000000; -5.21297815932332; 25.6808988764045; -12.0000000000000; -13.6000000000000;
#4: 12.5000000000000; -5.21297815932332; 25.6808988764045; 5.21297815932332; 25.6808988764045;
4 組の解が得られるが,2 番目のものが適解である。
中円の半径を表す数式はかなり長いものであり,SymPy ではうまく簡約化できない。
res[2][2] |> println
(R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2)
外円,甲円,乙円の直径がそれぞれ 64 寸,52 寸 9 分,36 寸 1 分のとき,中円の直径は 24 寸である。
その他のパラメータは以下のとおり。
外円の直径 = 64; 甲円の直径 = 52.9; 乙円の直径 = 36.1
大円の直径 = 25; 小円の直径 = 10.426; r3 = 12.5; r4 = 12; y4 = -13.6; r5 = 5.21298; y5 = 25.6809
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r2) = (640, 529, 361) .// 20
# (r3, r4, y4, r5, y5) = (25/2, 12, -68/5, 41292/7921, 11428/445)
t0 = r1 + r2
t1 = 2sqrt(R*r1*r2*(-R + t0))/t0
t2 = R*(r1 - r2)/t0
(r3, r4, y4, r5, y5) = (
t0 - R,
(R*(r1 - r2) + (t2 - t1)*(t0 - 2R))/t0,
t2 - t1,
-(R*(r1 - r2) + (t2 + t1)*(t0 - 2R))/t0,
t2 + t1)
@printf("外円の直径 = %g; 甲円の直径 = %g; 乙円の直径 = %g\n", 2R, 2r1, 2r2)
@printf("中円の直径 = %g; 大円の直径 = %g; 小円の直径 = %g; r3 = %g; r4 = %g; y4 = %g; r5 = %g; y5 = %g\n",
2r4, 2r3, 2r5, r3, r4, y4, r5, y5)
plot()
circle(0, 0, R, :black)
circle(0, r1 - R, r1)
circle(0, R - r2, r2, :blue)
circle(0, r1 - r2, r3, :green)
circle2(r4, y4, r4, :magenta)
circle2(r5, y5, r5, :orange)
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, r1 - R, " 甲円:r1,(0,r1-R)", :red, :left, :vcenter)
point(0, R - r2, "乙円:r2 \n(0,R-r2) ", :blue, :right, :vcenter)
point(0, r1 - r2, " 大円:r3\n (0,r1-r2)", :green, :left, :vcenter)
point(r4, y4, "中円:r4,(r4,y4)", :magenta, :center, delta=-delta)
point(r5, y5, "小円:r5\n(r5,y5)", :orange, :center)
point(R, 0, " R", :black, :left, :bottom, delta=delta/2)
end
end;