算額(その104)
愛媛県四国中央市川之江町 大岡八幡神社 文化元年(1804)4月
http://www.wasan.jp/ehime/kawanoe2.html
外円の中に甲円,乙円が1個ずつ,丙円,丁円,戊円が2つずつ互いに接して入っている。
甲円,乙円の直径をそれぞれ66寸,44寸としたとき,戊円の直径はいくつか。
甲円,乙円の半径をそれぞれ 33,22 とする。丙円,丁円,戊円の中心座標と半径を (x1, y1, r1), (x2, y2, r2), (x3, y3, r3) として,以下の9本の方程式を解く。
条件式と未知数の数が等しいので,解けるはずだが solve() では解けない。
using SymPy
@syms x1::positive, y1::positive, r1::positive, x2::positive, y2::positive, r2::positive, x3::positive, y3::positive, r3::positive;
@syms x1, y1, r1, x2, y2, r2, x3, y3, r3
eq1 = x3^2 + (y3 - 33)^2 - (33 + r3)^2 |> expand
eq2 = x1^2 + (y1 - 33)^2 - (33 + r1)^2 |> expand
eq3 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2 |> expand
eq4 = (x1 - x2)^2 + (y1 - y2)^2 - (r1 + r2)^2 |> expand
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq6 = x2^2 + (y2 - 88)^2 - (22 + r2)^2 |> expand
eq7 = x3^2 + (88 - y3)^2 - (22 + r3)^2 |> expand
eq8 = x2^2 + (y2 - 55)^2 - (55 - r2)^2 |> expand
eq9 = x1^2 + (y1 - 55)^2 - (55 - r1)^2 |> expand
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9]);
今までの経験で,nlsolve() で解くことにする。
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println
-r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3
-r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1
-r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2
-r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2
-r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2
-r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260
-r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260
-r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2
-r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1
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)
(x1, y1, r1, x2, y2, r2, x3, y3, r3) = u
return [
-r3^2 - 66*r3 + x3^2 + y3^2 - 66*y3,
-r1^2 - 66*r1 + x1^2 + y1^2 - 66*y1,
-r1^2 - 2*r1*r3 - r3^2 + x1^2 - 2*x1*x3 + x3^2 + y1^2 - 2*y1*y3 + y3^2,
-r1^2 - 2*r1*r2 - r2^2 + x1^2 - 2*x1*x2 + x2^2 + y1^2 - 2*y1*y2 + y2^2,
-r2^2 - 2*r2*r3 - r3^2 + x2^2 - 2*x2*x3 + x3^2 + y2^2 - 2*y2*y3 + y3^2,
-r2^2 - 44*r2 + x2^2 + y2^2 - 176*y2 + 7260,
-r3^2 - 44*r3 + x3^2 + y3^2 - 176*y3 + 7260,
-r2^2 + 110*r2 + x2^2 + y2^2 - 110*y2,
-r1^2 + 110*r1 + x1^2 + y1^2 - 110*y1
];
end;
iniv = [40.0, 57, 14,
33, 82, 12,
20, 68, 8];
res = nls(H, ini=iniv)
([40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996], false)
収束していないが,ftol=1e-14 が厳しすぎるのである。図に描いてみると十分な精度である。
戊円の直径は 7.5*2 = 15 である。元の単位で言えば 15 寸である。
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=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, y1, r1, x2, y2, r2, x3, y3, r3) = [40.58178048548881, 57.391304347826086, 14.347826086956518, 33.33503397022296, 82.5, 11.785714285714292, 21.21320343559642, 67.5, 7.499999999999996]
plot()
circle(0, 55, 55)
circle(0, 33, 33)
circle(0, 88, 22)
circle(x1, y1, r1, :green)
circle(-x1, y1, r1, :green)
circle(x2, y2, r2, :blue)
circle(-x2, y2, r2, :blue)
circle(x3, y3, r3, :magenta)
circle(-x3, y3, r3, :magenta)
if more
a = @sprintf("丙円 (%.1f, %.1f, %.1f)\n", x1, y1, r1)
b = @sprintf("丁円 (%.1f, %.1f, %.1f)\n", x2, y2, r2)
c = @sprintf("戊円 (%.1f, %.1f, %.1f)", x3, y3, r3)
point(-20, 22, a * b * c, :black, mark=false)
point(0, 33, " 甲 (0,33,33)", :red)
point(0, 88, " 乙 (0,88,22)", :red)
point(x1, y1, " 丙\n (x1,y1;r1)", :green, :top)
point(x2, y2, " 丁\n (x2,y2;r2)", :blue, :top)
point(x3, y3, " 戊\n (x3,y3;r3)", :magenta, :top)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます