算額（その1335）

算額（その1335）

京都府京都市上京区北野馬喰町 北野天満宮 明治12年(1879)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円5個（円欠）
#Julia, #SymPy, #算額, #和算

外円を水平な弦で切り取り，大円 1 個，小円 3 個を容れる。大円，小円の直径が 4 寸，1 寸のとき，外円の直径はいかほどか。

外円を切り取る水平な弦と y 軸の交点座標を (0, y); y < 0
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x, r1 + y)
小円の半径と中心座標を r2, (x2, y2), (x + a, y + r2), (x - a, y + r2)
大円と小円の中心の水平距離を a = 2sqrt(r1*r2)
とおき，以下の連立方程式を解く。

SymPy では r1, r2 についての一般解（数式解）を求めることができないので，特定の値を指定して解を求める。

include("julia-source.txt");

using SymPy
@syms R::positive, y::negative, a::positive,
     r1::positive, x::negative, r2::positive,
     x2::positive, y2::positive;
(r1, r2) = (4, 1) ./ Sym(2)
a = 2sqrt(r1*r2)
eq1 = x^2 + (r1 + y)^2 - (R - r1)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = (x - a)^2 + (y + r2)^2 - (R - r2)^2
eq4 = (x2 - x)^2 + (r1 + y - y2)^2 - (r1 + r2)^2
eq5 = (x2 - x - a)^2 + (y2 - y - r2)^2 - 4r2^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (R, y, x, x2, y2))[1]

   (-6/19 + 24*sqrt(6)/19, 98/95 - 88*sqrt(6)/95, 44/95 - 24*sqrt(6)/95, 1094/475 - 6*sqrt(6)/475, 1569/950 - 288*sqrt(6)/475)

大円，小円の直径が 4 寸, 1 寸のとき，外円の直径は 12(4√6 - 1)/19 = 5.556605665978554 寸である。

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

   r1 = 2;  r2 = 0.5;  R = 2.7783;  y = -1.23742;  a = 2;  x = -0.155661;  x2 = 2.27222;  y2 = 0.166415

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (4/2, 1/2)
   (R, y, a, x, x2, y2) = (-6/19 + 24*sqrt(6)/19, 98/95 - 88*sqrt(6)/95, 2*sqrt(299263 - 3282*sqrt(6))*(3*sqrt(6) + 547)/299155, 4*(311 - 171*sqrt(6))*sqrt(299263 - 3282*sqrt(6))/1495775, sqrt(1197052/225625 - 13128*sqrt(6)/225625), 1569/950 - 288*sqrt(6)/475)
   @printf("大円，小円の直径が %g, %g のとき，外円の直径は %g である。\n", 2r1, 2r2, 2R)
   @printf("r1 = %g;  r2 = %g;  R = %g;  y = %g;  a = %g;  x = %g;  x2 = %g;  y2 = %g\n", r1, r2, R, y, a, x, x2, y2)
   θ = atand(y, sqrt(R^2 - y^2))
   plot()
   circle(0, 0, R, :green, beginangle=θ, endangle=180-θ, n=1000)
   segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y)
   circle(x, r1 + y, r1)
   circle(x2, y2, r2, :blue)
   circle(x + a, y + r2, r2, :blue)
   circle(x - a, y + r2, r2, :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)
       point(0, R, "R", :green, :center, :bottom, delta=delta/2)
       point(0, 0, "外円:R,(0,0)", :green, :center, delta=-delta/2)
       point(x, r1 + y, "大円:r1,(x,r1+y)", :red, :center, delta=-delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :blue, :center, :bottom, delta=delta/2)
       point(x + a, y + r2, "小円:r2\n(x+a,y+r2)", :blue, :center, :bottom, delta=delta/2)
       point(x - a, y + r2, "小円:r2\n(x-a,y+r2)", :blue, :center, :bottom, delta=delta/2)
   end  
end;

draw(true)

算額（その1334）

算額（その1334）

七二 加須市大字外野 棘脱地蔵堂 明治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)