□演習 8.2 (擬点法) 擬点法は f(a1)f(b1)<0
となるような近似解の組 x=a1, x=b1 から出発して、漸化式 xk = bkf(ak)−akf(bk)
f(ak)−f(bk) (k= 1,2,3, ...) (ak+1, bk+1) =
( (ak, xk) if f(ak)f(xk)<0 (xk, bk) otherwise
y=f(x)
x b ak xk
k
y
により、次々とより良い近似解x=xk(k = 1,2,3, ...) を求めていこうというものである。
f(xk) = 0 となった時点、あるいは x1, x2, x3, ... が十分に収束したと判断できる時点で、
この計算は終了する。 擬点法を用いて方程式 f(x) =x−cosx= 0 の近似解を小数点以 下15桁まで求めるCプログラムを作成せよ。但し、ここではa1 = 0, b1 = 1 とせよ。
8.2. 数値積分 —Simpsonの公式— 123
'
&
$
% 補足：
• n は1以上の整数で、十分大きく選ぶ。
このnで指定された個数に積分区間[a, b]は等分割され、各々の小区間にお いてf(x)が2次多項式で近似され定積分の近似値が求められることになる。
• 特に f(x) = 4
1 +x2 の場合は、
f′(x) = −8x (1 +x2)2
f′′(x) = 8(3x2−1) (1 +x2)3
f′′′(x) = 96x(1−x2) (1 +x2)4
f(4)(x) = 96(5x4−10x2+ 1) (1 +x2)5
f(5)(x) = −960(x2−3)(3x2−1) (1 +x2)6 であるから、
maxn
|f(4)(ξ)|
0≤ξ≤1o
= maxn
|f(4)(ξ)|
ξ= 0, ξ= 1または(f(5)(ξ) = 0かつ0≤ξ≤1)o
= max
|f(4)(0)|, |f(4)(1)|, f(4)( 1
√3)
= max
|96|, | −12|, −81
2
= 96
となる。 それゆえ、
Z 1 0
4
1 +x2dxの近似値を求めるのにシンプソンの公式 を使った場合、小区間の個数がn= 10では
|誤差| ≤ (1−0)5
2880×104×96 +|計算に伴う誤差|
≈ 3.3×10−6+|計算に伴う誤差| となり、また n= 1000では
|誤差| ≤ (1−0)5
2880×10004×96 +|計算に伴う誤差|
≈ 3.3×10−14+|計算に伴う誤差| となる。
• 真値は
Z 1 0
4
1 +x2dx = π = 3.141592653589793238462643383279· · ·
(プログラミング) 与えられた式に従って計算するだけである。累算する箇所が2つあ
るが、そのうち4( f(a1) +f(a3) +· · ·+f(a2n−3) +f(a2n−1) ) の累算値を保持するため に sum4 という名前の double 型変数を、2( f(a2) +f(a4) +· · ·+f(a2n−2) ) の累算値 を保持するために sum2 という名前の double 型変数を用意してプログラムを構成した。
n = 1000としてこの計算を行うCプログラムと、これをコンパイル/実行している様子を
次に示す。(下線部はキーボードからの入力を表す。)
[motoki@x205a]$ nl numerical-integral-by-simpson.c Enter
1 /* 定積分 \int_0^1 \frac{4}{1+x^2} dx の値を Simpsonの公式 */
2 /* により数値的に求めて出力するCプログラム */
3 #include <stdio.h>
4 #define N (1000) 5 #define A (0.0) 6 #define B (1.0)
7 #define f(x) (4.0/(1.0+(x)*(x))) 8 int main(void)
9 {
10 double ans, sum4, sum2;
11 int i;
12 sum4 = 0.0;
13 for (i=2*N-1; i>=1; i-=2)
14 sum4 += f(A+(B-A)*((double)i)/(2.0*(double)N));
15 sum2 = 0.0;
16 for (i=2*N-2; i>=2; i-=2)
17 sum2 += f(A+(B-A)*((double)i)/(2.0*(double)N));
18 ans = (f(A)+f(B)+4.0*sum4+2.0*sum2)*(B-A)/(6.0*(double)N);
19 printf("定積分値(近似)は %22.16g\n", ans);
20 return 0;
21 }
[motoki@x205a]$ gcc numerical-integral-by-simpson.c Enter [motoki@x205a]$ ./a.out Enter
定積分値(近似)は 3.141592653589791
[motoki@x205a]$ տ下線部は真値と一致する箇所
ここで、
• プログラムの7行目 は引数付きマクロを定義したプリプロセッサ指令である。これに よって、以降(14行目,17〜18行目)に現れるf(x) という形の文字列は全てコンパイ
ル前に (4.0/(1.0+(x)*(x))) という文字パターンに置き換えられることになる。
• プログラムの13行目, 16行目 は、それぞれ
for (i=1; i<=2*N-1; i+=2), for (i=2; i<=2*N-2; i+=2) と書いた方が分かり易いかも知れない。 しかし、こう書くと次の繰り返しを進めるか どうかの判定の度に 2*N-1 もしくは 2*N-2 の計算を行わなければならないので、
与えられた式とは逆の順番に累算することにした。
□演習 8.4 (Simpsonの公式) Simpsonの公式を使えば、n= 10,100 の時に実際にどの
8.2. 数値積分 —Simpsonの公式— 125
程度の精度で Z 1
0
4
1 +x2dx の値が求まるか調べよ。
□演習 8.5 次の定積分の値を数値的に求めて出力するCプログラムを作成せよ。 [いず れの定積分値も π になることが、解析学の演習書の中に載っている。]
(1) Z 2
0
√4−x2dx
(2) Z 0.5
0
(12√
1−x2−√ 27)dx
□演習 8.6 (台形公式) Newton-Cotesの積分公式の中で、最も単純な(,そして最も近似 精度の悪い)ものは台形公式と呼ばれている。 台形公式によれば、一般に、定積分値は
Z b a
f(x)dx ≈ h 2 h
f(a0) +f(an)
+2( f(a1) +f(a2) +· · ·+f(an−2) +f(an−1) ) i 但し、h= b−na
ai =a+ih
と近似でき、この近似の際の誤差は
|誤差| ≤ (b−a)5
12n3 maxn
|f(2)(ξ)|
a≤ξ≤bo
と見積もられる。 この台形公式を用いて Z 1
0
4
1 +x2dx の値を数値的に求めて出力する Cプログラムを作成せよ。
9 配列
• 1次元配列を用いた計算,
• 整列化,
• 2次元配列,
• 付録 配列のまとめ
複数のものに添字番号付きの名前を付けると、それらを系統的に表せることがある。例 えば数学で、一般的な数列をa1, a2, a3, . . .という風に表したりする。この節では、この様 な「添字番号付きの名前」を使って大量のデータ領域のそれぞれに系統的な名前を付け、
それらの領域を規則的に扱う方法を紹介する。