前回資料の訂正

2

スプライン補間のアルゴリズム
入力: (x0, y0), (x1, y1), … , (xN, yN),x0< x1 < ∙∙∙ < xN

出力: Sj(x) = aj(x - xj)3 + bj (x - xj)2 + cj (x - xj) + dj
(j = 0, … , N - 1) （もしくは (aj , bj , cj , dj ) の組）
1. u0 ← 0; uN ← 0; h0 ← x1 - x0;
2. for j ∈ [1..N-1] do

3

計算機数学II (2018)
4: 積分
照井 章（筑波大学 数理物質系 数学域）
Akira Terui (Institute of Mathematics, University of Tsukuba)4第4章の内容● 数値積分
○ 台形則○ シンプソン則とロンバーグ積分法

54-1 微分と積分

6

積分の数値計算の必要性

微分原始関数

導関数
不定積分

7

初等関数の不定積分は常に初等関数になるとは限らない● 初等関数 … 次の関数、もしくはそれらの有限回の合成によって定まる関数○ 代数関数（多項式）○ 指数関数 (exp(x))○ 対数関数 (log(x))○ 三角関数 (sin(x), cos(x), … )○ 逆三角関数 (Arcsin(x), Arccos(x), … )8初等関数の不定積分は常に初等関数になるとは限らない
例1:
導関数:

(4.1)

不定積分は初等関数で表現できない(4.2)
Erf(x): 誤差関数
9

初等関数の不定積分は常に初等関数になるとは限らない
例2: f(x) = (1 - k2 sin2 x)-1/2
導関数:

(4.3)

不定積分: 第1種楕円積分

(4.4)

10

定積分の数値計算=数値積分● 初等関数の導関数は初等関数● その逆は必ずしも成り立つとは限らない● しかし（不定）積分を計算する必要性が生ずる● そんなとき、定積分を求めることができれば、不定積分（原始関数）の計算の役に立つ

11定積分の数値計算=数値積分● 数値積分（法） (numerical integration):（定）積分の値を数値計算で求める● 数値積分の条件
a. 被積分関数 f(x) が与えられている
b. （積分区間内の）x に対し、f(x) の値を求めることができる
12

積分の問題は常微分方程式の問題に帰着させることが可能● 計算目的:● 1階常微分方程式(4.6)
を y (a) = 0 の条件下で解き、y(b) の値を求める● 式 (4.6) の両辺をaからbまで積分すると
(4.7)
13

積分の問題は常微分方程式の問題に帰着させることが可能● さらに条件 y(a) = 0 を加えると(4.8)● 常微分方程式(4.6) は(4.9)に含まれる
14しかし、積分の問題は常微分方程式の問題とは分けて論じる●は右辺に y を含んでいるため、任意の x に対しては直ちには計算できない● そこで、数値積分は常微分方程式と区別して論じる

15

4-2 台形則

16台形則の原理(4.10)
y = f(x)+ab-17

台形則の原理
積分領域を N 個の台形で近似 (0, … , N-1)
各台形の横幅 h は等しい: h = (b - a) / N
f(a+jh)b-haa+hba+jh18

台形則の原理
積分領域を N 個の台形で近似 (0, … , N-1)
左から j 番目の台形の上底と下底の長さ:それぞれ f(a + j h), f(a + (j + 1) h)f(a+jh)b-haa+hba+jh19

台形則の原理
積分領域を N 個の台形で近似 (0, … , N-1)
左から j 番目の台形の上底と下底の長さ:それぞれ f(a + j h), f(a + (j + 1) h)
このとき、全部の台形の面積の和 T は(4.11)20台形則の原理● N = 1 の場合：T = h {f(a)+f(b)}/2● 台形則 (trapezoidal rule)● 各区間の左側の x 座標 xj = a + j h :分点 (abscissa)21

台形則の計算例● (a)(b)● 教科書: 表 4-1● 積分区間の分割数 h を2倍ずつ増やし
誤差 |T - I| を調べる
👉 分割を細かくすれば T が I の値に近づく
2● 誤差 ∝ 1/N (4.12)22台形則の誤差と分割数の関係● 区間 [xj, xj+1] において、関数 f(x) を折れ線関数
F(x) で近似していると見ることができる● F(x) は 2点 (xj, f(xj)) と (xj+1, f(xj+1)) を結ぶ線分
👉 1次のLagrange補間に他ならない(4.13)
23

台形則の誤差と分割数の関係
f(x) と F(x) の誤差: 式 (3.11) より(4.14)
よって，xj < x < xj+1 においてと評価できる(4.15)
24

台形則の誤差と分割数の関係[xj , xj+1] における積分値と台形の面積の誤差は(4.16)（ただし h = xj+1 - xj）
25

台形則の誤差と分割数の関係[a, b] における積分値と台形の面積の誤差はNを大きくすると1/N2に比例して小さくなる(4.17)
26

台形則の計算終了条件と計算手順● 台形則を用いて積分値をp桁の精度で求めるとする(4.17)●の右辺の max |f’’(ξ)| を見積もらなければならないが、これは一般に煩雑● そこで、分割数 N を増加させ、先頭から p 桁の値が変化しなくなるまで計算を続ける
27

台形則の計算終了条件と計算手順● 台形則の計算のほとんどは f(x) の計算● 分割数Nを2のべき乗にし、前のステップで求めた分点における f(x) の値をそれ以降の新たなステップでも利用する

28

台形則の計算終了条件と計算手順

Tn: N = 2n における台形則の値
このとき、台形の幅 h は (b-a)/2n29台形則の計算終了条件と計算手順(4.18)
より(4.19)30

台形則の計算終了条件と計算手順
ただし T0 = (b-a){f(a) + f(b)} / 2● 教科書 表 4-1 (a), (b) どちらの場合も精度6桁で積分値を求めるには N=1024, すなわち f(x) の計算が1025回必要

31

台形則のアルゴリズム

入力: f(x): 被積分関数, a < b: 積分区間,p: 精度（10進桁数）
出力:の近似値（精度 p 桁）

1. ε ← 10-p; N ← 1; h ← b - a; T ← h{f(a) + f(b)}/2;

32

台形則のアルゴリズム
2. while (1)a. N ← 2N; h ← h/2; s ← 0;
b. for i = 1, 3, … , N-1 doi. s ← s + f(a + i h);c. new_T ← T/2 + h∙s;d. if (| new_T - T | < ε | new_T |) then break;e. T ← new_T
3. return new_T;33

4-3 シンプソン則とロンバーグ積分法

34シンプソン則の原理● （複合）シンプソン則 (Simpson’s rule)● 区間の分割数 N に対し、誤差が 1/N4 に比例● 2次のLagrange補間多項式を利用

35

シンプソン則の原理● 被積分関数: f(x)● 積分区間: [a, b] (a < b)
● 分割数: N● 隣り合う分点の間隔: h = (b - a) / N● 分点: xj = a + j h (j = 0, 1, … , N) （N は偶数）36

シンプソン則の原理● 隣り合う3つの分点 xj, xj+1, xj+2 を含む区間[xj, xj+2] において、f(x) を近似する2次のLagrangeの補間多項式 P(x) を作る● P(xj) = f(xj), P(xj+1) = f(xj+1), P(xj+2) = f(xj+2)37

2次のLagrange補間多項式による補間の例
P(xj) = f(xj), P(xj+1) = f(xj+1), P(xj+2) = f(xj+2)

y = f(x)y = P(x)xj

xj+1xj+238

シンプソン則の原理(4.20)

39

シンプソン則の原理
区間 [xj, xj+2] での f(x) の積分を P(x) の積分で近似する(4.21)
ただし x = xj + t h, xj = a + j h40

シンプソン則の原理
区間 [a, b] における積分:● [a, b] を [x0, x2], [x2, x4], … , [xN-2, xN] の小区間に分ける（N は偶数とする）● 各小区間で式 (4.21) に基づいて得られる定積分の近似値を合計する● 合計値: S41

シンプソン則の原理シンプソン則 (Simpson’s rule)(4.22)

42

シンプソン則の原理の中括弧内において● f(x0), f(xN) の係数は1● f(x1)からf(xN-1)までの係数は 4, 2, 4, 2, … , 4 のように4と2を繰り返す● N = 2 の場合: S = (h/3){f(x0) + 4 f(x1) + f(x2)}43シンプソン則の計算例● (a)(b)● 教科書: 表 4-2● (a) N = 16, (b) N = 32 の際にそれぞれ精度6桁の答えを得られる（台形則では (a), (b) とも N = 1024 のとき）
● 例 (a), (b) においては誤差 ∝ 1/N

444

シンプソン則の誤差評価● 式 (4.22) の S と真の積分値の誤差は(4.23)● ただし、台形公式の誤差評価式 (4.17) との大小関係は一概には比べられない● あるNに対してシンプソン則の方の精度が悪くても、 Nを大きくすれば台形公式より速く精度が向上する可能性がある
45

シンプソン則と台形則の関係● シンプソン則を用いて積分値を p 桁の精度で求めるとする● 式 (4.23) の右辺の max |f(4)(ξ)| を見積もらなければならないが、これは一般に煩雑● そこで、台形則同様、分割数 N = 2n (n = 1, 2, 3, … )を増加させ、先頭から p 桁の値が変化しなくなるまで計算を続ける46

シンプソン則と台形則の関係● 台形則で得られた答えを使ってシンプソン則による積分値を計算する● 積分区間 [a, b] の分割数: N = 2n とする○ Tn: 台形則による計算結果○ Sn: シンプソン則による計算結果● このとき(4.24)
47

シンプソン則と台形則の関係【証明】 h = (b - a)/(2n+1) とすると次式が成り立つ:(4.25)48

シンプソン則のアルゴリズム
入力: f(x): 被積分関数, a < b: 積分区間,p: 精度（桁数）
出力:の近似値（精度 p 桁）

1. ε ← 10-p; N ← 2; h ← (b - a)/2;
T ← h{f(a) + 2f((a +b)/2) + f(b)}/2;
S ← h{f(a) + 4f((a +b)/2) + f(b)}/3;
49

シンプソン則のアルゴリズム
2. while (1)a. N ← 2N; h ← h/2; s ← 0;
b. for i = 1, 3, … , N-1 doi. s ← s + f(a + i h);c. new_T ← T/2 + h∙s; // 台形則の更新
d. new_S ← (4∙new_T - T)/3; // 式 (4.24) の更新
e. if (| new_S - S | < ε | new_S |) then break;f. T ← new_T; S ← new_S;
3. return new_S;50

ロンバーグ積分法● 台形則からシンプソン則を導いた(4.24)
をより一般化し、繰り返し計算してより速く誤差が小さくなるような計算法

51

ロンバーグ積分法● 分割数: N = 2n (n = 0, 1, 2, … )● Tn(0): 台形則によって計算された積分値● Tn(k) を求める漸化式 (n ≧ k, k = 1, 2, 3, …):(4.26)
k=1のとき(4.27)（シンプソン則による積分近似値）

52

ロンバーグ積分法
Tn(k)を計算する手順（図 4-7）

53

ロンバーグ積分法の計算例● (a)● 教科書: 表 4-3● T4(0) の計算における分点数 (= 被積分関数の評価回数) = 17
● T4(4)の計算値における誤差 ≦ 1.0 × 10-13● 少ない関数計算でよい精度の結果を得られる
54

第4章のまとめ● 数値積分
○ 台形則○ シンプソン則とロンバーグ積分法

55

次回の内容● 第5章: 常微分方程式○ 5-1: 常微分方程式○ 5-2: 微分と差分○ 5-3: 差分方程式

56