算額（その114）

算額（その114）

福島県田村郡三春町平沢担橋 諏訪神社 大正15年(1926)8月
http://www.wasan.jp/fukusima/miharusuwa.html

三角形の中に菱形，甲円，乙円が入っている。大斜（一番長い辺）が 9 寸，甲円の径が 2 寸，乙円の径が 1 寸のとき，菱形の辺の長さはいかほどか。

以下に，求めるものは異なるが同じ問題がある。

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（222)
長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)

方程式を書きやすくするために左右反転させて，大斜を 18，甲円，乙円の半径を 2, 1，菱形の辺の長さを a として下図のように記号を定める。

算額の問のように，「菱形の辺の長さ」だけを求めるならば，三角形の相似関係を使えば簡単に答えを得ることができる。⊿EFC，⊿ODE，⊿OBC は相似（相似比 1:2:3）

ここでは，解を確認するための図を描くのに必要な座標をすべて求めることにする。

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))^2  # 二乗距離を返す！！
end;

@syms a, x, y, x1, x2, y2, x3, y3;
eq0 = (18 - a) + (18 - a)/2 - 18
eq1 = distance(0, 0, x, y, x1, 2) - 2^2;
eq2 = distance(0, 0, x, y, x2, y2) - 1^2;
eq3 = distance(x3, y3, 18 - a, 0, x1, 2) - 2^2;
eq4 = distance(x, y, 18, 0, x2, y2) - 1^2;
eq5 = y3^2 + (18 - a - x3)^2 - 6^2;
eq6 = 2x - 3x3;
eq7 = 2y - 3y3;

例により，solve() では解けないので，nlsolve() で数値解を求めることにする。

# res = solve([eq0, eq1, eq2, eq3, eq4, eq5, eq6, eq7]);

eq0 |> println
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println

   9 - 3*a/2
   (-x*(x*x1 + 2*y)/(x^2 + y^2) + x1)^2 + (-y*(x*x1 + 2*y)/(x^2 + y^2) + 2)^2 - 4
   (-x*(x*x2 + y*y2)/(x^2 + y^2) + x2)^2 + (-y*(x*x2 + y*y2)/(x^2 + y^2) + y2)^2 - 1
   (x1 - (x1*x3^2 - 24*x1*x3 + 144*x1 + 2*x3*y3 + 12*y3^2 - 24*y3)/(x3^2 - 24*x3 + y3^2 + 144))^2 + (-y3*(x1*x3 - 12*x1 - 12*x3 + 2*y3 + 144)/(x3^2 - 24*x3 + y3^2 + 144) + 2)^2 - 4
   (x2 - (x2*(x^2 - 36*x + y^2 + 324) + y*(x*y2 - x2*y + 18*y - 18*y2))/(x^2 - 36*x + y^2 + 324))^2 + (-y*(x*x2 - 18*x - 18*x2 + y*y2 + 324)/(x^2 - 36*x + y^2 + 324) + y2)^2 - 1
   y3^2 + (12 - x3)^2 - 36
   2*x - 3*x3
   2*y - 3*y3

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)
   (a, x, y, x1, x2, y2, x3, y3) = u
   return [
       9 - 3*a/2,
       (-x*(x*x1 + 2*y)/(x^2 + y^2) + x1)^2 + (-y*(x*x1 + 2*y)/(x^2 + y^2) + 2)^2 - 4,
       (-x*(x*x2 + y*y2)/(x^2 + y^2) + x2)^2 + (-y*(x*x2 + y*y2)/(x^2 + y^2) + y2)^2 - 1,
       (x1 - (x1*x3^2 - 24*x1*x3 + 144*x1 + 2*x3*y3 + 12*y3^2 - 24*y3)/(x3^2 - 24*x3 + y3^2 + 144))^2 + (-y3*(x1*x3 - 12*x1 - 12*x3 + 2*y3 + 144)/(x3^2 - 24*x3 + y3^2 + 144) + 2)^2 - 4,
       (x2 - (x^2*x2 - 36*x*x2 + x*y*y2 + 324*x2 + 18*y^2 - 18*y*y2)/(x^2 - 36*x + y^2 + 324))^2 + (-y*(x*x2 - 18*x - 18*x2 + y*y2 + 324)/(x^2 - 36*x + y^2 + 324) + y2)^2 - 1,
       y3^2 + (12 - x3)^2 - 36,
       2*x - 3*x3,
       2*y - 3*y3,
   ]
end;

iniv = [5.0, 12.0, 7, 8, 12, 5, 8, 4]
res = nls(H, ini=iniv)
println(res)

   ([6.0, 12.125711439226585, 6.8185580517266695, 7.63711610345335, 11.902365677877736, 5.5457053678177815, 8.083807626151057, 4.54570536781778], true)

a = 6 となるが，元の単位では 3 寸である。

using Plots
using Printf

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")
   (a, x, y, x1, x2, y2, x3, y3) = res[1]
   @printf("a  = %8.5f,  x = %8.5f,  y = %8.5f, x1 = %8.5f\n", a, x, y, x1)
   @printf("x2 = %8.5f, y2 = %8.5f, x3 = %8.5f, y3 = %8.5f\n", x2, y2, x3, y3)
   plot([0, 18, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, 2, 2)
   circle(x2, y2, 1)
   plot!([x3+a, x3, 18-a], [y3, y3, 0], color=:black, lw=0.5)
   if more == true
       point(x1, 2, "(x1,2),2", :green, :top)
       point(x3, y3, "E:(x3,y3) ", :green, :right, :bottom)
       point(x3 + a, y3, "  F:(x3+a, y3)")
       point(x, y, "  C:(x,y),1", :green, :center, :bottom)
       point(x2, y2, "(x2,y2),1", :green, :top)
       point(18 - a, 0, " D:18-a")
       point(18, 0, "  B:18")
       point(0, 0, "O ", :green, :right)
       hline!([0], color=:black, lw=0.5, xlims=(-1, 20))
       vline!([0], color=:black, lw=0.5, ylims=(-1, 7))
   else
       plot!(showaxis=false)
   end
end;

   a  =  6.00000,  x = 12.12571,  y =  6.81856, x1 =  7.63712
   x2 = 11.90237, y2 =  5.54571, x3 =  8.08381, y3 =  4.54571