裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その692)

2024年02月11日 | Julia

算額(その692)

神壁算法 關流藤田貞資門人 龜井隠岐守家士 堀田人助泉尹 天明八年
藤田貞資(1789):神壁算法巻上

http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

直径 97 寸 5 分の大円の中に累円と,大円と甲円と累円に挟まれる円がある。
末円の直径が 1 分のとき,初円から末円まで何個の挟円があるか(つまり,末円は初円から数えて何番目か)。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (x1, y1); x1 = 0, y1 = r1
乙円の半径と中心座標を r2, (x2, y2); x2 = r0 - r2, y2 = 0
乙円以降 i 番目の円の半径と中心座標を ri, (xi, yi)
初円の半径と中心座標を R1, X1, Y1
初円以降 i 番目の挟円の半径と中心座標を Ri, (Xi, Yi)
とおき,順次パラメータを決定していく。

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

using SymPy

@syms r0::positive, r1::positive,
     r2::positive, x2::positive, y2::real,
     r3::positive, x3::positive, y3::positive;

まず,乙円については,半径 r2 = r0/3 になる。

r1 = r0/2
x2 = r0 - r2
y2 = 0
eq1 = x2^2 + r1^2 - (r1 + r2)^2
res2 = solve(eq1, r2)
res2[1] |> println

   r0/3

次に,丙円については,半径,中心座標は (r0/6, 2*r0/3, r0/2) になる。

r2 = r0/3
x2 = r0 - r2
y2 = 0
eq2 = x3^2 + y3^2 - (r0 - r3)^2
eq3 = x3^2 + (r1 - y3)^2 - (r1 + r3)^2
eq4 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
res3 = solve([eq2, eq3, eq4], (r3, x3, y3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/6, 2*r0/3, r0/2)

更に,丁円については,半径,中心座標は (r0/11, 6*r0/11, 8*r0/11) になる。

@syms r4, x4, y4
(r3, x3, y3) = res3[1]
eq5 = x4^2 + y4^2 - (r0 - r4)^2
eq6 = x4^2 + (r1 - y4)^2 - (r1 + r4)^2
eq7 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
res4 = solve([eq5, eq6, eq7], (r4, x4, y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/11, 6*r0/11, 8*r0/11)
    (r0/3, 2*r0/3, 0)

@syms r5, x5, y5
(r4, x4, y4) = res4[1]
eq8 = x5^2 + y5^2 - (r0 - r5)^2
eq9 = x5^2 + (r1 - y5)^2 - (r1 + r5)^2
eq10 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
res5 = solve([eq8, eq9, eq10], (r5, x5, y5))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/18, 4*r0/9, 5*r0/6)
    (r0/6, 2*r0/3, r0/2)

@syms r6, x6, y6
(r5, x5, y5) = res5[1]
eq11 = x6^2 + y6^2 - (r0 - r6)^2
eq12 = x6^2 + (r1 - y6)^2 - (r1 + r6)^2
eq13 = (x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2
res6 = solve([eq11, eq12, eq13], (r6, x6, y6))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/27, 10*r0/27, 8*r0/9)
    (r0/11, 6*r0/11, 8*r0/11)

@syms r7, x7, y7
(r6, x6, y6) = res6[1]
eq14 = x7^2 + y7^2 - (r0 - r7)^2
eq15 = x7^2 + (r1 - y7)^2 - (r1 + r7)^2
eq16 = (x6 - x7)^2 + (y6 - y7)^2 - (r6 + r7)^2
res7= solve([eq14, eq15, eq16], (r7, x7, y7))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/38, 6*r0/19, 35*r0/38)
    (r0/18, 4*r0/9, 5*r0/6)

@syms r8, x8, y8
(r7, x7, y7) = res7[1]
eq17 = x8^2 + y8^2 - (r0 - r8)^2
eq18 = x8^2 + (r1 - y8)^2 - (r1 + r8)^2
eq19 = (x7 - x8)^2 + (y7 - y8)^2 - (r7 + r8)^2
res8= solve([eq17, eq18, eq19], (r8, x8, y8))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0/51, 14*r0/51, 16*r0/17)
    (r0/27, 10*r0/27, 8*r0/9)

累円のパラメータ,半径,中心座標を表す数列には規則性があり,以下のような一般項になる。
甲円の半径 r1,中心座標 (x1, y1) は r0/2, (0, r1)
乙円の半径 r2,中心座標 (x2, y2) は r0/3, (r0 - r2, 0)
丙円以降は i ≧ 3 として,ri, (xi, yi) は以下のようになる。
ri = r0/(i^2 - 2i + 3)
xi = ri*(2i - 2)
yi = ri*(i^2 - 2i)

次に,互いに接する,外円と 2 個の累円に挟まれる円のパラメータを計算する。

@syms R1, X1
X1 = r0 - 2r2 - R1
eq20 = X1^2 + r1^2 - (r1 + R1)^2
res11 = solve(eq20, R1)
res11[1] |> println

   r0/15

R1 = r0/15
X1 = r0 - 2r2 - R1
(R1, X1)

   (r0/15, 4*r0/15)

@syms R2, X2, Y2
(x1, y1, x2, y2) = (0, r1, r0 - r2, 0)
(R1, X1) = (r0/15, 4*r0/15)
eq21 =  X2^2       + (Y2 - y1)^2 - (R2 + r1)^2
eq22 = (X2 - x2)^2 + (Y2 - y2)^2 - (R2 + r2)^2
eq23 = (X2 - x3)^2 + (Y2 - y3)^2 - (R2 + r3)^2
res12 = solve([eq21, eq22, eq23], (R2, X2, Y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-r0, 0, 0)
    (r0/23, 12*r0/23, 8*r0/23)

@syms R3, X3, Y3
(R2, X2, Y2) = res12[2]
eq24 =  X3^2       + (Y3 - y1)^2 - (R3 + r1)^2
eq25 = (X3 - x3)^2 + (Y3 - y3)^2 - (R3 + r3)^2
eq26 = (X3 - x4)^2 + (Y3 - y4)^2 - (R3 + r4)^2
res13 = solve([eq24, eq25, eq26], (R3, X3, Y3))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-r0, 0, 0)
    (r0/39, 20*r0/39, 8*r0/13)

@syms R4, X4, Y4
#(x1, y1, x2, y2) = (0, r1, r0 - r2, 0)
(R3, X3, Y3) = res13[2]
eq27 =  X4^2       + (Y4 - y1)^2 - (R4 + r1)^2
eq28 = (X4 - x4)^2 + (Y4 - y4)^2 - (R4 + r4)^2
eq29 = (X4 - x5)^2 + (Y4 - y5)^2 - (R4 + r5)^2
res14 = solve([eq27, eq28, eq29], (R4, X4, Y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-r0, 0, 0)
    (r0/63, 4*r0/9, 16*r0/21)

@syms R5, X5, Y5
#(x1, y1, x2, y2) = (0, r1, r0 - r2, 0)
(R4, X4, Y4) = res14[2]
eq30 =  X5^2       + (Y5 - y1)^2 - (R5 + r1)^2
eq31 = (X5 - x5)^2 + (Y5 - y5)^2 - (R5 + r5)^2
eq32 = (X5 - x6)^2 + (Y5 - y6)^2 - (R5 + r6)^2
res14 = solve([eq30, eq31, eq32], (R5, X5, Y5))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-r0, 0, 0)
    (r0/95, 36*r0/95, 16*r0/19)

挟円のパラメータは
初円は
R1 = r0/15, X1 = 4*r0/15, Y1 = 0
二円以降 i ≧ 2 においては
半径 Ri = r0/(4i^2 - 4i + 15)
x座標 Xi = Ri * (8i - 4)
y座標 Yi = Ri * (4i^2 - 4i)

問は,「末円の径が一分であるとき,それは何番目の挟円か」ということである。
以下の方程式を解けば,16番目の挟円であることがわかる。

@syms i::positive
r0 = 975/20  # 直径 = 97.5
eq = r0/(4i^2 - 4i + 15) - 1//20  # 直径 = 0.1 寸
solve(eq, i)[1] |> println

   16.0000000000000

一般的な解を求めるには,末円の半径を Rn とすれば

@syms i::positive, r0::positive, Rn::positive
eq = r0/(4i^2 - 4i + 15) - Rn
res = solve(eq, i)[2]
res |> println

   1/2 + sqrt(-14*Rn + r0)/(2*sqrt(Rn))

SymPy ではこれ以上簡単にできないが,手計算で,(1 + sqrt(r0/Rn - 14))/2 となる。

外円の半径を末円の半径で割り,14を引いて,開平し,1を加えて,2で割る」となり,「術」と一致する。
なお,当たり前であるが,術は,半径を直径と読み替えても成り立つ。

ちなみに,20番目の挟円の半径は 0.031759 ほどであり,数値を丸めるとたしかに20番目である事がわかる。

(sqrt(97.5/(2*0.031759) - 14) + 1)/2

   19.99998688035496

それぞれの円のパラメータは以下通り。

   番号:   累円 xi,   累円 yi,   累円 ri;    挟円 Xi,   挟円 Yi,    挟円Ri
      1:  0.000000, 24.375000, 24.375000;  13.000000,  0.000000,  3.250000
      2: 32.500000,  0.000000, 16.250000;  25.434783, 16.956522,  2.119565
      3: 32.500000, 24.375000,  8.125000;  25.000000, 30.000000,  1.250000
      4: 26.590909, 35.454545,  4.431818;  21.666667, 37.142857,  0.773810
      5: 21.666667, 40.625000,  2.708333;  18.473684, 41.052632,  0.513158
      6: 18.055556, 43.333333,  1.805556;  15.888889, 43.333333,  0.361111
      7: 15.394737, 44.901316,  1.282895;  13.852459, 44.754098,  0.266393
      8: 13.382353, 45.882353,  0.955882;  12.238494, 45.690377,  0.203975
      9: 11.818182, 46.534091,  0.738636;  10.940594, 46.336634,  0.160891
     10: 10.572289, 46.987952,  0.587349;   9.880000, 46.800000,  0.130000
     11:  9.558824, 47.316176,  0.477941;   9.000000, 47.142857,  0.107143
     12:  8.719512, 47.560976,  0.396341;   8.259669, 47.403315,  0.089779
     13:  8.013699, 47.748288,  0.333904;   7.629108, 47.605634,  0.076291
     14:  7.412281, 47.894737,  0.285088;   7.086137, 47.765814,  0.065612
     15:  6.893939, 48.011364,  0.246212;   6.614035, 47.894737,  0.057018
     16:  6.442731, 48.105727,  0.214758;   6.200000, 48.000000,  0.050000
     17:  6.046512, 48.183140,  0.188953;   5.834089, 48.087035,  0.044198
     18:  5.695876, 48.247423,  0.167526;   5.508475, 48.159806,  0.039346
     19:  5.383436, 48.301380,  0.149540;   5.216920, 48.221258,  0.035249
     20:  5.103306, 48.347107,  0.134298;   4.954397, 48.273616,  0.031759

function circle4f(x, y, r, color=:red)
  circlef(x, y, r, color)
  circlef(x, -y, r, color)
  circlef(-x, y, r, color)
  circlef(-x, -y, r, color)
end;

function draw(more=false)
   pyplot(size=(800, 800), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 20
   names = Char["甲乙丙丁戊己庚辛壬癸"...]
   r = Vector{Float64}(undef, n)
   x = Vector{Float64}(undef, n)
   y = Vector{Float64}(undef, n)
   r0 = 975//20
   r1 = r0/2
   r2 = r0/3
   x2 = r0 - r2
   (r[1], x[1], y[1]) = (r1, 0, r1)
   (r[2], x[2], y[2]) = (r2, x2, 0)
   for i in 3:n
       factor = r0/(i^2 - 2i + 3)
       r[i] = factor
       x[i] = factor*(2i - 2)
       y[i] = factor*(i^2 - 2i)
   end
   R = Vector{Float64}(undef, n)
   X = Vector{Float64}(undef, n)
   Y = Vector{Float64}(undef, n)
   (R[1], X[1], Y[1]) = (r0/15, 4*r0/15, 0)
   for i in 2:n
       factor = r0/(4i^2 - 4i + 15)
       R[i] = factor
       X[i] = factor*(8i - 4)
       Y[i] = factor*(4i^2 - 4i)
   end
   plot()
   circle(0, 0, r0)
   @printf("%4s: %9s, %9s, %9s;  %9s, %9s, %9s\n", "番号", "累円 xi", "累円 yi", "累円 ri", "挟円 Xi", "挟円 Yi", "挟円Ri")
   for i in 1:n
       @printf("%4i: %9.6f, %9.6f, %9.6f;  %9.6f, %9.6f, %9.6f\n", i, x[i], y[i], r[i], X[i], Y[i], R[i])
       circle4(x[i], y[i], r[i], i)
       circle4f(X[i], Y[i], R[i], i)
       i < 11 && point(x[i], y[i], names[i], :black, :center, :vcenter, mark=false)
       i < 9 && point(X[i]-R[i]-1, Y[i], string(i), :black, :center, :vcenter, mark=false)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(X[1], Y[1], "初円", :white, :left, :bottom, delta=delta/2)
       segment(X[16], Y[16], X[16], Y[16] - 1)
       point(X[16], Y[16] - 1, "末円", :red, :center, :top, delta=-delta/4, mark=false)
       plot!(xlims=(-2, 50), ylims=(-R[1], 50))
   end
end;

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村