算額（その2134）

算額（その2134）

七二 加須市大字外野 棘脱地蔵堂 明治7年(1874)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：正方形5個
#Julia, #SymPy, #算額, #和算

正方形の中に，東西南北4個の正方形を容れる。南正方形（一番大きい正方形）の一辺の長さが 350 寸のとき，北正方形（一番小さい正方形）の一辺の長さはいかほどか。

図に示すような a 〜 f の長さの相互関係に関する連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     d::positive, e::positive, f::positive, x::positive;
e = a + c - x
eq1 = (a - x)/(x - d) - c/a
eq2 = d/(b - x) - c/a
eq3 = e/(a + b + x - c - a) - c/a
eq4 = e^2 + (a + b - e)^2 - f^2
eq5 = sqrt(c^2 + a^2) + sqrt((a - x)^2 + (x - d)^2) - f;

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   #(a, x, b, c, d, e) = u
   (a, x, b, c, d) = u
   return [
       (a - x)/(-d + x) - c/a,  # eq1
       d/(b - x) - c/a,  # eq2
       (a + c - x)/(b - c + x) - c/a,  # eq3
       -f^2 + (a + c - x)^2 + (b - c + x)^2,  # eq4
       -f + sqrt(a^2 + c^2) + sqrt((a - x)^2 + (-d + x)^2),  # eq5
   ]
end;

f = 350
iniv = BigFloat[454, 350, 550, 171, 75]
iniv = BigFloat[205, 157, 247, 77, 34]
res = nls(H, ini=iniv)
res |> println
res[1][1] + res[1][2] + res[1][3] |> println

   ([204.8607611954352, 157.11193591084313, 249.43871128176957, 81.15124262513827, 36.57329253880388], true)
   611.4114083880479

南正方形の一辺の長さが 350 寸のとき，北正方形の一辺の長さは 157.112 寸である。

全てのパラメータは以下のとおりである。

   f = 350;  a = 204.861;  x = 157.112;  b = 249.439;  c = 81.1512;  d = 36.5733;  e = 128.9

「答」は 165 寸となっている。数値解が安定しないのは，重要な条件を見逃しているのだと思う。

function draw(f, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, x, b, c, d) = res[1]
   e = a + c - x
   axb = a + x + b
   @printf("南正方形の一辺の長さが %g のとき，北正方形の一辺の長さは %g である。\n", f, x)
   @printf("f = %g;  a = %g;  x = %g;  b = %g;  c = %g;  d = %g;  e = %g\n", f, a, x, b, c, d, e)
   plot([0, a + x - d, e + a + x - d, e, 0], [a + c, x, axb - c - a +x, axb, a + c], color=:green, lw=0.5)
   rect(0, 0, axb, axb, :red)
   rect(0, 0, a, a, :red)
   rect(a + x, 0, a + x + b, b, :red)
   rect(a, 0, a + x, x, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       dimension_line(0, 0, a, 0, "a", :black, :center, delta=-delta)
       dimension_line(a, 0, a + x, 0, "x", :black, :center, delta=-delta)
       dimension_line(a + x, 0, a + x + b, 0, "b", :black, :center, delta=-delta)
       dimension_line(a + x, x, a + x -d, x, "d", :black, :center, delta=-delta)
       dimension_line(0, a, 0, a + c, "c ", :black, :right, delta=-delta)
       dimension_line(0, axb, e, axb, "e", :black, :center, :bottom, delta=delta)
       dimension_line(0, a + c, e, axb, " f", :black, :left)
       plot!(xlims=(-5delta, axb + 5delta), ylims=(-5delta, axb + 5delta))
   end  
end;

draw(350, true)