算額（その312）
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県菰野町 広幡神社 嘉永5年(1852)
問題
半円の中に円弧があり，半円と円弧に挟まれた領域にいくつかの円が挟まれている。半円の直径，丙円，丁円の直径が与えられたとき，戊円の直径はいくつか。
円弧の半径と中心座標を r0, (0, y0)
半円の半径と中心座標を r, (0, 0)
丙円の半径と中心座標を r1, (x1, y1)
丁円の半径と中心座標を r2, (x2, y2)
戊円の半径と中心座標を r3, (x3, y3)
として，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r0::positive, y0::positive, r::positive,
r1::positive, x1::positive, y1::negative,
r2::positive, x2::positive, y2::negative,
r3::positive, x3::positive, y3::negative;
(r, r1, r2) = (50, 10, 5)
eq1 = x1 ^2 + (y0 - y1)^2 - (r0 + r1)^2
eq2 = x2 ^2 + (y0 - y2)^2 - (r0 + r2)^2
eq3 = x3 ^2 + (y0 - y3)^2 - (r0 + r3)^2
eq4 = x1^2 + y1^2 - (r - r1)^2
eq5 = x2^2 + y2^2 - (r - r2)^2
eq6 = x3^2 + y3^2 - (r - r3)^2
eq7 = (x2 - x1)^2 + (y2 - y1)^2 - (r2 + r1)^2
eq8 = (x3 - x2)^2 + (y3 - y2)^2 - (r3 + r2)^2
eq9 = r^2 + y0^2 - r0^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])
連立方程式が複雑すぎて solve() では解けないので，nlsolve() で数値解を求める。
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
println(eq9, ", # eq9")
x1^2 - (r0 + 10)^2 + (y0 - y1)^2, # eq1
x2^2 - (r0 + 5)^2 + (y0 - y2)^2, # eq2
x3^2 - (r0 + r3)^2 + (y0 - y3)^2, # eq3
x1^2 + y1^2 - 1600, # eq4
x2^2 + y2^2 - 2025, # eq5
x3^2 + y3^2 - (50 - r3)^2, # eq6
(-x1 + x2)^2 + (-y1 + y2)^2 - 225, # eq7
-(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2, # eq8
-r0^2 + y0^2 + 2500, # eq9
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(r0, y0, x1, y1, x2, y2, r3, x3, y3) = u
return [
x1^2 - (r0 + 10)^2 + (y0 - y1)^2, # eq1
x2^2 - (r0 + 5)^2 + (y0 - y2)^2, # eq2
x3^2 - (r0 + r3)^2 + (y0 - y3)^2, # eq3
x1^2 + y1^2 - 1600, # eq4
x2^2 + y2^2 - 2025, # eq5
x3^2 + y3^2 - (50 - r3)^2, # eq6
(-x1 + x2)^2 + (-y1 + y2)^2 - 225, # eq7
-(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2, # eq8
-r0^2 + y0^2 + 2500, # eq9
]
end;
iniv = [big"66.0", 45.0, 31.2, -24.5, 42, -13, 3, 47, -7]
res = nls(H, ini=iniv);
println(res);
(BigFloat[76.12612612612612612612612611757250855206487617275733389169965733209449470577765, 57.40372007954591399580886043401747772776331606928505070760935723414311831179716, 33.42516087186933536638372029599284981641171160639776150802648249557181481914727, -21.97176872010205673632683968491622029305014289094789997653090282550024009759404, 43.6384044716071878394454126067690155978138981377281139789361700897187230488509, -10.9858843600510283681634198397671466786061193950162591863412958011058978766673, 2.226463104325699745547073791346257958834702049836892015992527235142629533435189, 47.52241383500506650373257942562149702505122375947591050606909041806671507003448, -4.8919332392084731919302760073024864943304938091114520230459338837681434939186], true)
円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
半円の直径 = 100.000000, 中心座標 = (0.000000, 0.000000)
円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
丙円の直径 = 20.000000, 中心座標 = (33.425161, -21.971769)
丁円の直径 = 10.000000, 中心座標 = (43.638404, -10.985884)
戊円の直径 = 4.452926, 中心座標 = (47.522414, -4.891933)
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r, r1, r2) = (50, 10, 5)
(r0, y0, x1, y1, x2, y2, r3, x3, y3) = res[1]
@printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
@printf("半円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r, 0, 0)
@printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
@printf("丙円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r1, x1, y1)
@printf("丁円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r2, x2, y2)
@printf("戊円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r3, x3, y3)
plot()
circle(0, y0, r0, :black, beginangle=180, endangle=360)
circle(0, 0, r, beginangle=180, endangle=360)
circle(x1, y1, r1, :blue)
circle(x2, y2, r2, :orange)
circle(x3, y3, r3, :magenta)
if more
point(0, y0, " y0", :black)
point(r, 0, "r ", :red, :right, :bottom)
point(x1, y1, " 丙円:r1,(x1,y1)", :blue, :left, :bottom)
point(x2, y2, " 丁円:r2,(x2,y2)", :orange, :left, :bottom)
point(x3, y3, " 戊円:r3,(x3,y3)", :magenta, :left, :bottom)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
