算額（その320）

算額（その320）
解析解は 算額（その834）を参照のこと
長野県長野市戸隠 戸隠山中院権現堂 安永5年(1776)
中村信弥「改訂増補 長野県の算額」(26)
http://www.wasan.jp/zoho/zoho.html

外円の中に中円1個，小円2個が入っている。外円の面積から中円・小円の面積を引くと 120 歩である。また，中円と小円の直径の差は 5 寸である。小円の直径を求めよ。

 

以下の連立方程式を解く。
eq4 の 3 は，この算額で用いる円周率の近似値である(120 はこれに依存している）。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, y::negative;
@syms R, r1, r2, y;

eq1 = r1 - r2 - 5//2
eq2 = r2^2 + (R - r1 - y)^2 - (r1 + r2)^2
eq3 = r2^2 + y^2 - (R - r2)^2
eq4 = 3*(R^2 - (r1^2 + 2r2^2)) - 120
# eq4 = r1 - 2r2
res = solve([eq1, eq2, eq3, eq4], (R, r1, r2, y))

   2-element Vector{NTuple{4, Sym}}:
    (-sqrt(185)/2, 5/2, 0, -sqrt(185)/2)
    (sqrt(185)/2, 5/2, 0, sqrt(185)/2)

この連立方程式は正しいにも関わらず，小円の半径 = 0 という不適切な解を与える。

nlsolve() で数値解を求めると解が求まる。

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

   r1 - r2 - 5/2,  # eq1
   r2^2 - (r1 + r2)^2 + (R - r1 - y)^2,  # eq2
   r2^2 + y^2 - (R - r2)^2,  # eq3
   3*R^2 - 3*r1^2 - 6*r2^2 - 120,  # 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)
   (R, r1, r2, y) = u
   return [
       2*r1 - 2*r2 - 5,  # eq1
       r2^2 - (r1 + r2)^2 + (R - r1 - y)^2,  # eq2
       r2^2 + y^2 - (R - r2)^2,  # eq3
       3*R^2 - 3*r1^2 - 6*r2^2 - 120,  # eq4
   ]
end;

iniv = [big"10.0", 6, 4, -5]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[10.56270582162731798080457930499016663309319075011358623115618193017898504536659, 6.40671194097832904571536744155326822289140318621162475377480880617077430209802, 3.90671194097832904571536744155326822289140318621162475377480880881194616675602, -5.388864105677014011100110408852971086316399261553975564921545695903341835297999], true)

   R = 10.562706;  r1 = 6.406712;  r2 = 3.906712;  y = -5.388864
   小円の直径 = 7.813423881956658

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2, y) = res[1]
   @printf("R = %.6f;  r1 = %.6f;  r2 = %.6f;  y = %.6f\n", R, r1, r2, y)
   println("小円の直径 = $(Float64(2r2))")
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r1, r1)
   circle(r2, y, r2, :blue)
   circle(-r2, y, r2, :blue)
   if more
       point(0, R - r1, " R-r1")
       point(R, 0, " R")
       point(r2, y, "(r2,y)")
       annotate!(0.8, -4, text("小円の直径 = " * @sprintf("%.5f", 2r2), :left, 9))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;