算額（その2002）

算額（その2002）

(7) 滋賀県 錦織良夫氏宅　p.114
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，楕円，算法助術公式96

楕円の中に 4 個の等円を容れる。楕円の長径，短径が与えられたとき，等円の直径を求めよ。

1.　算法助術の公式96

楕円の長半径，短半径を a, b，等円の半径を r とおく。

算法助術では，楕円の長軸 p = 2a，短軸 q = 2b，円の直径 d = 2r を使っているので，変換する。

include("julia-source.txt");

using SymPy

@syms p, q, d, a, b, r
p = 2a
q = 2b
d = 2r
eq1 = p^2*q^2*(p^2 - q^2)^2 - 4d^2*(p^2 + q^2)*(p^2 - q^2)^2 + 12d^4*(p^2 - q^2)^2 + 4d^4*p^2*q^2 - 12d^6*(p^2 + q^2) + 4d^8;

eq2 = eq1 |> simplify |> factor |> x -> x/256;
eq2 |> println

   a^6*b^2 - 4*a^6*r^2 - 2*a^4*b^4 + 4*a^4*b^2*r^2 + 12*a^4*r^4 + a^2*b^6 + 4*a^2*b^4*r^2 - 20*a^2*b^2*r^4 - 12*a^2*r^6 - 4*b^6*r^2 + 12*b^4*r^4 - 12*b^2*r^6 + 4*r^8

r の 8次式なので r について解くと 8 個の解が得られるが，其のうちの　6 番目が適解である。
式は長いので省略する。

res = solve(eq2, r)[6];

長半径，短半径が 10, 4 のとき，等円の半径は 1.95789337209971 である。

res(a => 10, b => 4).evalf() |>  println

   1.95789337209971

2.　連立方程式を解く

楕円の長半径，短半径を a, b
等円の半径と中心座標を r, (r, r)
楕円と等円の接点座標を (x0, y0)
とおき，以下の連立方程式を解き数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive,
     x0::positive, y0::positive
@syms a, b, x0, y0, r
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/(y0 - r)
eq3 = (x0 - r)^2 + (y0 - r)^2 - r^2;

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)
   (r, x0, y0) = u
   return [
       -1 + y0^2/b^2 + x0^2/a^2,  # eq1
       (-r + x0)/(-r + y0) - b^2*x0/(a^2*y0),  # eq2
       -r^2 + (-r + x0)^2 + (-r + y0)^2,  # eq3
   ]
end;

(a, b) = (10, 4)
iniv = BigFloat[1.943, 2.13, 3.877]
res = nls(H, ini=iniv)

   ([1.9578933720997147, 2.127795322101097, 3.9084009429380644], true)

楕円の長径，短径が 20，8 のとき，等円の直径は 3.9157867441994294 である。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0) = res[1]
   @printf("a = %g;  b = %g;  r = %g;  x0 = %g;  y0 = %g\n", a, b, r, x0, y0)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle4(r, r, r, :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(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b, "b", :green, :center, :bottom, delta=delta)
       point(r, r, "等円:r\n(r,r)", :blue, :center, delta=-delta)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
   end
end;

draw(10, 4, true)

 

算額（その2001）

算額（その2001）

(7) 滋賀県 吉田博氏宅　p.116
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，長方形，斜線2本

長辺 50 間，短辺 30 間の長方形の田がある。上，中，下の（円形の）田と接する2本の斜線を引くとき，中田の直径はいかほどか。

長方形の長辺，短辺を 2a, b
上田の半径と中心座標を r1, (0, b - r1)
中田の半径と中心座標を r2, (a - r2, r2)
下田の半径と中心座標を r3, (0, r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive, r3::positive
eq1 = (a - r2)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq2 = dist2(0, 2r3, a, b, 0, b - r1, r1)
eq3 = dist2(0, 2r3, a, b, a - r2, r2, r2)
res = solve([eq1, eq2, eq3], (r1, r2, r3))[4]

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

2res[1](a => 50/2, b => 30).evalf() |> println
2res[2](a => 50/2, b => 30).evalf() |> println
2res[3](a => 50/2, b => 30).evalf() |> println

   16.6666666666667
   20.0000000000000
   11.2500000000000

上田，中田，下田の直径はそれぞれ 50/3, 20, 45/4 である。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(a^2 + 4b^2)
   u = 2a - 3b
   (r1, r2, r3) =(-a*sqrt(b^2*u^2 + (-a^2 + a*b + a*t + 2b^2 - b*t)^2)/(b*u) + a*(-a^2 + a*b + a*t + 2*b^2 - b*t)/(b*u), a/2 + b - t/2, (a^2 - a*b + a*t + 2b^2 - b*t)/(8b))
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, b, r1, r2, r3)
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:green, lw=0.5)
   circle(0, b - r1, r1)
   circle2(a - r2, r2, r2, :blue)
   circle(0, r3, r3, :magenta)
   segment(0, 2r3, a, b, :orange)
   segment(0, 2r3, -a, b, :orange)
   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(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b, "b", :green, :center, :bottom, delta=delta)
       point(0, b - r1, "上:r1,(0,b-r1)", :red, :center, delta=-delta)
       point(a - r2, r2, "中:r2,(a-r2,r2)", :blue, :center, delta=-delta)
       point(0, r3, "下:r3,(0,r3)", :magenta, :center, delta=-delta/2)
   end
end;

draw(50/2, 30, true)

算額（その2000）

算額（その2000）

(9) 滋賀県新旭町大田 大田神社 慶応4年(1868)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円5個，外円，デカルトの円定理

外円の中に，大円，中円，小円，内円を容れる。大円，中円，小円の直径がそれぞれ 69 寸，46 寸，23 寸のとき，内円，外円の直径はいかほどか。

1.　デカルトの円定理を用いる場合

外円，大円，中円，小円，内円 の半径を r0, r1, r2, r3, r4 とおき，以下により内円，外円の半径を得る。

include("julia-source.txt");

using SymPy
@syms R::positiveb, r1::positive, r2::positive,
     r3::positive, r4::positive
k4 = 1/r1 + 1/r2 + 1/r3 + 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r4 = 1/k4
k0 = 1/r1 + 1/r2 + 1/r3 - 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r0 = -1/k0;

大円，中円，小円の直径がそれぞれ 69 寸，46 寸，23 寸のとき，内円，外円の直径は 6 寸と 138 寸である。
なお，答えの「外円径一二八寸」は誤記。

2r4(r1 => 69//2, r2 => 46//2, r3 => 23//2) |> println
2r0(r1 => 69//2, r2 => 46//2, r3 => 23//2) |> println

   6
   138

2.　2 円の相互関係から中心座標も求める場合

外円の半径と中心座標を R
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (x3, y3)
内円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。

@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4
(r1, r2, r3) = (69, 46, 23) .// 2
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq4 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq5 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2
(R, x2, y2, x3, y3) = solve([eq1, eq2, eq3, eq4, eq5], (R, x2, y2, x3, y3))[4]
#
eq6 = x4^2 + (R - r1 - y4)^2 - (r1 + r4)^2
eq7 = (x4 - x2)^2 + (y4 - y2)^2 - (r2 + r4)^2
eq8 = (x4 - x3)^2 + (y4 - y3)^2 - (r3 + r4)^2
(r4, x4, y4) = solve([eq6, eq7, eq8], (r4, x4, y4))[2]
(R, x2, y2, x3, y3,   r4, x4, y4)

   (69, 46, 0, 46, 69/2, 3, 36, 24)

内円，外円の直径は 138 寸, 6 寸である。
算額に描かれている図とはずいぶん違う。

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

  r1 = 34.5;  r2 = 23;  r3 = 11.5;  R = 69;  x2 = 46;  y2 = 0;  x3 = 46;  y3 = 34.5

function draw(r1, r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x2, y2, x3, y3) = (69, 46, 0, 46, 69/2)
   (r4, x4, y4) = (3, 36, 24)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  R = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n", r1, r2, r3, R, x2, y2, x3, y3)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x2, y2, r2, :magenta)
   circle(x3, y3, r3, :green)
   circle(x4, y4, r4, :orange)
   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(R, 0, " R", :red, :left, :bottom, delta=delta)
       point(0, R, "R", :red, :center, :bottom, delta=delta)
       point(0, R - r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, y2, "中円:r2,(x2,r2)", :magenta, :center, delta=-delta)
       point(x3, y3, "小円:r3\n (x3,r3)", :green, :center, delta=-delta/2)
       point(x4, y4, "内円:r4,(x4,r4)  ", :black, :right, :vcenter)
   end
end;

draw(69/2, 46/2, 23/2, true)

算額（その1199）

算額（その1199）

(9) 滋賀県マキノ町海津 天神社 明治8年(1875)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円5個，半円1個

半円の中に，甲円，乙円，丙円を容れる。甲円，丙円の直径が 18寸，2 寸のとき，乙円の直径はいかほどか。

この図形においては，「和算の心（その003）」に示すように，甲円，乙円，丙円の半径は 1:1/2:1/9 である。さらに甲円の中心は半円の直径上にあり，x1 = 0 である。

半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2); x2 < 0
丙円の半径と中心座標を r3, (x3, r3); x3 < 0
とおく。

この図形においては，「和算の心（その003）」に示すように，算額の図とは全く異なり左右対称で，甲円，乙円，丙円の半径は 1:1/2:1/9 である。さらに甲円の中心は半円の直径上にあり，x1 = 0 である。

include("julia-source.txt");

using SymPy
@syms R::positiveb, r1::positive, x1,
     r2::positive, x2::negative,
     r3::positive, x3::negative
x1 = 0
eq1 = x1^2 + r1^2 - (R - r1)^2
eq2 = x2^2 + r2^2 - (R - r2)^2
eq3 = x3^2 + r3^2 - (R - r3)^2
eq4 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2;

res = solve([eq1, eq2, eq3, eq4], (R, r2, x2, x3))[1]

   (2*r1, r1/2, -sqrt(2)*r1, -2*sqrt(r1)*sqrt(r1 - r3))

甲円，丙円の直径が 18寸，2 寸のとき，乙円の直径は甲円の直径の 1/2 である。すなわち，9 寸である。

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

  r1 = 9;  r3 = 1;  R = 18;  x1 = 0;  r2 = 4.5;  x2 = -12.7279;  x3 = -16.9706

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x1 = 0
   (R, r2, x2, x3) = (2*r1, r1/2, -sqrt(2)*r1, -2*sqrt(r1)*sqrt(r1 - r3))
   @printf("r1 = %g;  r3 = %g;  R = %g;  x1 = %g;  r2 = %g;  x2 = %g;  x3 = %g\n", r1, r3, R, x1, r2, x2, x3)
   plot()
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle2(x2, r2, r2, :magenta)
   circle2(x3, r3, r3, :green)
   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(R, 0, " R", :red, :left, :bottom, delta=delta)
       point(0, R, "R", :red, :center, :bottom, delta=delta)
       point(x1, r1, "甲円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "乙円:r2,(x2,r2)", :magenta, :center, delta=-delta)
       point(x3, r3, " 丙円:r3,(x3,r3)", :green, :left, :vcenter)
   end
end;

draw(18/2, 2/2, true)