数値積分
数値積分
•数値積分のプログラムを作る
–但し，
積分したい関数 f(x) = e
積分区間: [a, b]
区間数: n
-x2
定積分
区間 [a, b] で，連続関数 f, x軸,
x=a, x=b で囲まれた面積
f(x)
abx
b
adxxf )(
定積分：
区間 [a, b] の小区間への分割
•n 個の等間隔な小区間に分割
–幅： h = (b-a) / n
–小区間： [x0, x1], [x1, x2], ..., [xn-1, xn]
但し，x0= a, xi= x0+ i ×h
f(x)
x0= a xn= b
x
x1x2xn-1
小長方形
•小長方形の面積は
f(x)
x0= a xn= b
x
xixi+1
小長方形
f(xi)
hxf i)(
小台形
•小台形の面積は
f(x)
x0= a xn= b
x
xixi+1
小台形
f(xi)f(xi+1
)
))()((
21
ii xfxf
h
台形則 trapezoidal rule
•小台形の面積の和は
•定積分 I を，この和 Snで近似 ＝ 台形則
f(x)
x0= a xn= b
x
))()((
21
1
0
ii
n
i
nxfxf
h
S
x1x2xn-1
x3xn-2
台形則による数値積分
•区間［a，b］を n 等分 (1区間の幅h=(b-a)/n)
•n 個の台形を考え, その面積の和 Snで，定積分 I
を近似
–f(x) が連続関数のときは，n を無限大に近づければ，
SnはI に近づく
•但し，単純に「n を大きくすればよい」とは言えない
–n を大きくすると ⇒計算時間の問題，丸め誤差の問
題
n
n
ii
n
i
nSxfxf
h
I
lim))()((
2
lim 1
1
0
数値積分の意味
•式で書いてある関数の積分
⇒ごく限られた関数しか積分できない
•定積分の「数値積分」が重要になる
数値積分の精度
•分割幅を小さくするほど高精度
•近似式の次数が高いほど高精度
•台形則より次数が１高い近似の手法
Simpson則
区間を２次曲線により近似
台形則：区間を直線（１次式）により近似
Simpson則
• 積分値∫ f(x)dx を計算(関数f(x)をある区間
[a,b] で積分)
– 積分区間を単位(幅2h)に分割
– 各単位を二次曲線で近似
• 座標 (-h,f(-h)), (0,f(0)), (h,f(h)) の3点を通る二次曲線
–区間 -h～h までの積分を計算
– それらの合計により全体の積分値を求める
a
b
0-h h
f(0) f(h)f(-h)
x
f(x)
Simpsonの公式
•3点を通る二次曲線 f(x)=Ax2+Bx+C の区間
[-h,h] での積分
∫f(x) dx =∫ (Ax2+Bx+C)dx
= h(2Ah2+6C)/3
= h[(Ah2-Bh+C)+4C+(Ah2+Bh+C)] / 3
= h(f(-h)+4f(0)+f(h)) / 3
h
-h h
-h
Simpsonの公式
積分値∫ f(x)dx の計算手順
•x=a から x=b までを n 等分し，その幅を2h(=(b-a)/n)と
する
∫f(x)dx = h[f(a)+4f(a+h)+f(a+2h)
+f(a+2h)+4f(a+3h)+f(a+4h)
+ ・・・
+f(a+(n-2)h)+4f(a+(n-1)h)+f(a+nh)]/3
=h[f(a)+f(a+nh)
+4(f(a+h)+f(a+3h)+・・・+f(a+(n-1)h))
+2(f(a+2h)+f(a+4h)+・・・+f(a+(n-2)h))]/3
これで積分値が求まる
a
b
a
b
Simpson公式のプログラム
• 積分区間と分割数を読み込んで，Simpson
公式による数値積分を行う
–関数f(x) はプログラム中に記述
#include <stdio.h>
#include <math.h>
double f(double x)
{return exp(-x*x);
}
main()
{int i, num;
double S, a, b, x, d;
char buf[100];
printf("積分区間(a～b) : ");
fgets(buf, 100, stdin);
sscanf(buf, "%lf %lf", &a, &b);
printf("分割数(偶数) : ");
fgets(buf,100,stdin);
sscanf(buf, "%d", &num);
d = (b - a) / num; /* 分割幅 */
S = f(a) + f(b);
for(i = 1 ; i < num ; i = i + 2){
x = a + d * i;
S = S + 4 * f(x);
}
for(i = 2 ; i < num ; i = i + 2){
x = a + d * i;
S = S + 2 * f(x);
}
S = S * d / 3; /* 面積 */
printf("面積 = %lf\n", S);
}
実行結果
f(x) = exp(-x2)
積分区間(a～b) : 0 1
分割数 : 1000
面積 = 0.746824
積分区間(a～b) : 0 10
分割数 : 1000000
面積 = 0.886277
練習問題１
•台形則を使って，次を計算せよ
•計算結果を，以下と比較せよ
1
01
2log dx
x
x
e
log 2 は，およそ 0.6931471805599453
練習問題２
•練習問題１について，区間数 n を，n = 4, 8,
16, 32, 64, 128 と変えて計算を行い，刻み
幅と誤差の関係を調べよ （A）
–区間数 n と誤差の関係を書いたグラフを作成
せよ
–この場合，おおよそ次の関係が成り立ってい
ることを確認せよ
2
n
c
ISn
（c は定数）
