算額（その300）

算額（その300）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（80)
長野県長野市若槻東條 雁田八幡宮 文化3年(1806) 文面一部欠損
県内の算額（274)
長野県長野市北尾張部 連開庵 奉納年不明

二等辺三角形に全円が内接し，斜線で区切られた三角形にそれぞれ等円が内接している。全円の直径が 3 寸，中鉤（高さ）が 4 寸のとき，等円の直径を求めよ。

5元連立方程式は SymPy の solve() では解けないが，

include("julia-source.txt");

using SymPy

@syms h::positive, a::positive, w::positive, x0::positive,
     r0::positive, r1::positive, x1::positive, x2::positiv;

h = 4
r0 = 3//2
x = a/2
l = sqrt(x^2 + h^2)
m = sqrt((a - x)^2 + h^2)
eq1 = (2w + l)r0 - w*h
eq2 = (a + 2l)r1 + (2w - a + l)r1 - w*h
eq3 = (w - x)^2 + h^2 - w^2
eq4 = distance(x, h, w, 0, x2, r1) - r1^2
eq5 = distance(x, h, w, 0, x0, r0) - r0^2;
@syms d
eq4 = apart(eq4, d) |> numerator
eq5 = apart(eq5, d) |> numerator;

# solve([eq1, eq2, eq3, eq4, eq5], (a, w, x0, x2, r1))

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

   -w + 3*sqrt(a^2/4 + 16)/2,  # eq1
   r1*(a + 2*sqrt(a^2/4 + 16)) + r1*(-a + 2*w + sqrt(a^2/4 + 16)) - 4*w,  # eq2
   -w^2 + (-a/2 + w)^2 + 16,  # eq3
   16*a*r1*w - 16*a*r1*x2 - 64*r1^2 - 32*r1*w^2 + 32*r1*w*x2 + 64*w^2 - 128*w*x2 + 64*x2^2,  # eq4
   24*a*w - 24*a*x0 + 16*w^2 - 80*w*x0 + 64*x0^2 - 144,  # eq5

[eq1, eq3] は他とは独立しているので，これをまず解く。

solve([eq1, eq3], (a, w))

   1-element Vector{Tuple{Sym, Sym}}:
    (2*sqrt(2), 9*sqrt(2)/2)

a, w を既知として方程式 [eq2, eq4, eq5] を組み直す。

@syms h::positive, x0::positive,
     r0::positive, r1::positive, x1::positive, x2::positiv;

h = 4
r0 = 3//2
a = 2sqrt(Sym(2))    # 既知となった
w = 9sqrt(Sym(2))/2  # 既知となった
x = a/2
l = sqrt(x^2 + h^2)
m = sqrt((a - x)^2 + h^2)
eq1 = (2w + l)r0 - w*h
eq2 = (a + 2l)r1 + (2w - a + l)r1 - w*h
eq3 = (w - x)^2 + h^2 - w^2
eq4 = distance(x, h, w, 0, x2, r1) - r1^2
eq5 = distance(x, h, w, 0, x0, r0) - r0^2;
@syms d
eq4 = apart(eq4, d) |> numerator
eq5 = apart(eq5, d) |> numerator;

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

   0,  # eq1
   18*sqrt(2)*r1 - 18*sqrt(2),  # eq2
   0,  # eq3
   -32*r1^2/81 + 56*sqrt(2)*r1*x2/81 - 56*r1/9 + 32*x2^2/81 - 32*sqrt(2)*x2/9 + 16,  # eq4
   32*x0^2/81 - 68*sqrt(2)*x0/27 + 52/9,  # eq5

解は 4 組得られる。

res = solve([eq2, eq4, eq5], (x0, x2, r1))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (3*sqrt(2)/2, 5*sqrt(2)/2, 1)
    (3*sqrt(2)/2, 19*sqrt(2)/4, 1)
    (39*sqrt(2)/8, 5*sqrt(2)/2, 1)
    (39*sqrt(2)/8, 19*sqrt(2)/4, 1)

names = ("x0", "x2", "r1")
for i in 1:length(res)
   println("result $i")
   for (j, name) in enumerate(names)
       @printf("　　%s = %.10f\n", name, res[i][j])
   end
end

   result 1
   　　x0 = 2.1213203436
   　　x2 = 3.5355339059
   　　r1 = 1.0000000000
   result 2
   　　x0 = 2.1213203436
   　　x2 = 6.7175144213
   　　r1 = 1.0000000000
   result 3
   　　x0 = 6.8942911166
   　　x2 = 3.5355339059
   　　r1 = 1.0000000000
   result 4
   　　x0 = 6.8942911166
   　　x2 = 6.7175144213
   　　r1 = 1.0000000000

x0 , x2 < 6 < w ゆえ，1 組目のものが適解で，等円の半径は 1 ゆえ，直径は 2 寸。
a = 2*sqrt(2)
w = 9*sqrt(2)/2
x0 = 3*sqrt(2)/2
x2 = 5*sqrt(2)/2
r1 = 1

数値解は nlsolve() で求めることもできる。

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, w, x0, x2, r1) = u
   return [
       -w + 3*sqrt(a^2/4 + 16)/2,  # eq1
       r1*(a + 2*sqrt(a^2/4 + 16)) + r1*(-a + 2*w + sqrt(a^2/4 + 16)) - 4*w,  # eq2
       -w^2 + (-a/2 + w)^2 + 16,  # eq3
       16*a*r1*w - 16*a*r1*x2 - 64*r1^2 - 32*r1*w^2 + 32*r1*w*x2 + 64*w^2 - 128*w*x2 + 64*x2^2,  # eq4
       24*a*w - 24*a*x0 + 16*w^2 - 80*w*x0 + 64*x0^2 - 144,  # eq5
   ]
end;

r3 = 1
h = 4
r0 = 3//2
iniv = [big"2.8", 6.3, 2.15, 3.5, 0.98]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[2.828427124746190097603377448419396157139390481511502720472215083707481854103854, 6.363961030678927719607599258943641353563445321336646285922337820157100414824585, 2.121320343559642573202533086283528300895624077402693816025820728389215128629599, 3.535533905932737622004221810468008036286461682938140594670174743019746380147886, 1.000000000000000000000000000000000000000569878492122589984801745785913706379155], true)

names = ("a", "w", "x0", "x2", "r1")
for (i, name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

    a = 2.828427
    w = 6.363961
   x0 = 2.121320
   x2 = 3.535534
   r1 = 1.000000

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, r0) = (4, 3//2)
   (a, w, x0, x2, r1) = (2*sqrt(2), 9*sqrt(2)/2, 3*sqrt(2)/2, 5*sqrt(2)/2, 1)
   x = a/2
   plot([0, w, x, 0], [0, 0, h, 0], color=:black, lw=0.5)
   circle(x0, r0, r0)
   circle(x, r1, r1, :blue)
   circle(x2, r1, r1, :blue)
   segment(x, h, a, 0)
   if more
       point(x0, r0, " 全円:r0,(x0,r0)\n", :red, :left, :bottom)
       point(x, r1, " 等円:r1\n (x,r1)", :blue)
       point(x2, r1, " 等円:r1\n (x2,r1)", :blue)
       point(x, h, " (x,h)", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
       point(w, 0, " w", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;