算額(その47)
第2版を公開しました。
https://blog.goo.ne.jp/r-de-r/e/e79c4d483dbfe55f6973570fa73a0c04
岩手県東磐井郡藤沢町 藤勢寺 元治 2 年
http://www.wasan.jp/iwate/fujise.html
1 辺が 1 寸の正方形に2種類の円が収まっている。それぞれの径を求めよ。
下図のように記号を定め解を求める。
using SymPy
function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))
end;
using SymPy
function intersection(x1, y1, x2, y2, x3, y3, x4, y4)
p1, p2, p3, p4 = sympy.Point(x1, y1), sympy.Point(x2, y2), sympy.Point(x3, y3), sympy.Point(x4, y4)
l1 = sympy.Line(p1, p2)
l2 = sympy.Line(p3, p4)
l1.intersection(l2)
end;
@syms r1::positive, r2::positive;
x1, x2 は半径 r1, r2 で表すことができる。
x1 = (1 + sqrt(Sym(2))) * r1
x2 = (1 - sqrt(Sym(2)) * r2)/2
(x1, x2) |> println
(r1*(1 + sqrt(2)), -sqrt(2)*r2/2 + 1/2)
以下の連立方程式を解くが,solve() では解けないので,nlsolve() で解く。
eq1 = distance(r1, x1, 1, 0, x1, r1) - r1;
eq2 = distance(r1, x1, 1, 0, x2, x2) - r2;
# solve([eq1, eq2], (r1, r2))
eq1 |> println
-r1 + sqrt((r1 - r1*(14*r1^2 + 10*sqrt(2)*r1^2 - 10*r1 - 7*sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2 + (r1*(1 + sqrt(2)) - r1*(4*sqrt(2)*r1^2 + 6*r1^2 - 2*r1 - sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2)
eq2 |> println
-r2 + sqrt((-sqrt(2)*r2/2 + 1/2 - (-6*r1^2*r2 - 4*sqrt(2)*r1^2*r2 + 18*r1^2 + 13*sqrt(2)*r1^2 + 5*sqrt(2)*r1*r2 + 8*r1*r2 - 4*sqrt(2)*r1 - 5*r1 - 2*r2 - sqrt(2)*r2 + 1 + sqrt(2))/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))))^2 + (-r1*(-10*sqrt(2)*r1*r2 - 14*r1*r2 + 4*r1 + 3*sqrt(2)*r1 + 4*r2 + 3*sqrt(2)*r2 + 2*sqrt(2) + 3)/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))) - sqrt(2)*r2/2 + 1/2)^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 = u
return [
-r1 + sqrt((r1 - r1*(14*r1^2 + 10*sqrt(2)*r1^2 - 10*r1 - 7*sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2 + (r1*(1 + sqrt(2)) - r1*(4*sqrt(2)*r1^2 + 6*r1^2 - 2*r1 - sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2),
-r2 + sqrt((-sqrt(2)*r2/2 + 1/2 - (-6*r1^2*r2 - 4*sqrt(2)*r1^2*r2 + 18*r1^2 + 13*sqrt(2)*r1^2 + 5*sqrt(2)*r1*r2 + 8*r1*r2 - 4*sqrt(2)*r1 - 5*r1 - 2*r2 - sqrt(2)*r2 + 1 + sqrt(2))/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))))^2 + (-r1*(-10*sqrt(2)*r1*r2 - 14*r1*r2 + 4*r1 + 3*sqrt(2)*r1 + 4*r2 + 3*sqrt(2)*r2 + 2*sqrt(2) + 3)/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))) - sqrt(2)*r2/2 + 1/2)^2)
];
end;
iniv = [0.1, 0.2];
res = nls(h, ini=iniv)
([0.09990042250462969, 0.18946869098150596], 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; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1, r2 = [0.09990042250462969, 0.18946869098150596]
x1, x2 = r1*(1 + sqrt(2)), -sqrt(2)*r2/2 + 1/2
res = intersection(0, 1, x1, r1, r1, x1, 1, 0)
cros = res[1].coordinates[1] # 3 直線の交点
p = plot([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:black, lw=0.5)
circle(r1, x1, r1)
circle(x1, r1, r1)
circle(1-r1, 1-x1, r1)
circle(1-x1, 1-r1, r1)
circle(x2, x2, r2, :red)
circle(sqrt(2)*r2 + x2, sqrt(2)*r2 + x2, r2, :red)
plot!([0, 1], [1, 0], color=:black, lw=0.5)
plot!([cros, 0], [cros, 1], color=:black, lw=0.5)
plot!([cros, 1], [cros, 0], color=:black, lw=0.5)
plot!([cros, 0], [cros, 0], color=:black, lw=0.5)
plot!([1-cros, 0], [1-cros, 1], color=:black, lw=0.5)
plot!([1-cros, 1], [1-cros, 0], color=:black, lw=0.5)
plot!([1-cros, 1], [1-cros, 1], color=:black, lw=0.5)
if more
point(x1, r1, "(x1,r1)", :green, :center)
point(r1, x1, "(r1,x1)", :green, :center)
point(x2, x2, "(x2,x2)", :green, :center)
end
display(p)
savefig("fig1.png")
end
draw(true)