算額（その567）

算額（その567）

群馬の算額 29-3 山名八幡宮 文化11年
http://takasakiwasan.web.fc2.com/gunnsann/g029-3.html

楕円の中に長径に沿って等円を描く。両端の円は楕円の頂点で楕円に接している。楕円の長径，短径（差渡し径）が 15 寸，5 寸のとき，等円の個数の最小値はいくつか。

長径，短径が a, b の楕円の曲率円（長径の端の 1 点で内接する半径が最大の円）の半径は，b^2/a である。

等円の個数の最小値は「楕円の長径/曲率円の半径」である。
a = 15/2, b = 5/12 のとき 9 個である。

a = 15//2
b = 5//2
println("曲率円の半径 = $(b^2/a)")
a/(b^2/a)

   曲率円の半径 = 5//6

   9//1

   a = 7.5;  b = 2.5;  曲率円の半径 = 0.833333

include("julia-source.txt");

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (15, 5) .// 2
   r = b^2/a
   @printf("a = %g;  b = %g;  曲率円の半径 = %g\n", a, b, b^2/a)
   plot()
   ellipse(0, 0, a, b, color=:red)
   for x in -a + r:2r:a - r
       circle(x, 0, r, :blue)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, b, " b", :red, :left, :bottom, delta=delta)
       point(a, 0, " a", :red, :left, :bottom, delta=delta)
       point(r - a, 0, "r-a", :blue, :center, :bottom, delta=delta)
       point(a - r, 0, "a-r", :blue, :center, :bottom, delta=delta)
   end
end;

算額（その566）

算額（その566）

群馬の算額 56-4改 八幡宮 天保5年
http://takasakiwasan.web.fc2.com/gunnsann/g056-4kai.html

正方形の中に 2 本の斜線を引き，大小 5 個の円を入れる。直角三角形の鈎（直角を挟む 2 辺のうち，短い方の辺）の長さが 1 寸のとき，大円の直径はいかほどか。

正方形の中心を原点とし，一辺の長さを 2a とする。
斜線と正方形の一辺が交差する座標を (b, -a)
大円の半径と中心座標を r1, (x1, a - r1)
小円の半径と中心座標を r2, (a - r2, r2 - a)
として，以下の連立方程式を解く

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, k:: positive,
     r1::positive, x1::positive, r2::positive

eq1 = (a - b) - k
eq2 = x1^2 + (a - r1)^2 - (r1 - r2)^2
eq3 = 2a + k - sqrt(4a^2 + k^2) - 2r2
eq4 = a*(a + b) - r1*sqrt(4a^2 + k^2)  # 中央の平行四辺形の面積を二通りの表現で
eq5 = distance(b, -a, a, a, x1, a - r1) - r1^2|>simplify|>numerator;
eq5 = a^3 - a^2*r1 - 2*a^2*x1 + a*b*r1 - a*r1^2 + a*r1*x1 + a*x1^2 - b*r1*x1;
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, r1, x1, r2))

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 v, r.f_converged
end;

function H(u)
   (a, b, r1, x1, r2) = u
   return [
a - b - k,  # eq1
x1^2 + (a - r1)^2 - (r1 - r2)^2,  # eq2
2*a + k - 2*r2 - sqrt(4*a^2 + k^2),  # eq3
a*(a + b) - r1*sqrt(4*a^2 + k^2),  # eq4
a^3 - a^2*r1 - 2*a^2*x1 + a*b*r1 - a*r1^2 + a*r1*x1 + a*x1^2 - b*r1*x1,  # eq5
   ]
end;

k = 1
iniv = BigFloat[1.95, 0.945, 1.3, 0.236, 0.532]
res = nls(H, ini=iniv)

  (BigFloat[1.550697252562826473616653568524307307723291771022322587402262886914361142019489, 0.5506972525628264736166535685243073077232917710223225874022628869143611420202746, 0.9999999999999999999999999999999999999999999999667428781642976015715212111819139, 0.1775643993863705460811874612232267136362404546140429319571589691343659605198328, 0.4213839097383412741278163330894292515538227879203389202145579779541173446366729], true)

大円の直径は「鈎」の 2 倍である。

鈎 = k = 1 のときのパラメータの値。
a = 1.5507;  b = 0.550697;  r1 = 1;  x1 = 0.177564;  r2 =0.421384
鈎 = 1;  大円の直径 = 2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r1, x1, r2) = res[1]
   k = a - b
   @printf("a = %g;  b = %g;  r1 = %g;  x1 = %g;  r2 =%g\n", a, b, r1, x1, r2)
   @printf("鈎 = %g;  大円の直径 = %g\n", k, 2r1)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle(x1, a - r1, r1)
   circle(-x1, r1 - a, r1)
   circle(0, 0, r2, :green)
   circle(a - r2, r2 - a, r2, :green)
   circle(r2 - a, a - r2, r2, :green)
   segment(b, -a, a, a, :magenta)
   segment(-a, -a, -b, a, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, a - r1, "大円:r1,(x1,a-r1)", :red, :center, :bottom, delta=delta)
       point(a - r2, r2 - a, "小円:r2\n(a-r2,r2-a)", :green, :center, delta=-delta)
       point(a, 0, "a ", :blue, :right)
       point(0, -a, "-a ", :blue, :right, :bottom, delta=delta)
   end
end;

算額（その565）

算額（その565）

群馬の算額 51-3 稲荷神社 文政11年
http://takasakiwasan.web.fc2.com/gunnsann/g051-3.html

直角三角形内に全円と正方形が内接している。この直角三角形の弦を一辺とし，それと直角をなす辺の長さを前述の正方形の一辺の長さとする直角三角形を作る。この直角三角形に内接する円の直径が 5 寸，全円の直径が 7 寸のとき，もとの直角三角形の三辺（鈎，股，弦）の和を求めよ。

元の直角三角形の三辺を「鈎」，「股」，「弦」とする。
全円の半径と中心座標を r1, (r1, r1)
正方形の一辺の長さを a
上の直角三角形に内接する円の半径と中心座標を r2, (x2, y2)，上の頂点を (x, y)
として，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     鈎2::positive, 股2::positive, 弦2::positive,
     r1::positive, r2::positive, a::positive,
     x::positive, y::positive,
     x2::positive, y2::positive

鈎2 = a  # sqrt(x^2 +(y - 鈎)^2)
x = sqrt(a^2 - (y - 鈎)^2)
弦2 = sqrt((股 - x)^2 + y^2)
股2 = 弦

eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = (鈎 - a)/a - 鈎/股
eq4 = 鈎2^2 + 股2^2 - 弦2^2
eq5 = 鈎2 + 股2 - 弦2 - 2r2
eq6 = distance(0, 鈎, 股, 0, x2, y2) - r2^2 |> numerator;
eq7 = distance(x, y, 股, 0, x2, y2) - r2^2 |> numerator;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (鈎, 股, 弦, a, y, x2, y2))

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 v, r.f_converged
end;

function H(u)
   (鈎, 股, 弦, a, y, x2, y2) = u
   return [
       -弦^2 + 股^2 + 鈎^2,  # eq1
       -2*r1 - 弦 + 股 + 鈎,  # eq2
       -鈎/股 + (-a + 鈎)/a,  # eq3
       a^2 - y^2 + 弦^2 - (股 - sqrt(a^2 - (y - 鈎)^2))^2,  # eq4
       a - 2*r2 + 弦 - sqrt(y^2 + (股 - sqrt(a^2 - (y - 鈎)^2))^2),  # eq5
       -r2^2 + (x2 - 股*(x2*股 - y2*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y2 - 鈎*(-x2*股 + y2*鈎 + 股^2)/(股^2 + 鈎^2))^2,  # eq6
       -r2^2 + (x2 - (x2*(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) - y*(x2*y - y*股 + y2*股 - y2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2)))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2))^2 + (-y*(-x2*股 + x2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) + y*y2 + 股^2 - 股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) + y2)^2,  # eq7
   ]
end;

(r1, r2) = (7, 5) .// 2
iniv = BigFloat[9.41379, 16.89655, 19.7931, 6.27586, 15.2069, 3.45, 10.55]
res = nls(H, ini=iniv)

  (BigFloat[10.49999999999999999999999999999999999999999129919479365547439836695548116931581, 14.00000000000000000000000000000000000000001100960233399594221790231656113690163, 17.5000000000000000000000000000000000000000023087971276514166162692720423047526, 5.99999999999999999999999999999999999999999995377737351386752851972765331489346, 15.29999999999999999999999999999999999999999475386066794051341170073747364389192, 3.500000000000000000000000000000000000000000057743557772276350570950054497150227, 10.99999999999999999999999999999999999999999392286748897320175788883023235161591], true)

元の直角三角形の鈎，股，弦の和は 42 寸である。

その他のパラメータは以下の通り。
鈎 = 10.5;  股 = 14;  弦 = 17.5;  a = 6;  x = 3.6;  y = 15.3; , x2 = 3.5;  y2 =11

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (7, 5) .// 2
   (鈎, 股, 弦, a, y, x2, y2) = Float64.(res[1])
   x = sqrt(a^2 - (y -鈎)^2)
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  a = %g;  x = %g;  y = %g; , x2 = %g;  y2 =%g\n", 鈎, 股, 弦, a, x, y, x2, y2)
   @printf("元の直角三角形の鈎，股，弦の和 = %g\n", 鈎 + 股 + 弦)
   plot([0, 0, 股, x, 0, 股], [鈎, 0, 0, y, 鈎, 0], color=:blue, lw=0.5)
   plot!([0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       # hline!([0], color=:black, lw=0.5)
       # vline!([0], color=:black, lw=0.5)
       point(x, y, " (x,y)", :blue, :left, :vcenter)
       point(0, 鈎, " 鈎", :blue, :left, :vcenter, delta=delta/2)
       point(股, 0, "股", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(0, a, " a", :magenta, :left, :bottom, delta=delta/2)
       point(r1, r1, "全円:r1,(r1,r1)", :red, :center, :top, delta=-delta/2)
       point(x2, y2, "円:r2,(x2,y2)", :green, :center, :top, delta=-delta/2)
   end
end;