算額（その38） - 算額あれこれ

新潟県糸魚川市天津神社寛政12年(1800)

http://www.wasan.jp/niigata/tensin.html
キーワード：円8個，外円
#Julia #SymPy #算額 #和算 #数学

図のように，外円の中に 5 種類の径（木，火，土，金，水）を持つ円が収まっている。土円，金円の径はそれぞれ 9分9釐（厘），一寸9分9釐（厘）である。水円の径を求めよ。

図のように記号を定め，方程式を解く。水円の径は \(r_3\) である。（プログラムでは，すべて直径ではなく半径を使っている。）

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms　r1::positive, r2::positive, y2::positive, r0::positive, y3::positive, r3::positive, r4::positive, x4::positive, y4::positive, r5::positive, x5::positive, y5::positive;

r1 = 99
r2 = 198
y3 = r3
x5 = r5
x4 = r4
eq1 = x4^2+(2r0-r1-y4)^2-(r1+r4)^2
eq2 = x4^2+(y4-y2)^2-(r2+r4)^2
eq3 = x5^2+(y2-y5)^2-(r2+r5)^2
eq4 = x5^2+(y5-y3)^2-(r3+r5)^2
eq5 = x4^2+(y4-r0)^2-(r0-r4)^2b
eq6 = x5^2+(r0-y5)^2-(r0-r5)^2
eq7 = (x5 - x4)^2 + (y4 - y5)^2 - (r4 + r5)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (y2, r0, r3, r4, y4, r5, y5))

solve() では解が求まらないので，nlsolve() で求める。

function h(u)
y2, r0, r3, r4, y4, r5, y5 = u
return [
r4^2 - (r4 + 99)^2 + (2*r0 - y4 - 99)^2
r4^2 - (r4 + 198)^2 + (-y2 + y4)^2
r5^2 - (r5 + 198)^2 + (y2 - y5)^2
r5^2 + (-r3 + y5)^2 - (r3 + r5)^2
r4^2 + (-r0 + y4)^2 - (r0 - r4)^2
r5^2 - (r0 - r5)^2 + (r0 - y5)^2
(-r4 + r5)^2 - (r4 + r5)^2 + (y4 - y5)^2
];
end;

iniv = [2260., 1521, 580, 326, 2670, 752, 1679];

res = nls(h, ini=iniv)

y2, r0, r3, r4, y4, r5, y5 = res[1]

7-element Vector{Float64}:
2260.2616848711004
1521.243310737389
580.0010348927411
326.36779179808775
2670.68298737871
752.3774162147455
1679.6194053038619

求める水円の径は 580.0010348927411 である。

描画関数プログラムのソースを見る

function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 99
r2 = 198
y2, r0, r3, r4, y4, r5, y5 = [2260.2616848711004, 1521.243310737389, 580.0010348927411, 326.36779179808775, 2670.68298737871, 752.3774162147455, 1679.6194053038619]
y3 = r3
x4 = r4
x5 = r5
@printf("y2 = %g; r0 = %g; r3 = %g; r4 = %g; y4 = %g; r5 = %g; y5 = %g\n", y2, r0, r3, r4, y4, r5, y5)
plot()
circle(0, r0, r0, :green)
circle(0, 2r0-r1, r1)
circle(0, y2, r2, :blue)
circle(0, y3, r3, :brown)
circle(x4, y4, r4, :magenta)
circle(-x4, y4, r4, :magenta)
circle(x5, y5, r5, :black)
circle(-x5, y5, r5, :black)
if more
point(0, r0, "r0 ", :green, :right)
point(0, 2r0-r1, "r1 ", :red, :right)
point(0, y2, "y2 ", :blue, :right)
point(0, y3, "r3 ", :brown, :right)
point(x4, y4, "(x4,y4)", :magenta, :center)
point(x5, y5, "(x5,y5)", :black, :center)
vline!([0], color=:black, linewidth=0.25)
hline!([0], color=:black, linewidth=0.25)
end
end;