算額（その143）
岐阜県大垣市西外側町 大垣八幡宮 天保年間
http://www.wasan.jp/gifu/ogakihatiman.html
第8問: 半円の中に，貞，亨，元，利，天，地，人の 7 個の円がある。それぞれの円の半径と中心の x 座標を求めよ。
半円の半径を 1，各円の半径を r1 ~ r7，x 座標を x1 ~ x7 と置き，式を立て，解く。なお 算額では 「地径，人径をもって天径を求めよ］とあるが，下記のように方程式は 13 本なので，未知数は 13 でよいので，人円の半径だけが既知であるとすればよい。
貞: 半径 = r1, 中心の x 座標 = x1
亨: 半径 = r2, 中心の x 座標 = x2
元: 半径 = r3, 中心の x 座標 = x3
利: 半径 = r4, 中心の x 座標 = x4
天: 半径 = r5, 中心の x 座標 = x5
地: 半径 = r6, 中心の x 座標 = x6
人: 半径 = 既知, 中心の x 座標 = x7
solve() ではなく，nlsolve() による。
using SymPy
@syms r1, r2, r3, r4, r5, r6, r7, x1, x2, x3, x4, x5, x6, x7
r7 = 0.04
eq01 = (1 - x1)^2 + r1^2 - (1 - r1)^2
eq02 = (1 - x2)^2 + r2^2 - (1 - r2)^2
eq03 = (x3 - 1)^2 + r3^2 - (1 - r3)^2
eq04 = (x4 - 1)^2 + r4^2 - (1 - r4)^2
eq12 = (x2 - x1)^2 + (r2 - r1)^2 - (r1 + r2)^2
eq17 = (x7 - x1)^2 + (r1 - r7)^2 - (r1 + r7)^2
eq23 = (x3 - x2)^2 + (r3 - r2)^2 - (r2 + r3)^2
eq25 = (x5 - x2)^2 + (r2 - r5)^2 - (r2 + r5)^2
eq27 = (x2 - x7)^2 + (r2 - r7)^2 - (r2 + r7)^2
eq34 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq35 = (x3 - x5)^2 + (r3 - r5)^2 - (r3 + r5)^2
eq36 = (x6 - x3)^2 + (r3 - r6)^2 - (r3 + r6)^2
eq46 = (x4 - x6)^2 + (r4 - r6)^2 - (r4 + r6)^2;
# res = solve([eq01, eq02, eq03, eq04, eq12, eq17, eq23, eq25, eq27, eq34, eq35, eq36, eq46])
println(eq01, ",")
println(eq02, ",")
println(eq03, ",")
println(eq04, ",")
println(eq12, ",")
println(eq17, ",")
println(eq23, ",")
println(eq25, ",")
println(eq27, ",")
println(eq34, ",")
println(eq35, ",")
println(eq36, ",")
println(eq46, ",")
r1^2 - (1 - r1)^2 + (1 - x1)^2,
r2^2 - (1 - r2)^2 + (1 - x2)^2,
r3^2 - (1 - r3)^2 + (x3 - 1)^2,
r4^2 - (1 - r4)^2 + (x4 - 1)^2,
(-r1 + r2)^2 - (r1 + r2)^2 + (-x1 + x2)^2,
(r1 - 0.04)^2 - (r1 + 0.04)^2 + (-x1 + x7)^2,
(-r2 + r3)^2 - (r2 + r3)^2 + (-x2 + x3)^2,
(r2 - r5)^2 - (r2 + r5)^2 + (-x2 + x5)^2,
(r2 - 0.04)^2 - (r2 + 0.04)^2 + (x2 - x7)^2,
(r3 - r4)^2 - (r3 + r4)^2 + (-x3 + x4)^2,
(r3 - r5)^2 - (r3 + r5)^2 + (x3 - x5)^2,
(r3 - r6)^2 - (r3 + r6)^2 + (-x3 + x6)^2,
(r4 - r6)^2 - (r4 + r6)^2 + (x4 - x6)^2,
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
println(params...)
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)
(r1, r2, r3, r4, r5, r6, x1, x2, x3, x4, x5, x6, x7) = u
return [
r1^2 - (1 - r1)^2 + (1 - x1)^2,
r2^2 - (1 - r2)^2 + (1 - x2)^2,
r3^2 - (1 - r3)^2 + (x3 - 1)^2,
r4^2 - (1 - r4)^2 + (x4 - 1)^2,
(-r1 + r2)^2 - (r1 + r2)^2 + (-x1 + x2)^2,
(r1 - 0.04)^2 - (r1 + 0.04)^2 + (-x1 + x7)^2,
(-r2 + r3)^2 - (r2 + r3)^2 + (-x2 + x3)^2,
(r2 - r5)^2 - (r2 + r5)^2 + (-x2 + x5)^2,
(r2 - 0.04)^2 - (r2 + 0.04)^2 + (x2 - x7)^2,
(r3 - r4)^2 - (r3 + r4)^2 + (-x3 + x4)^2,
(r3 - r5)^2 - (r3 + r5)^2 + (x3 - x5)^2,
(r3 - r6)^2 - (r3 + r6)^2 + (-x3 + x6)^2,
(r4 - r6)^2 - (r4 + r6)^2 + (x4 - x6)^2,
]
end;
iniv = [0.08, 0.35, 0.44, 0.17, 0.10, 0.10, 0.01, 0.48, 1.31, 1.81, 0.83, 1.56, 0.21]
res = nls(H, ini=iniv)
println(res)
([0.09120478381765836, 0.3506438413377939, 0.4657568995543326, 0.16582757446843538, 0.10052335452190926, 0.06504524638060978, 0.09579292616973895, 0.45345419466945663, 1.2616986833962582, 1.8175236088720186, 0.8289426508101116, 1.6098094260341749, 0.21659344525316562], true)
算額では，天 = 8 * 地 * 人 / (2 * sqrt( 地 * 人) + 地 + 人) としている。
地円の半径は r6，人は 0.04 と与えたので，天円の半径は，以下のようになる。
(r1, r2, r3, r4, r5, r6, x1, x2, x3, x4, x5, x6, x7) = res[1]
8 * r6 * 0.04 / (2 * sqrt( r6 * 0.04) + r6 + 0.04)
0.1005233545219093
r5 # 求められた天円の半径
0.10052335452190926
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r7 = 0.04
(r1, r2, r3, r4, r5, r6, x1, x2, x3, x4, x5, x6, x7) = res[1]
plot()
circle(1, 0, 1, :black, beginangle=0, endangle=180)
circle(x1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :brown)
circle(x4, r4, r4, :green)
circle(x5, r5, r5, :magenta)
circle(x6, r6, r6, :orange)
circle(x7, r7, r7, :purple)
if more == true
println("貞: 半径 = $r1, 中心の x 座標 = $x1")
println("亨: 半径 = $r2, 中心の x 座標 = $x2")
println("元: 半径 = $r3, 中心の x 座標 = $x3")
println("利: 半径 = $r4, 中心の x 座標 = $x4")
println("天: 半径 = $r5, 中心の x 座標 = $x5")
println("地: 半径 = $r6, 中心の x 座標 = $x6")
println("人: 半径 = $r7, 中心の x 座標 = $x7")
point(x1, r1, "貞")
point(x2, r2, "亨", :blue)
point(x3, r3, "元", :brown)
point(x4, r4, "利", :green)
point(x5, r5, "天", :magenta)
point(x6, r6, "地", :orange)
point(x7, r7, "人", :purple)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
segment(0, 0, 2, 0, :black)
end;
貞: 半径 = 0.09120478381765836, 中心の x 座標 = 0.09579292616973895
亨: 半径 = 0.3506438413377939, 中心の x 座標 = 0.45345419466945663
元: 半径 = 0.4657568995543326, 中心の x 座標 = 1.2616986833962582
利: 半径 = 0.16582757446843538, 中心の x 座標 = 1.8175236088720186
天: 半径 = 0.10052335452190926, 中心の x 座標 = 0.8289426508101116
地: 半径 = 0.06504524638060978, 中心の x 座標 = 1.6098094260341749
人: 半径 = 0.04, 中心の x 座標 = 0.21659344525316562
