前回の記事
「FORTRAN77数値計算プログラミング」のプログラムをCとPythonに移植してみる（その1）
FORTRAN77数値計算プログラミング、1.2のコードより。
#台形則における誤差の累積
##丸め誤差の累積
P.5の終わり〜P.6プログラムの手前まで引用
丸め誤差の累積を示す一つの例として、積分
I = \int_{0}^{1}\frac{1}{1+x^2}dx = \frac{\pi}{4} = 0.7853982\tag{1.13}
>を、きざみ幅$h$を次第に小さくしながら、**台形則**
>```math
I_h = h\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)\biggr], \quad h=\frac{b-a}{n}\tag{1.14}
a=0, \quad b=1, \quad f(x)=\frac{1}{1+x^2} \tag{1.15}
によって数値積分してみる。
##Fortranのコード
(1.14),(1.15)の数値積分を行うプログラム。
* Accumulation of round off error
PROGRAM STRAP
PARAMETER (ONE = 1.0)
*
EV = ATAN(ONE)
*
WRITE (*,2000)
2000 FORMAT (' ---- TRAP ----')
*
DO 10 K = 1, 13
*
N = 2**K
H = 1.0 / N
*
S = 0.5 * (1.0 + 0.5)
DO 20 I = 1, N - 1
S = S + 1 / (1 + (H * I)**2)
20 CONTINUE
S = H * S
*
ERR = S - EV
IF (ERR .NE. 0.0) THEN
ALERR = ALOG10 (ABS(ERR))
ELSE
ALERR = -9.9
END IF
*
WRITE (*,2001) N, H, S, ERR, ALERR
2001 FORMAT (' N=',I6,' H=',1PE9.2,' S=',E13.6,
$ ' ERR=',E8.1,' L(ERR)=',0PF4.1)
*
10 CONTINUE
*
END
実行結果は
---- TRAP ----
N= 2 H= 5.00E-01 S= 7.750000E-01 ERR=-1.0E-02 L(ERR)=-2.0
N= 4 H= 2.50E-01 S= 7.827941E-01 ERR=-2.6E-03 L(ERR)=-2.6
N= 8 H= 1.25E-01 S= 7.847471E-01 ERR=-6.5E-04 L(ERR)=-3.2
N= 16 H= 6.25E-02 S= 7.852354E-01 ERR=-1.6E-04 L(ERR)=-3.8
N= 32 H= 3.12E-02 S= 7.853574E-01 ERR=-4.1E-05 L(ERR)=-4.4
N= 64 H= 1.56E-02 S= 7.853879E-01 ERR=-1.0E-05 L(ERR)=-5.0
N= 128 H= 7.81E-03 S= 7.853957E-01 ERR=-2.4E-06 L(ERR)=-5.6
N= 256 H= 3.91E-03 S= 7.853976E-01 ERR=-5.4E-07 L(ERR)=-6.3
N= 512 H= 1.95E-03 S= 7.853978E-01 ERR=-4.2E-07 L(ERR)=-6.4
N= 1024 H= 9.77E-04 S= 7.853981E-01 ERR=-6.0E-08 L(ERR)=-7.2
N= 2048 H= 4.88E-04 S= 7.853982E-01 ERR= 6.0E-08 L(ERR)=-7.2
N= 4096 H= 2.44E-04 S= 7.853983E-01 ERR= 1.2E-07 L(ERR)=-6.9
N= 8192 H= 1.22E-04 S= 7.853981E-01 ERR=-6.0E-08 L(ERR)=-7.2
になります。あら、$N=128$あたりからERRとL(ERR)の出力が本と違います…。
##コードの理解
そもそもなんでこのプログラムの名前がSTRAPなんでしょうか。と疑問はさておき、積分てなんでしたっけねー、としばし固まる私。
気を取り直して、本に書いていることとコードを確認していきます。
本に書いてある、積分の式を台形則によって数値積分するっていうのは、元の積分で表されるグラフの面積のところにたくさんの台形をおいて、その台形の面積の和を積分グラフの面積の近似値にするということ。
ので、プログラムがやってることは(1.14)(1.15)の式を、$h$をだんだん小さくしながら実行するという作業。それが（1.13）で求められる答えとどのくらい誤差があるか、というのが主題。
流れ図描きました。
まず、※1でやっているのは $arctan 1.0$、逆三角関数アークタンジェントを使っています。アークタンジェントを定積分で表すと、
arctan\,x = \int_0^x \frac 1 {z^2 + 1}\,dz
という式になります。ここで$x$を$1$、$z$を$x$とおくと（1.13）の式と同じになるので、ここでEVには(1.13）と同じ値、つまり$I$が代入されます。
次に※2では、$N$ に $2^K$ を代入していて、$N$は式(1.14)で $$\sum_{j=1}^{n-1}$$
に使われている順序数、まあ、ざっくり言うとだんだん増えていく値ってことで。
ここで使っている$K$はループ1のカウンタで、誤差の挙動を2の13乗まで試しましょ、というのがこのループ1というわけ。
そして※3ですが、これはもともと（1.14）で $$h=\frac{b-a}{n}$$ と表されていて、（1.15）で$a=0, \quad b=1$ となっているので、代入すると$$h=\frac{1}{n}$$ そのままですね。
で、（1.14）の１つ目の式$$I_h = h\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)\biggr]$$ にまず２つ目の式$$h=\frac{1}{n}$$ を当てはめると
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+j\frac{1}{n})+\frac{1}{2}f(b)\biggr]
となり、さらに（1.15）から$a=0$と$b=1$で置き換えると
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(0)+\sum_{j=1}^{n-1}f(0+\frac{j}{n})+\frac{1}{2}f(1)\biggr]
そして$$f(x)=\frac{1}{1+x^2}$$ なので、
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}\times\frac{1}{1+0^2}+\sum_{j=1}^{n-1}\frac{1}{1+\big(\frac{j}{n}\big)^2} +\frac{1}{2}\times\frac{1}{1+1^2}\biggr]
から、
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n}{1+\frac{j^2}{n^2}}+\frac{1}{4}\biggr]\tag{r1}
さらに
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}+\frac{1}{4}\biggr]
で、
I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{3}{4}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}\biggr]\tag{r2}
と整理できます。
※４は、$0.5\times(1.0+0.5)$ を $S$に代入していますが、これを$$\frac{1}{2}\times\big(\frac{2}{2}+\frac{1}{2}\big)$$ と考えてみると、$$\frac{1}{2}\times\frac{3}{2}=\frac{3}{4}$$ となって、式（r2）に出てくる値になってますね。
ループ2のカウンタIは $$\sum_{j=1}^{n-1}$$ で書かれている$j$ですね。Jにしといてくれたらわかりやすいのに。
そして※5では、$$S+1\div(1+(H\times I)^2)$$ を $S$に代入していますが、この式の$$1\div(1+(H\times I)^2)$$ の部分を整理していきます。まずは
$$1\div1+{\big(\frac{1}{n}\times j\big)}^2$$ なので
$$\frac{1}{1+{\big(\frac{1}{n}\times j\big)}^2}$$
から
$$\frac{1}{1+\frac{ｊ^2}{n^2}}$$ となって式(r1)に出てくる形になりました。
というわけで、ループ2は（1.14）の
\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)
の部分の計算を行っているものになります。
※6で（1.14）の最初の式が完成して、※7では、（1.13）の式との差分、つまりもともと積分計算で求めた計算結果と、台形則を使った近似値との差を計算しています。
※8では、※7で求めた差が0と等しくなかった場合に、差の絶対値を求め（ABS(ERR））、その値の常用対数を求め（ALOG10()）ています。
常用対数とは、$$x=10^a$$ となる$a$の値のことで、数式では$$a=\log_{10}x$$ と表されるもの。
この値を使って、次のようなグラフがプロットできます。
このグラフが表しているのは、$h^2$に比例して累積誤差が減少していくということ、ただし、誤差が計算機イプシロンに近づくにつれて減少しなくなるということ。
###本との違い
比べると、$S$の値が$N=128$からごく僅かですが増えています。累積誤差の減少が止まる地点が少しあとになってるんですね。
本が使っている計算機（あえてコンピュータと書かないw）がMV4000 AOS/VSということで、16ビットOSなんですよね。
で、私が出したグラフ、本で参考として出ている倍精度浮動小数点数を使って計算したグラフに似ています。あくまでも似ているだけでちょっと違うんですが、本の環境では浮動小数点演算が6桁丸め演算とのこと。
ということで、これは64ビットOS（環境は32ビットかな）で切り捨て演算を行った結果による違いであろうと推測。
識者の方いたら教えて下さい。
追記：ほんとに識者の方（@cure_honey）にコメントで教えていただきました。ありがとうございます。
##Cのコード
#include <stdio.h>
#include <math.h>
int main(void)
{
const float ONE = 1.0f;
float ev = atanf(ONE);
int k, i;
int n;
float h, s, err, alerr;
printf("%14.7E\n", ev);
printf("--- TRAP ---\n");
for(k=1; k<=13; k++)
{
n = pow(2,k);
h = 1.0f / n;
s = 0.5f * (1.0f + 0.5f);
for(i=1; i<=n-1; i++)
{
s = s + 1 / (1+ powf(h*i,2));
}
s = h * s;
err = s - ev;
if(err != 0.0f)
{
alerr = log10(fabsf(err));
}
else
{
alerr = -9.9f;
}
printf("N=%6d H=%9.2E S=%13.6E ERR=%8.1E L(ERR)=%4.1F\n",n, h, s, err, alerr);
}
return 0;
}
実行結果は変わりません。
##Pythonのコード
import numpy as np
import matplotlib.pyplot as plt
ONE = np.array((1.0),dtype=np.float32)
ev = np.arctan(ONE,dtype=np.float32)
h, s, err, alerr = np.array([0.0, 0.0, 0.0, 0.0],dtype=np.float32)
# 計算用 1.0, 0.5, 0.0, -9.9の単精度浮動小数点数
tmp = np.array([1.0, 0.5, 0.0, -9.9],dtype=np.float32)
# グラフ描画用リスト
x = []
y = []
print("--- TRAP ---")
for k in range(13):
n = np.power(2,k+1)
h = tmp[0] / n.astype(np.float32)
s = tmp[1] * (tmp[0] + tmp[1])
for i in range(n-1):
s = s + tmp[0] / (tmp[0] + np.square(h*(i+1),dtype=np.float32))
s = s * h
err = s - ev
if err != tmp[2]:
alerr = np.log10(np.abs(err))
else:
alerr = tmp[3]
print("N=%6d H=%9.2E S=%13.6E ERR=%8.1E L(ERR)=%4.1F" % (n, h, s, err, alerr))
# グラフ用変数セット
x.append(k+1)
y.append(alerr)
# グラフ描画
fig, ax = plt.subplots()
ax.scatter(x,y)
ax.set_xlabel("n (1/2^n)")
ax.set_ylabel("l (10^l)")
plt.show()
出力結果は同じです。あと、グラフ出しました。
##まとめ
Pythonのコードはもっときれいになりそうな気がします。
そして、単精度にするためにnumpy書き散らすのもどうかと思う。Fortranのコードの方を倍精度にするか……いやしかしそれだと実行結果が本と照らし合わせられないし。むう。