読者です 読者をやめる 読者になる 読者になる

桁落ちさせてみた

3回生配当の「数値解析」という授業を取りこぼしていたので、4回生ながら履修することにしました。
で、第一回の授業が今日あって、桁落ちの例が紹介されたので実際にPythonでやってみる。

内容

3a_{n+1}+5a_n-2a_{n-1}=0
a_0=1,a_1=1/3
という漸化式を考えると一般解は
a_n=\frac{1}{3^n}
となるが
\right a_n=\frac{2a_{n-1}-5a_n}{3}
に順次数値を代入していくと引き算があるため、桁落ちが起こり結果がおかしくなってしまう

コード

 # a0 = 1
 first = 1.0
 # a1 = 1/3
 second = 1.0/3.0
 
 # 正しい値はans^nで求められる。
 ans = 1.0/3.0
  
 for i in range(2, 60):
     # an+1を計算
     third = (2*first - 5*second)/3
     # 次のループのために値を移動させる
     first, second = second, third
     # 結果を出力。左から n 計算結果 正しい答え
     print '%s\t%s\t%s'%(i,third, ans**i)

やってみる

2       0.111111111111          0.111111111111
3       0.037037037037          0.037037037037
4       0.0123456790123         0.0123456790123
5       0.00411522633745        0.00411522633745
6       0.00137174211248        0.00137174211248
7       0.000457247370826       0.000457247370828
8       0.000152415790279       0.000152415790276
9       5.08052634194e-05       5.08052634253e-05
10      1.69350878203e-05       1.69350878084e-05
11      5.6450292458e-06        5.64502926948e-06
12      1.88167647051e-06       1.88167642316e-06
13      6.27225379688e-07       6.27225474386e-07
14      2.09075347525e-07       2.09075158129e-07
15      6.96913405842e-08       6.96917193763e-08
16      2.32313307095e-08       2.32305731254e-08
17      7.74200920696e-09       7.74352437514e-09
18      2.58420512807e-09       2.58117479171e-09
19      8.54330924519e-10       8.60391597238e-10
20      2.98918544516e-10       2.86797199079e-10
21      7.13563754855e-11       9.55990663597e-11
22      8.03517372017e-11       3.18663554532e-11
23      -8.63486450125e-11      1.06221184844e-11
24      1.97482233155e-10       3.54070616147e-12
25      -3.86702818601e-10      1.18023538716e-12
26      7.76159519771e-10       3.93411795719e-13
27      -1.55140107869e-09      1.3113726524e-13
28      3.10310814432e-09       4.37124217466e-14
29      -6.206114293e-09        1.45708072489e-14
30      1.24122625845e-08       4.85693574962e-15
31      -2.48245138362e-08      1.61897858321e-15
32      4.96490314501e-08       5.39659527735e-16
33      -9.9298061641e-08       1.79886509245e-16
34      1.98596123702e-07       5.99621697484e-17
35      -3.97192247264e-07      1.99873899161e-17
36      7.94384494574e-07       6.66246330538e-18
37      -1.58876898913e-06      2.22082110179e-18
38      3.17753797827e-06       7.40273700597e-19
39      -6.35507595654e-06      2.46757900199e-19
40      1.27101519131e-05       8.22526333997e-20
41      -2.54203038261e-05      2.74175444666e-20
42      5.08406076523e-05       9.13918148886e-21
43      -0.000101681215305      3.04639382962e-21
44      0.000203362430609       1.01546460987e-21
45      -0.000406724861218      3.38488203291e-22
46      0.000813449722437       1.12829401097e-22
47      -0.00162689944487       3.76098003657e-23
48      0.00325379888975        1.25366001219e-23
49      -0.00650759777949       4.1788667073e-24
50      0.013015195559          1.3929555691e-24
51      -0.026030391118         4.64318523033e-25
52      0.0520607822359         1.54772841011e-25
53      -0.104121564472         5.15909470036e-26
54      0.208243128944          1.71969823345e-26
55      -0.416486257888         5.73232744485e-27
56      0.832972515775          1.91077581495e-27
57      -1.66594503155          6.3692527165e-28
58      3.3318900631            2.12308423883e-28
59      -6.6637801262           7.07694746278e-29

結果

n=15くらいまでは正しく計算できているが、その後で桁落ちの影響でどんどん正しい値から離れていくのが分かる。
見やすくするとa_{59}の時点で
正解a_{59}=\frac{1}{3^{59}}=7.07694746278\times10^{-29}
桁落ちした方a_{59}=-6.6637801262
桁落ち恐ろしや