算額（その241）

算額（その241）

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

2 個の大円が交わり，その中に長方形と小円が 4 個入っている。小円の径を知って，大円の径を求めよ。

大円に接する長方形の右上の頂点座標を (x, y) とおく。
大円，小円の半径と中心座標を以下のようにおく。

大円: r1, (x1, 0)　ただし与えられた条件より r1 = 1 とする。
小円: r2, (0, y + r2) および (x + r2, 0)

以下の連立方程式を立てるが，問では条件が一つ足りない。そこで，x と y の間に x = ny のような制約条件を追加して解く（以下では x = 2y とした）。

なお，nlsolve() ではこの制約条件を追加しなくても一応の解は得られる。

include("julia-source.txt");

using SymPy

@syms x::positive, y::positive, r1::positive, x1::positive, r2::positive;

r2 = 1
x = 2y
eq1 = (x - x1)^2 + y^2 - r1^2 |> expand
eq2 = x - x1 + 2r2 - r1 |> expand
eq3 = x1^2 + (y + r2)^2 - (r1 - r2)^2 |> expand
eq4 = ((2r1 - 2r2)^2 + y^2) + (4r2^2 + y^2) - (2r1)^2 |> expand;

res = solve([eq1, eq2, eq3], (x1, r1, y))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (-5*sqrt(2)*(8*sqrt(2) + 71)^(2/3)/289 + sqrt(2)*(8*sqrt(2) + 71)^(1/3)/34 + 21/16 + 99*(8*sqrt(2) + 71)^(1/3)/272 + 421*(8*sqrt(2) + 71)^(2/3)/4624, -3*sqrt(2)*(8*sqrt(2) + 71)^(2/3)/289 - sqrt(2)*(8*sqrt(2) + 71)^(1/3)/34 + 173*(8*sqrt(2) + 71)^(1/3)/272 + 715*(8*sqrt(2) + 71)^(2/3)/4624 + 59/16, 3/2 + 17/(4*(sqrt(2) + 71/8)^(1/3)) + (sqrt(2) + 71/8)^(1/3))

    x1 = 4.33657; r1 = 8.92148; y = 5.62902; x = 11.25805

小円の半径を 1 としたので，8.92148 は大円の半径である。

using Printf
function draw(x1, r1, x, y, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1
   @printf("x1 = %.5f; r1 = %.5f; y = %.5f; x = %.5f\n", x1, r1, y, x)
   plot([x, x, -x, -x, x], [-y, y, y, -y, -y], color=:black, lw=0.5)
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle(0, y + r2, r2, :green)
   circle(0, -y - r2, r2, :green)
   circle(x + r2, 0, r2, :green)
   circle(-x - r2, 0, r2, :green)
   if more == true
       point(x1, 0, "x1 ", :red, :right)
       point(x + r2, 0, "x+r2 ", :green, :right)
       point(0, y + r2, "    y+r2", :green, :left, :bottom)
       point(x + 2r2, 0, " x+2r2", :red)
       point(x, y, " (x,y)", :red, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

(x1, r1, y) = res[1]

---

条件不足のまま 4元連立方程式を nlsolve() で x1, r1, x, y について解く。

using SymPy

@syms x::positive, y::positive, r1::positive, x1::positive, r2::positive;

r2 = 1
eq1 = (x - x1)^2 + y^2 - r1^2 |> expand
eq2 = x - x1 + 2r2 - r1 |> expand
eq3 = x1^2 + (y + r2)^2 - (r1 - r2)^2 |> expand
eq4 = ((2r1 - 2r2)^2 + y^2) + (4r2^2 + y^2) - (2r1)^2 |> expand;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")

   -r1^2 + x^2 - 2*x*x1 + x1^2 + y^2,  # eq1
   -r1 + x - x1 + 2,  # eq2
   -r1^2 + 2*r1 + x1^2 + y^2 + 2*y,  # eq3
   -8*r1 + 2*y^2 + 8,  # eq4

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 = 1
   (x1, r1, x, y) = u
   return [
       -r1^2 + x^2 - 2*x*x1 + x1^2 + y^2,  # eq1
       -r1 + x - x1 + 2,  # eq2
       -r1^2 + 2*r1 + x1^2 + y^2 + 2*y,  # eq3
       -8*r1 + 2*y^2 + 8,  # eq4
      ]
end;

iniv = [big"4.0", 9, 10, 6]
res2 = nls(H, ini=iniv);
println((res2));

   (BigFloat[3.183012652585522591832434702426538376856770104448615704443237602882606475256619, 8.076648273506633325075325119429150245928992397056981109398659352229911092241916, 9.259660926092155916907759821855688622785762501505596813841896955962457886997144, 5.320394073189178089410508073953356392935683449161830208448592642794610637849292], true)

   x1 = 3.18301; r1 = 8.07665; y = 5.32039; x = 9.25966