数値計算入門
武尾英哉1. 離散数学と数値計算
数学的解法の中には理論計算では求められないものもある．例えば，定積分は，まず は積分（被積分関数の原始関数をみつけること）できなければ値を得ることはできない． また，ある関数の所定の値における微分値を得るには，まずその関数の微分ができなけ ればならない．さらに代数方程式の解を得るためには，解析的に代数方程式を解く必要 がある．ところが，これらは必ずしも解析的に導けるとは限らない．そこで本章では， 計算機を使って，離散的な数値（とびとびの数値を扱う）を扱って，近似的に解を導く 手法について説明する．このことを離散数学という． 本稿では，連続的な数値による数学（これまでの解析的数学）を，離散的な数値によ る数学として近似的に扱うことで解を求めることができる手法について説明する．この ように，離散数学を用いて種々の手法による近似計算により解を求めることを数値 計算とよぶ．2. 数値代数方程式（ニュートン法）
数値計算を使って代数方程式を解く代表的な手法であるニュートン法について説明す る． 代数方程式を解く（f x
( )
0
を満足するx
を求める）ということは，グラフにおいて 関数y
f x
( )
のx
軸切片を求めることと等価である．すなわち，図 1 におけるx
a
を 求めることである．以下，ニュートン法による手順について説明する． ①まず初めに，予想される真の解に近いと思われる値をひとつとる（x
の初期値）．x
の初期値をx
0とする． ②次に，グラフの接線を考え，そのx
軸切片を計算する．図 1 において，x
x
0にお ける接線の方程式はy
f x
'( )(
0x x
0)
f x
( )
0 であるので，この接線がx
軸と交差す る点は，x
が 0 0 0(
)
'(
)
f x
x
f
x
のときである．このx
を新たにx
1とする． 0 1 0 0(
)
'(
)
f x
x
x
f
x
ただし，f x
'( )
0
0
とする． このx
軸切片の値は，予想される真の解により近いものとなるのが一般である． ③②の手法を繰り返すと，限りなくx
a
に近い値となる．すなわち，一般的にx
x
nから次のステップの
x
x
n1を求める式は次のようになる． 1(
)
'(
)
n n n nf x
x
x
f
x
ただし，f x
'( )
n
0
とする． (1) ④この繰り返し計算を希望する精度の解を得るまで行う．すなわち，実際には，次の 収束判定条件を計算し，満足するまで繰り返すことで，解の近似解を得る． 1 1 1 n n nx
x
x
またはx
n1
x
n
2 (2) ここで，
1，
2は収束条件である（例えば，0.01 など微小な値を設定する）．なお， 収束条件の詳細に関しては，数値計算や数値解析等の関連書籍を参照されたい． ニュートン法は，一般的に早く収束するが，初期値x
0のとり方によっては収束しない こともある．また，解が複数ある場合には，初期値の設定により別の解に収束する場 合がある．よって，予めおおまかな解を初期値としてうまく設定する必要がある．y
y
f x
'( )(
0x x
0)
f x
( )
0y
f x
( )
x
図 1: ニュートン法の原理 [例1]2
を計算で求める場合として，f x
( )
x
2
2
0
の解をニュートン法で求める．'( )
2
f
x
x
であるので，式(1)は 2 12
2
n n n nx
x
x
x
と書き表せる． ここで，初期値を例えばx
0
1
とおくと，x
n1は2
に収束し，一方，x
0
1
と おくと，x
n1は
2
に収束していく．a
0x
1x
3. 数値微分（差分法）
数値計算を使って関数の所定値の微分値を求める代表的な手法である差分法について 説明する．数値計算では無限小を扱うことができないので，小さいが有限の値を用いた 微分演算の近似である差分を用いる．ここで，微分はその点での傾き 0(
)
( )
'( )
lim
hf x
h
f x
f
x
h
(3) で定義されるが，差分では有限距離離れた 2 点の傾き(
)
( )
f x
h
f x
h
(4) である．差分は下記のように，前進差分，後退差分，中央差分の 3 つがある（図 2 参 照）．( )
y
f x
のx
x
0における微分係数は，微小区間h
について近似的に下記のように求 めることができる． （１）前進差分 3 2 0 0 0(
)
(
)
'(
)
f
f
f x
h
f x
f
x
h
h
(5) （２）後退差分 0 0 2 1 0(
)
(
)
'(
)
f
f
f x
f x
h
f
x
h
h
(6) （３）中央差分 3 1 0 0 0(
)
(
)
'(
)
2
f
f
f x
h
f x
h
f
x
h
h
(7) 0x
h
x
0x
0
h
図 2: 数値微分の原理y
( )
y
f x
1f
2f
3f
x
4. 数値積分（台形法とシンプソン法）
定積分の数値計算法としては，x
をn
分割し各区間の両端のf x
( )
を直線で結んで台形 の面積の総和を求める方法（台形法）と，2 区間を 2 次式で近似して面積の総和を求め る方法(シンプソン法)などがある．解析的に積分が不可能な場合の定積分を求めるとき に有効であり，数値積分とよばれる． （１）台形法 図 3 に示すように，区間[a,b]をn
等分して，x
x x
k,
k1における関数f x
( )
の値 1( ), (
k k)
f x
f x
を利用した微小台形を作成する．この台形の面積は，下式で計算され， 微小領域の面積を近似する． 1 1 1(
)
(
)
(
)
(
)
(
)
2
2
k k k k k kf x
f x
b
a
f x
f x
x
x
n
(8) ここで，x
k1
x
kは台形の高さ，f x
( ), (
kf x
k1)
は台形の上底と下底を表している． 一方，微小領域の面積
S
は定積分を用いて下記のように表せるので，式(8)と併せて 次式が成り立つ． 1(
)
(
1)
( )
2
k k x k k xf x
f x
b
a
S
f x dx
n
(9) よって，区間[a,b]の面積S
は，
1 0 1 1 2 2 1 1 0(
)
( )
( )
(
)
(
)
(
)
(
)
(
)
2
n n n n kb a
S
S
f x
f x
f x
f x
f x
f x
f x
f x
n
(
0)
2
( )
1(
2)
(
1)
(
)
2
n nb
a
f x
f x
f x
f x
f x
n
(10) となり，定積分 b( )
af x dx
は次式で近似できる．これを台形法とよぶ．
0 1 2 1
( )
(
)
2
( )
(
)
(
)
(
)
2
b n n ab
a
f x dx
f x
f x
f x
f x
f x
n
(11)図 3: 台形法の原理 図 4: シンプソン法の原理 （２）シンプソン法 図 4 に示すように区間[a,b]を
2n
等分して，微小区間x
2k
x
x
2k2で関数f x
( )
を
x
2k, (
f x
2k) ,
x
2k1, (
f x
2k1) ,
x
2k2, (
f x
2k2)
の3 点を通る 2 次曲線で補間するこ とを考える． 2 次曲線を 2( )
g t
p qt
rt
と定義し，微小区間[
x
2k,
x
2k2]
における面積
S
は下 記のように表せる．
2 2 2 21
21
3( )
( )
2
3
k k h x h h x h h hS
f x dx
g t dt
p
qt
rt dt
pt
qt
rt
2 3 2 3 3 21
1
1
1
2
2
6
2
2
3
2
3
3
3
h
ph
qh
rh
ph
qh
rh
ph
rh
p
rh
(12) ここで， 2 2(
k)
(
)
f x
g
h
p
qh
rh
2 1(
k)
(0)
f x
g
p
2 2 2(
k)
( )
f x
g h
p
qh
rh
から， 2 2 2 1 2 26
p
2
rh
f x
(
k)
4 (
f x
k)
f x
(
k)
(13) 式(13)を式(12)に代入して，
2 2 1 2 2
( )
(
)
4 (
)
(
)
3
h k k k hh
g t dt
f x
f x
f x
x
y
( )
y
f x
0x
a
x
kx
2k1x
2k2y
x
( ) y f x 1 kx
x
n
b
( k) f x 2kx
1 ( k ) f x 0x
a
x
2n
b
2n
等分n
等分2
b a
h
n
であるので， 2 2 2 2 2 1 2 2(
)
4 (
)
(
)
( )
( )
2
3
k k x h k k k x hf x
f x
f x
b a
S
f x dx
g t dt
n
よって，
1 1 2 2 1 2 2 0 0( )
(
)
4 (
)
(
)
6
n n b k k k a k kb a
S
f x dx
S
f x
f x
f x
n
(
0)
4 ( )
1(
2)
(
2)
4 ( )
3(
4)
6
b
a
f x
f x
f x
f x
f x
f x
n
0 1 3 2 1{ ( ) 4
( )
( )
(
)
6
nb a
f x
f x
f x
f x
n
2 4 2 2
22
f x
( )
f x
( )
f x
(
n)
f x
(
n)}
(14) となる．これをシンプソン法とよぶ． ただし，シンプソン法では区間を2n
等分するので，偶数分割することが前提となる．5. 数値微分方程式（オイラー法）
本節では，微分方程式の数値解析の一つであるオイラー法について説明する． 次式の微分方程式を数値解析によって解くこととする．( , )
dy
f x y
dx
ただし，初期値x
0のときy
0とする（y x
( )
0
y
0） (15) 微分の定義から，(
)
( )
lim
hdy
y x
h
y x
dx
h
(16) これがf x y
( , )
と等しい． 一方，h
が十分に小さいとき，式(16)の右辺はy x
(
h
)
y x
( )
h
に近似することから，(
)
( )
( , )
y x
h
y x
f x y
h
(17) すなわち，(
)
( )
( , )
y x
h
y x
hf x y
(18) が成り立つ．初期値を
x y
0,
0
y x
( )
0 とすれば， 1 0,
1( )
1 0( ,
0 0),
2 1,
2( )
2 1( ,
1 1),
x
x
h y
y x
y
hf x y
x
x
h y
y x
y
hf x y
と順 次にy
i
y x
( )
i を計算できる．微分方程式の数値解析とは，これらの点( ,
x y
i i)
の集まり を求めることであり，さらに平面上に( ,
x y
i i)
をプロットすると求める関数y
のグラフが 描ける．このようにして，微分方程式の解（関数y
）がグラフ上に得られることになる． オイラー法をまとまると，以下の手順になる． 初期値を与え，漸化式y
k1
y x
( )
k
hf x y
( ,
k k)
を順に計算すればよいので， ①初期値x y
0,
0
y x
( )
0 を決める． ②y
1
y x
( )
1
y x
(
0
h
)
y x
( )
0
hf x y
( ,
0 0)
y
0
hf x y
( ,
0 0)
を計算する． ③y
2
y x
( )
2
y x
(
1
h
)
y x
( )
1
hf x y
( ,
1 1)
y
1hf x y
( ,
1 1)
を計算する． ④y
3
y x
( )
3
y x
(
2
h
)
y x
( )
2
hf x y
( ,
2 2)
y
2
hf x y
( ,
2 2)
を計算する． ⑤順次，y
k1
y x
(
k1)
y x
(
k
h
)
y x
( )
k
hf x y
( ,
k k)
y
k
hf x y
( ,
k k)
を計算す る． [例2] 微分方程式dy
y
dx
，初期値をx
0
0,
y
0
1
，h
0.5
とするとき，オイラー法に より微分方程式の解を求める．( , )
f x y
y
，前記オイラー法の手順より， 1(0 0.5)
(0.5)
0( ,
0 0)
00.5
( ,
0 0) 1 0.5 1 1.5
y
y
y
y
hf x y
y
f x y
2(0.5 0.5)
(1)
1( ,
1 1)
10.5
( ,
1 1) 1.5 0.5 1.5
2.25
y
y
y
y
hf x y
y
f x y
3(1 0.5)
(1.5)
2( ,
2 2)
20.5
( ,
2 2)
2.25 0.5 2.25 3.375
y
y
y
y
hf x y
y
f x y
となり，(0)
1, (0.5)
1.5, (1)
2.25, (1.5)
3.375,
y
y
y
y
となる． 一方，与式の微分方程式の解析解は，y x
( )
e
xより，(0)
1, (0.5)
1.649, (1)
2.718, (1.5)
4.482,
y
y
y
y
となり，誤差は図5のように変化する．また，刻み幅
![図 3: 台形法の原理 図 4: シンプソン法の原理 （２）シンプソン法 図 4 に示すように区間[a,b]を 2n 等分して，微小区間 x 2 k x x 2 k 2 で関数 f x ( ) を x 2 k , ( f x 2 k ) , x 2 k 1 , ( f x 2 k 1 ) , x 2 k 2 , ( f x 2 k 2 ) の 3 点を通る 2 次曲線で補間するこ とを考える． 2 次曲線を g t ( ) p](https://thumb-ap.123doks.com/thumbv2/123deta/8515189.925858/5.892.185.715.153.438/シンプソンシンプソンx次曲線考える次曲線.webp)
