算額（その247）

算額（その247）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（229)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

菱形の中に日円，月円，星円がそれぞれ 1 個，3 個，2 個入っている。
日円，月円の径を知って，星円の径を求めよ。

菱形の中心を原点とする。x 軸上の頂点の x 座標を x，y 軸上の頂点の y 座標を y とする。
日円の半径，中心座標を r1, (0, y1) とする。
真ん中の月円の半径，中心座標を r2, (0, y2) とする。
右の月円の半径，中心座標を r3, (x3, y3) とする。r3 = r2 である。
星円の半径，中心座標を r4, (x4, r4) とする。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::positive, r2::positive, y2::negative,
   r3::positive, x3::positive, y3::positive,
   r4::positive, x4::positive, y4::positive,
   x::positive, y::positive;

r3 = r2
eq1 = x4^2 + (y1 - y4)^2 - (r1 + r4)^2
eq2 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq3 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2
eq4 = x3^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = distance(0, y, x, 0, 0, y1) - r1^2
eq6 = distance(0, y, x, 0, x4, y4) - r4^2
eq7 = distance(0, -y, x, 0, 0, y2) - r2^2
eq8 = distance(0, -y, x, 0, x3, y3) - r3^2
eq9 = (sqrt(x^2 + y^2)/x + 1) * (r1 + r2) - 2y;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
println(eq8, ",  # eq8")
println(eq9, ",  # eq9")

   x4^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq1
   -(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq2
   x3^2 - (r1 + r2)^2 + (y1 - y3)^2,  # eq3
   -4*r2^2 + x3^2 + (y2 - y3)^2,  # eq4
   -r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2,  # eq5
   -r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2,  # eq6
   -r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2,  # eq7
   -r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2,  # eq8
   -2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2),  # eq9

日円，月円の半径を 5, 4 とする。

using Printf
using Plots

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) = (5, 4)
   (r4, y1, y2, x3, y3, x4, y4, y, x) = u
   return [
       x4^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq1
       -(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq2
       x3^2 - (r1 + r2)^2 + (y1 - y3)^2,  # eq3
       -4*r2^2 + x3^2 + (y2 - y3)^2,  # eq4
       -r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2,  # eq5
       -r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2,  # eq6
       -r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2,  # eq7
       -r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2,  # eq8
       -2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2),  # eq9
   ]
end;

iniv = [big"1.0", 1.9, -3.1, 3.1, -1.3, 4, 1.9, 5.4, 9.3]
res = nls(H, ini=iniv);
println(res[2] ? "収束した" : "収束していない")

   収束した

names = ["r4", "y1", "y2", "x3", "y3", "x4", "y4", "y ", "x "]
for i in 1:length(names)
   println("$(names[i]) = $(Float64(res[1][i]))")
end

   r4 = 1.7470608602667996
   y1 = 3.9418436943485617
   y2 = -5.058156305651438
   x3 = 7.166451331820933
   y3 = -1.5026007500958825
   x4 = 6.740960675000237
   y4 = 4.228687607033678
   y  = 9.523406750862943
   x  = 19.195039966835868

星円の半径は 1.7470608602667996。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, y1, y2, x3, y3, x4, y4, y, x) = res[1]
   (r1, r2) = (5, 4)
   r3 = r2
   # y = 5/2 + 5*sqrt(3)/3
   # x = sqrt(3) * y
   @printf("r1 = %.5f;  r2 = %.5f;  \nr4 = %.5f;  y1 = %.5f;  y2 = %.5f;  \nx3 = %.5f;  y3 = %.5f;  x4 = %.5f;  \ny4 = %.5f;  y = %.5f;  x = %.5f\n", r1, r2, r4, y1, y2, x3, y3, x4, y4, y, x)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, y1, r1)
   circle(0, y2, r2, :blue)
   circle(x3, y3, r3, :blue)
   circle(-x3, y3, r3, :blue)
   circle(x4, y4, r4, :green)
   circle(-x4, y4, r4, :green)
   if more == true
       point(0, y1, " y1\n 日円", :red, :left, :bottom)
       point(0, y2, " y2\n 月円", :blue, :left, :bottom)
       point(x3, y3, "(x3,y3)\n月円", :blue, :left, :center)
       point(x4, y4, "(x4,y4)\n星円", :green, :left, :center)
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, "x", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;