算額（その811）

算額（その811）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

長方形に楕円と円が入っている。長方形の長辺と短辺が 16 寸，9 寸のとき，円の直径はいかほどか。

算法助術の公式90 そのものである。
楕円の長径，短径を a，b，円の直径を d とすれば，
d^2 +ab -2d√ab - 2d(a + b) = 0
が成り立つ。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, d::positive
formula90 = d^2 + a*b - 2d*sqrt(a*b) - 2d*(a + b);

d について解けば，a, b を与えて d を求める公式になる。

d = solve(formula90, d)[1]
d |> println

   sqrt(a)*sqrt(b) + a + b - sqrt(2*a^(3/2)*sqrt(b) + 2*sqrt(a)*b^(3/2) + a^2 + 2*a*b + b^2)

楕円の長径，短径を与えれば円の直径が求まる。

d(a => 16, b => 9) |> println

   2

関数にしてみる。

func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);

長径，短径が 16, 9 のとき，円の直径は 2 である。

func(16, 9)

   2.0

この公式，関数を用いるとき，長径，短径をあたえれば直径が得られ，長半径，短半径を与えれば半径が得られる。

func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (16, 9) .// 2
   r = func(a, b)
   @printf("円の直径 = %g\n", 2r)
   # 図形描画システムと統一を図るため長半径，短半径，半径にする
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   circle4(a - r, b - r, r, :green)
   if more == true
       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(a - r, b - r, "(a-r,b-r) ", :black, :right, :bottom, delta=delta)
       point(a, 0, " a", :red, :left, :bottom, delta=delta)
       point(0, b, " b", :red, :left, :bottom, delta=delta)        
   end
end;

算額（その810）

算額（その810）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直角三角形内に，大円，中円，小円を入れる。中円，小円の直径がそれぞれ 4635 寸，2060 寸のとき，大円の直径はいかほどか。

直角三角形の中に3個の円を入れる類題は多い。算額（その740），算額（その453），算額（その188）

大円の半径と中心座標を r1, (r1, r2)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (r3, y3)
直角を挟む短い方の辺と長いほうの辺の長さを a, b
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, x2::positive,
     r3::positive, y3::positive, a::positive,
     b::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = r1/(b - r1) - r2/(b - x2)
eq4 = r1/(a - r1) - r3/(a - y3)
eq5 = a + b - sqrt(a^2 + b^2) - 2r1;
#res = solve([eq1, eq2, eq3, eq4, eq5], (r1, x2, y3, a, b))

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)
   (r1, x2, y3, a, b) = u
   return [
       (-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq1
       (-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,  # eq2
       r1/(b - r1) - r2/(b - x2),  # eq3
       r1/(a - r1) - r3/(a - y3),  # eq4
       a + b - 2*r1 - sqrt(a^2 + b^2),  # eq5
   ]
end;

(r2, r3) = (4635, 2060) .// 2
iniv = BigFloat[3500, 9200, 7300, 8900, 20300]
res = nls(H, ini=iniv)

   ([3515.5003002014273, 9224.150559489277, 7321.26713972666, 8898.390201990349, 20267.383999142992], true)

大円の直径は 7031 寸である。

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

r1 = 3515.5;  x2 = 9224.15;  y3 = 7321.27;  a = 8898.39;  b = 20267.4

function draw(more)
    pyplot(size=(700, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (4635, 2060) .// 2
   (r1, x2, y3, a, b) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   @printf("r1 = %g;  x2 = %g;  y3 = %g;  a = %g;  b = %g\n", r1, x2, y3, a, b)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(r3, y3, r3, :magenta)
   if more == true
       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(r1, r1, "大円:r1\n(r1,r1)", :red, :center, :bottom, delta=3delta)
       point(x2, r2, " 中円:r2\n(x2,r2)", :green, :center, delta=-delta)
       point(r3, y3, "小円:r3,(r3,y3) ", :blue, :left, delta=-3delta)
       point(0, a, " a", :black, :left, :bottom)
       point(b, 0, " b", :black, :left, :bottom)
   end
end;

算額（その809）

算額（その809）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

不等辺三角形内に全円，大円，中円，小円を入れる。大円，中円，小円の直径がそれぞれ 9 寸，4 寸，1 寸のとき，全円の直径はいかほどか。

三角形の底辺の長さを a，頂点の座標を (x0, y0) 
全円の半径と中心座標を r1, (x1, r1)
大円の半径と中心座標を r2, (x2, r2)
中円の半径と中心座標を r3, (x3, r3)
小円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     x2::positive, r3::positive, x3::positive,
     r4::positive, x4::positive, y4::positive,
     a::positive, x0::positive, y0::positive,
     d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - x4)^2 + (r1 - y4)^2 - (r1 + r4)^2
# eq4 = r1/(a - x1) - r3/(a - x3)
eq5 = r2/x2 - r1/x1
eq6 = (a + sqrt(x0^2 + y0^2) + sqrt((a - x0)^2 + y0^2))*r1 - y0*a
eq7 = dist(0, 0, x0, y0, x2, r2) - r2^2
eq8 = dist(0, 0, x0, y0, x4, y4) - r4^2
eq9 = dist(a, 0, x0, y0, x3, r3) - r3^2
eq10 = dist(a, 0, x0, y0, x4, y4) - r4^2
eq7 = numerator(apart(eq7, d));
eq8 = numerator(apart(eq8, d));
eq9 = numerator(apart(eq9, d));
eq10 = numerator(apart(eq10, d));
# res = solve([eq1, eq2, eq3, eq5, eq6, eq7, eq8, eq9, eq10],
#     (r1, x1, x2, x3, x4, y4, a, x0, y0))

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)
   (r1, x1, x2, x3, x4, y4, a, x0, y0) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq2
       -(r1 + r4)^2 + (r1 - y4)^2 + (x1 - x4)^2,  # eq3
       -r1/x1 + r2/x2,  # eq5
       -a*y0 + r1*(a + sqrt(x0^2 + y0^2) + sqrt(y0^2 + (a - x0)^2)),  # eq6
       -r2^2*y0^2 - 2*r2*x0*x2*y0 + x2^2*y0^2,  # eq7
       -r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2,  # eq8
       -2*a^2*r3*y0 + a^2*y0^2 + 2*a*r3*x0*y0 + 2*a*r3*x3*y0 - 2*a*x3*y0^2 - r3^2*y0^2 - 2*r3*x0*x3*y0 + x3^2*y0^2,  # eq9
       -a^2*r4^2 + a^2*y0^2 - 2*a^2*y0*y4 + a^2*y4^2 + 2*a*r4^2*x0 + 2*a*x0*y0*y4 - 2*a*x0*y4^2 - 2*a*x4*y0^2 + 2*a*x4*y0*y4 - r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2,  # eq10
   ]
end;

(r2, r3, r4) = (9, 4, 1) .// 2
iniv = BigFloat[5.5, 54.7, 44.8, 61.4, 57, 11.1, 65.1, 57.21, 11.6]
res = nls(H, ini=iniv)

   ([5.5, 54.7243090408641, 44.7744346697979, 61.357558621574896, 56.97961389830577, 11.06, 65.14798695340964, 57.20514438404994, 11.616], true)

全円の直径は 11 寸である。

算額にはよくあることであるが，実際の図と算額に描かれている図には大きな違いがある。

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

r1 = 5.5;  x1 = 54.7243;  x2 = 44.7744;  x3 = 61.3576;  x4 = 56.9796;  y4 = 11.06;  a = 65.148;  x0 = 57.2051;  y0 = 11.616

function draw(more)
    pyplot(size=(800, 170), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (9, 4, 1) .// 2
   (r1, x1, x2, x3, x4, y4, a, x0, y0) = res[1]
   @printf("全円の直径 = %g\n", 2r1)
   @printf("r1 = %g;  x1 = %g;  x2 = %g;  x3 = %g;  x4 = %g;  y4 = %g;  a = %g;  x0 = %g;  y0 = %g\n", r1, x1, x2, x3, x4, y4, a, x0, y0)
   plot([0, a, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(x3, r3, r3, :magenta)
   circle(x4, y4, r4, :blue)
   if more == true
       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(x1, r1, "全円:r1\n(x1,r1)", :red, :center, :bottom, delta=3delta)
       point(x2, r2, "大円:r2\n(x2,r2)", :green, :center, delta=-3delta)
       point(x3, r3, "中円:r3,(x3,r3) ", :black, :right, :vcenter)
       point(x4, y4, "小円:r4,(x4,y4) ", :blue, :right, :bottom, delta=2delta)
       point(x0, y0, " (x0,y0)", :black, :left, :vcenter)
       point(a, 0, " a", :black, :left, :bottom)
   end
end;