桁落ちさせてみた - SELECT * FROM life;

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$
桁落ち恐ろしや