算額（その47）

算額（その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)