算額(その447)
河□狭山牛頭天王社 林助右衛門自弘 寛政7年乙卯11月(1795)
藤田貞資(1789):神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
図のように 12 個の円がある。甲円,乙円,丙円の直径がそれぞれ 43寸7分5厘,51寸0分3厘,141寸7分5厘のとき丁円の直径はいかほどか。
![](https://blogimg.goo.ne.jp/user_image/79/86/96dc4226260693061fea9d8f815e8643.png)
左右対称なので,右半分の図形で方程式を立てる。
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (r3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (r5, y5)
己円の半径と中心座標を r6, (r6, 0)
とおき,連立方程式の解を求める。
include("julia-source.txt")
using SymPy
@syms r1::positive, y1::positive,
r2::positive, y2::positive,
r3::positive, y3::positive,
r4::positive, x4::positive, y4::positive,
r5::positive, y5::positive,
r6::positive;
(r1, r2, r3) = (4375, 5103, 14175) .// 200
eq1 = (r1 - r3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (r1 - r5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq3 = (r2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = (r6 - r2)^2 + y2^2 - (r2 + r6)^2
eq5 = (r3 - r5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq6 = (r3 - r6)^2 + y3^2 - (r3 + r6)^2
eq7 = (x4 - r5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq8 = (r3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq9 = (x4 - r6)^2 + y4^2 - (r4 + r6)^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (y1, y2, y3, r4, x4, y4, r5, y5, r6))
3-element Vector{NTuple{9, Sym}}:
(1071/8, 5103/40, 1701/8, 189/920, 84861/920, 26649/184, 14175/128, 567/16, 5103/32)
(1071/8, 5103/40, 1701/8, 6048/215, 18333/860, 43659/344, 2025/224, 162, 5103/32)
(2331/8, 5103/40, 1701/8, 3267/40, 8883/40, 1863/8, 14175/128, 6237/16, 5103/32)
3 組目のものが適解である。
r1 = 21.875; r2 = 25.515; r3 = 70.875
y1 = 291.375; y2 = 127.575; y3 = 212.625; r4 = 81.675; x4 = 222.075; y4 = 232.875; r5 = 110.742; y5 = 389.812; y6 = 159.469
丁円の直径は 2*3267/40 = 163.35, 163寸3分5厘である。
2*3267/40
163.35
using Plots
function circle2(x, y, r, color)
circle(x, y, r, color)
circle(-x, y, r, color)
end
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (4375, 5103, 14175) .// 200
(y1, y2, y3, r4, x4, y4, r5, y5, r6) = [305, 138, 224, 85, 225, 252, 106.0, 403, 174.4]
(y1, y2, y3, r4, x4, y4, r5, y5, r6) = res[3]
@printf("r1 = %g; r2 = %g; r3 = %g\n", r1, r2, r3)
@printf("y1 = %g; y2 = %g; y3 = %g; r4 = %g; x4 = %g; y4 = %g; r5 = %g; y5 = %g; y6 = %g\n", y1, y2, y3, r4, x4, y4, r5, y5, r6)
@printf("丁円の直径 = %g\n", 2r4)
plot()
circle2(r1, y1, r1, :red)
circle2(r2, y2, r2, :blue)
circle2(r3, y3, r3, :green)
circle2(x4, y4, r4, :brown)
circle2(r5, y5, r5, :magenta)
circle2(r5, y5, r5, :magenta)
circle2(r6, 0, r6, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r1, y1, " 甲:r1,(r1,y1)", :red, :left, :bottom, delta=delta)
point(r2, y2, " 乙:r2,(r2,y2)", :black, :left, :vcenter)
point(r3, y3, " 丙:r3,(r3,y3)", :green, :center, :top, delta=-delta)
point(x4, y4, " 丁:r4,(r4,y4)", :brown, :center, :bottom, delta=delta)
point(r5, y5, " 戊:r5,(r5,y5)", :magenta, :center, :bottom, delta=delta)
point(r6, 0, "己:r6,(r6,0)", :orange, :center, :bottom, delta=delta)
hline!([0], color=:gray, lw=0.5)
vline!([0], color=:gray, lw=0.5)
else
plot!(showaxis=false)
end
end;