算額(その93)
千葉県君津市 神野寺 明治20(1887)年4月
http://www.wasan.jp/chiba/jinyaji2.html
甲円の下に乙円,丙円,丁円がある。各円の径を求めよ。
甲円の半径を 1 として,以下のように記号を定め,方程式を解く。
なお,solve() では解が求まらないので,nlsolve() を使用し,数値的に解を求める。
using SymPy
@syms r1::positive, r2::positive, r3::positive, x1::positive, x2::positive, x3::positive, x::positive, y::positive;
r = 1
x1 = r1
eq1 = x*x3 + y*(r3 - y) |> expand
eq2 = r3*x - y*(x - x3) |> expand
eq3 = x1^2 + (y - r1)^2 - (r + r1)^2 |> expand
eq4 = x2^2 + (y - r2)^2 - (r + r2)^2 |> expand
eq5 = x3^2 + (y - r3)^2 - (r + r3)^2 |> expand
eq6 = (x2 - x1)^2 + (r2 - r1)^2 - (r2 + r1)^2 |> expand
eq7 = (x3 - x2)^2 + (r3 - r2)^2 - (r3 + r2)^2 |> expand;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x2, x3, x, y))
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
r3*y + x*x3 - y^2
r3*x - x*y + x3*y
r1^2 - 2*r1*y - 2*r1 + y^2 - 1
-2*r2*y - 2*r2 + x2^2 + y^2 - 1
-2*r3*y - 2*r3 + x3^2 + y^2 - 1
r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2
-4*r2*r3 + x2^2 - 2*x2*x3 + x3^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
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, x2, x3, x, y) = u
return [
r3*y + x*x3 - y^2,
r3*x - x*y + x3*y,
r1^2 - 2*r1*y - 2*r1 + y^2 - 1,
-2*r2*y - 2*r2 + x2^2 + y^2 - 1,
-2*r3*y - 2*r3 + x3^2 + y^2 - 1,
r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2,
-4*r2*r3 + x2^2 - 2*x2*x3 + x3^2
];
end;
iniv = [0.1, 0.2, 0.3, 0.4, 0.7, 1.0, 1.1];
res = nls(H, ini=iniv)
([0.13839876084731936, 0.17954306798152128, 0.3287706665765487, 0.45366714256276236, 0.9395827489780465, 1.2683534155545952, 1.268353415554595], true)
using Plots
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; fontsize=10, mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, fontsize, vertical, position, color))
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")
(r1, r2, r3, x2, x3, x, y) = [0.13839876084731936, 0.17954306798152128, 0.3287706665765487,
0.45366714256276236, 0.9395827489780465, 1.2683534155545952, 1.268353415554595]
r = 1
plot()
circle(0, y, 1)
circle(r1, r1, r1, :green)
circle(-r1, r1, r1, :green)
circle(x2, r2, r2, :blue)
circle(-x2, r2, r2, :blue)
circle(x3, r3, r3, :magenta)
circle(-x3, r3, r3, :magenta)
hline!([0], color=:black, lw=0.5)
if more
point(0, y, "y ", :red, :right)
point(0, y-r, "y-r ", :red, :right)
point(r1, r1, "(r1,r1)", :green, :top)
point(x2, r2, "(x2,r2)", :green, :top)
point(x3, r3, "(x3,r3)", :magenta, :top)
point(x, 0, " x", :magenta, :bottom)
segment(0, y, x, 0, :black)
vline!([0], color=:black, lw=0.5)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます