算額（その447）
河□狭山牛頭天王社 林助右衛門自弘 寛政7年乙卯11月(1795)
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
図のように 12 個の円がある。甲円，乙円，丙円の直径がそれぞれ 43寸7分5厘，51寸0分3厘，141寸7分5厘のとき丁円の直径はいかほどか。
左右対称なので，右半分の図形で方程式を立てる。
甲円の半径と中心座標を 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;
