6

第1章

1.2

数値と誤差

1.2.1

数値の表現

数値の誤差と計算

コンピュータ内部では，数値は 10 進法ではなく 2 進法で表現されている．
その理由は，2 進法を用いると計算を行なう電子回路が単純になって処理が速くなること，数値の記憶方法が効率的になること，などの利点があるからである．一方，我々が普段扱う数は 10 進法で表現されているため，2 進法との変換が必要になる．ところが，この変換により，やっかいなことが生じる．
例えば，10 進法の 0.1 を 2 進法に変換すると 0.000110011001100 …という無限小数となる．しかしながら，コンピュータは有限の桁しか扱うことができないため，10 進数の 0.1 をコンピュータで扱うと必然的に誤差を生じる．
それゆえ , 簡単な四則演算であっても誤差が混入する可能性がある．電卓で1 ÷ 3 × 3 を計算すると，計算結果は 1 ではなくて 0.9999999 …と表示されることがあるが，同様のことはコンピュータでも起こっている．本書では簡単のため 10 進法を用いて説明するが，コンピュータ内部では 2 進法が用いられていることに注意してほしい．
コンピュータが扱う数値には，主として整数型と実数型の 2 種類があって，それぞれ表現の仕方が異なっている．
整数型は文字通り整数を扱うための形式で，ある大きさ以下の整数を誤差なく表現できる．例えば，2 進数 32 桁（32 ビット）で符号付きの整数を表現するときは，−231 から 231 − 1 までの整数を誤差なく表現できる．
無限小数や，整数型では表せないような大きな整数を表現するためには，実数型と呼ばれる形式で数値を表現する．例えば，小数点以下何桁まで表示
するかを指定して，実用上は十分な精度で数を表現する．有効数字の桁数を重視する場合には，数の正負を表す符号部，有効数字を表す仮数部，そして位取りを表す指数部の組で表現する．これを浮動小数点表示という．例えば，円周率 π を
＋0.31415926 × 101

1.2

数値と誤差

7

などと表現したりする．ここで，＋ が符号部，0.31415926 が仮数部，101が指数部である．普通，仮数部としては 0 と 1 の間の有限小数をとり，小数点以下の有効桁数をなるべく多くとるために，小数第 1 位を非零の数になるようにする．また，101 が指数部で，10 の整数乗の形をとる．
なお，浮動小数点表示の方法は，IEEE1) が定めた規格2) があり，それにしたがって設計されているコンピュータが多い．
高い精度を確保するために，桁数が多い数値計算を必要とするときには，
普通より仮数部の桁数を大きくとった多倍長実数計算のためのソフトウェア
パッケージを利用することも可能である．ただし，この場合には，より多くの計算機資源を使うので注意が必要である．
数が大きすぎると，コンピュータ上で扱える数値の範囲を超えてしまうことがある．範囲を超えるような大きすぎる数が生じることをオーバーフローという．逆に，0 とほとんど区別できないような小さな数が生じると，十分な有効桁数を確保できないことがあり，これをアンダーフローという．いずれの場合にも，コンピュータで数値計算を行なうときの障害となるので，このようなことが起こらないように注意が必要である．

1.2.2 誤 差
以下では，x によって近似値を，x （t は true の頭文字）によって真の値を表すとする．このとき，e：＝ x − x 
を誤差といい， e  を絶対誤差，
e
e
≃
x
x
を相対誤差という．ここで=≃?は左辺が右辺で近似されることを表す記号である．
1)

Institute for Electrical and Electronics Engineers

2)

IEEE Standard for Floating Arithmetic 754 - 2008

8

第1章

数値の誤差と計算

誤差そのものはわからないが，その大きさがある程度推定できるとき，例えば e  ≤ δを満たす数 δ がとれるとき，δ のことを誤差限界という．
例1
円周率 x  ＝ π ＝ 3.1415926 …の近似値として x ＝ 3.14 をとると，誤差，絶対誤
差，相対誤差，誤差限界は，それぞれ
誤差：

e ＝ − 0.0015926 …

絶対誤差：  e  ＝ 0.0015926 …
相対誤差：

e
＝ − 0.0005072 … ≃ 0.05 ％
x

誤差限界：

δ ＝ 0.0016

となる．

□

1.3

誤差の発生

1.3.1

誤差の分類

計算の過程において，誤差はいろいろな形で混入してくる．ここでは，誤差をその発生原因によって分類してみよう．
丸め誤差

実数を有限桁数の小数で表現したときに生じる誤差のことである．仮数部の長さはコンピュータによって決まっており，それ以上の長さの有限小数や無限小数をそのまま表現することはできない．このような場合には，数値の上位何桁のみを取り出し，それより下の桁を切り捨て，切り上げ，四捨五入（2 進法の場合は 0 捨 1 入）などの操作を行なうことによって近似的に数値を表現する．このように，長い桁数の数値を短くすることを丸めるといい，この際に生じる誤差を丸め誤差という．
打ち切り誤差無限級数で表現されている数値に対し，級数を有限の部分和で近似することによって生じる誤差を打ち切り誤差という．例えば，

34

第2章

非線形方程式＝ (K ＋ ε＋1)(K ＋ ε)(x − α)2−{(K ＋ ε＋1)(K ＋ ε)(x − α) − (K ＋ ε)(x − α)}(K ＋ ε＋1)(K ＋ ε)(x − α) − 2(K ＋ ε)(x − α) ＋ (x − α)＝ (K ＋ ε＋1)(K ＋ ε)(x − α)2−{(K ＋ ε＋1)(K ＋ ε) − (K ＋ ε)}(x − α)(K ＋ ε＋1)(K ＋ ε) − 2(K ＋ ε) ＋ 1 が得られる．これより
x − αx − α2＝ (K ＋ ε＋1)(K ＋ ε) −→ K2 −(K 2 − K)2
K 2 − 2K ＋ 1{(K ＋ ε＋1)(K ＋ ε) − (K ＋ ε)}(K ＋ ε＋1)(K ＋ ε) − 2(K ＋ ε) ＋ 1(n → ∞)＝ K2 − K2 ＝ 0
を得る．
注意 2

□

エイトケンの加速法は適用範囲が広く，1 次収束していない数列にも適用できる場合がある．一方で，もとの数列が収束しない場合であっても，エイトケンの加速法で得られた数列が収束することがあるが，このときには必ずしも解に収束するとは限らないため，注意を要する．

2.5

ニュートン法

2.5.1

ニュートン法

代入法は一般には f(x) の関数値だけを用いた手法であるが，微分係数の値を用いることによって，より収束の速い強力な方法が得られる．
(x0) の値が計算できた
いま x0 を近似解とし，その点において f (x0) と f ′と仮定しよう．このとき，x - y 平面において，点 (x0，f (x0)) の近くでは
y ＝ f(x) のグラフはこの点における接線で近似できる．そこで接線が x 軸0) とすると，図 2.5 のように，x1 は x0 よりも解 α に近と交わった点を (x1，

2.5yニュートン法35

y ＝（x）faO

xn ＋ 2 xn ＋ 1

xn

x

図 2.5

ニュートン法の原理

いことが期待できる．接線の方程式は
y ＝ f′(x0)(x − x0) ＋ f (x0)であるから，y ＝ 0 ， x ＝ x1 とすると，x1 は
x1 ＝ x0 −f (x0)
f′(x0)と求められる．
f(x1)) において接線を引き，x 軸との交点を (x2，0) とす
同様に，点 (x1，れば，x2 は x1 よりもさらに α に近くなる．この操作を繰り返し，数列 x1，x2，x3，… をx＋1 ＝ x −f (x)
f′(x)(n ＝ 0，
1，2，
…)

(2.5)

で計算すれば，この数列の極限として解 α が計算できることが期待される．
このようにして，解の逐次近似列を構成する方法をニュートン法という．
ニュートン法x＋1 ＝ x −f (x)
f′(x)(n ＝ 0，1，2，…)

例題 2.4方程式 x 5 ＋ x − 1 ＝ 0 の解をニュートン法で求めよ．

【解】 f(x) ＝ x 5 ＋ x − 1 とおくと，ニュートン法の公式は

36

第2章

非線形方程式x＋1 ＝ x −x5 ＋ x − 1
5x4 ＋ 1となる．f(0) ＝ − 1 < 0，f(1) ＝ 1 > 0 であるから，区間 (0，1) 内に少なくとも1 つの解がある．そこで，初期値として x0 ＝ 0.5 をとると，
x1 ＝ 0.5 −(0.5)5 ＋ 0.5 − 1＝ 0.85714286 …
5・(0.5)4 ＋ 1と計算でき，同様の計算を繰り返すと，表 2.2 のような結果が得られる．
表 2.2 ニュートン法による非線形方程式の解の計算
nx

f (x，y)

0
1
2
3
4

0.500000
0.8571430.770682
0.755283
0.754878−0.468750
0.319807
0.042562
0.001064
0.000001

これより，f (x) の値が急速に 0 に近づき，近似解は 0.754878 … に収束していくことが見てとれる．

2.5.2

□

ニュートン法の収束

ニュートン法は，代入法 x＋1 ＝ g(x) において g(x) ＝ x − f(x)  f ′(x)とした場合に相当する．実際，x ＝ α を方程式 f(x) ＝ 0 の解で f ′(α) ≠ 0を満たすものとするとg(α) ＝ α −

f (α)
＝α
f′
(α)

が成り立つ．この変換関数 g を用いて，解に収束する逐次近似列が構成されるための条件について調べてみよう．
定理 2.4 f(x) を R 上の C 2 級の関数とし，ある α に対して f(α) ＝ 0，f′(α) ≠ 0 を満たすとする．このとき，ある正数 δ > 0 がとれて，区間[α − δ，α ＋ δ] 内の任意の初期値 x0 に対し，ニュートン法によって得られる数列 { x } は方程式 f(x) ＝ 0 の解 α に一意に収束する．

50

3.1
3.1.1

第3章

代数方程式

代数方程式の性質
代数方程式

f(x) が x の n 次多項式のとき，方程式 f(x) ＝ 0 を n 次代数方程式といい，f(x) が多項式ではない非線形関数のとき，方程式 f(x) ＝ 0 を超越方程式という．代数方程式には超越方程式にはない特別な性質がいくつかあり，異なる取り扱いが可能となる．代数方程式においては，最高次の係数が 1 であるとしても一般性を失わないことに注意し，この章では主に
x  ＋ a1 x −1 ＋ … ＋ a−1 x ＋ a ＝ 0(3.1)の形の方程式について考えていくことにする．
超越方程式に対しては一般に解が全部でいくつあるかを特定することは難しいのに対し，n 次代数方程式には（重複を含めて）ちょうど n 個の複素数解が存在することがわかっている．これは代数学の基本定理と呼ばれており，史上最大の数学者といわれるガウスによって 18 世紀に証明された．
次に，代数方程式は因数分解により解を見つけることができ，また，逆に解 α が 1 つ見つかれば，f(x) を因数分解することができる．実際，因数定理により f(x) は x − α で割り切れ，これにより，1 つ次数の下がった n − 1
次代数方程式が導かれる．さらには，実数係数の代数方程式の場合，複素数解は必ず共役複素数との対として現れる．
このような性質を有効に用いることにより，代数方程式のすべての解を数値的に求めることが可能となる．

3.1.2

解の公式

まず，1 次方程式 ax ＋ b ＝ 0 の場合，a ≠ 0 であることをチェックした
後，単に x ＝ −ba とすれば解が求まる．次に，2 次方程式
ax 2 ＋ bx ＋ c ＝ 0の場合，解の公式より 2 つの解は(a ≠ 0)3.2

1 次因子の組立除法

x＝−b ±  b 2 − 4ac2a51と表せる．ただし，第 1 章の例 1 でも述べたように，b 2 ≫  4ac  の場合には
桁 落 ち が 生 じ る た め 注 意 が 必 要 で，こ の 場 合 に は 分 母 と 分 子 に −b ∓

 b 2 − 4ac を掛けて，
x＝2c−b ∓  b 2 − 4acのように分子の有理化を行なった公式を併用し，桁落ちが生じないように工夫する．
さらに，3 次方程式にはカルダノの公式，4 次方程式にはフェラーリの公式と呼ばれる公式が知られているので，係数の値を公式に代入することにより，原理的には解を計算できる．しかしながら，これらの公式は煩雑な形をしており，それほど実用的とはいえない．また，5 次以上の代数方程式に対しては，19 世紀の数学者アーベルとガロアによって，一般的な解の公式は存在しないことが示されている．もう少し詳しくいえば，方程式(3.1) に対
a に対する有限回の代数的操作（加減乗除および累乗根）
し，係数 a1，a2，…，では解を表現することはできないことがわかっている．
以上の理由により，数値計算の立場からは，3 次以上の代数方程式に対しては逐次近似法によって解を計算するのが普通である．以下では，代数方程式(3.1) のすべての解を，逐次近似法によって計算する方法について解説し
よう．

3.2

1 次因子の組立除法

3.2.1 複素ニュートン法
こ の 節 で は，a1，a2，…，a は 与 え ら れ た 複 素 数 と し，方 程 式 (3.1) の 解（複素数を含む）の計算法について解説する．
複素数を係数とする方程式を扱うことを明示するために，未知数として

52

第3章

代数方程式

x の代わりに z を用い，次の形の n 次代数方程式
f (z)：＝ z  ＋ a1 z −1 ＋ … ＋ a ＝ 0(3.2)

について考える．この方程式を解くために，まず 1 つの解 α1 を計算する．
そのためには，例えば（複素）ニュートン法などの逐次近似法で計算すればよい．複素数値のニュートン法の公式は実数の場合とまったく同じである．
すなわち，z を複素数値とし，反復z＋1 ＝ z −f (z)f′(z)( ＝ 0，1，2，…)(3.3)

を行なう．ただし，係数 a1，a2，…，
a がすべて実数の場合には，初期値 z0 を実数から選ぶと反復で得られる数列はすべて実数となり，複素数解が計算できないことになるので注意が必要である．
逐次近似法を行なう際には，まず適当な初期値 z0 を選ぶことになるが，うまい初期値を選ぶためには，解が存在する範囲をある程度押さえておくとよい．この目的のためには，次の補題が有用である．
補題

正の数 M > 0 を…， a  }M ：＝ max {  a1 ， a2 ，としたとき，方程式(3.2) の解 α の絶対値  α  は M ＋ 1 未満である．すな
わち，すべての解は複素平面上の{ z ∈ C ： z  < max {  a1 ， a2 ，…， a  } ＋ 1 }の範囲にある．
【証明】 方程式(3.2) の解 α が  α  ≥ M ＋ 1 を満たすと仮定して矛盾を導こう．方程式(3.2) より，α は−α  ＝ a1 α −1 ＋ a2 α −2 ＋ … ＋ aを満たす．この式の両辺の絶対値をとると，M の定義より
 α  ＝  a1 α −1 ＋ a2 α −2 ＋ … ＋ a ≤  a1   α ≤ Mα

−1

−1＋  a2   α ＋ Mα

−2

−2＋ … ＋  a 

＋…＋M

(3.4)

100

第4章

連立 1 次方程式

切な値を試行錯誤により選択する必要がある．
例題 4.4
連立 1 次方程式



2x1 ＋ 3x2 ＋ x3 ＝ 1
2x1 ＋ 4x2 − x3 ＝ −2
3x1 ＋ x2 − 3x3 ＝ −3をヤコビ法，ガウス - ザイデル法，SOR 法で解き，収束の結果を比較せよ．

【解】 この連立 1 次方程式の解は
x1 ＝ 0.434783…，

x2 ＝ −0.391304…，である．反復法で用いる初期値を x

(0)
1

＝ 0，x

(0)
2

x3 ＝ 1.304348…

＝ 0，x3(0) ＝ 0 とする．このとき，ヤコビ法では表 4.4 のような結果となり，ほぼ収束するまで 70 回程度の反復が必要である．
ガウス - ザイデル法での反復では表 4.5 のような結果となり，32 回の反復でほぼ収束している．
SOR 法（ω ＝ 0.8）での反復では表 4.6 のような結果となり，28 回の反復でほぼ収束している．
以上の結果，3 つの方法の中では SOR 法が最も収束が速いことがわかる．

表 4.4

ヤコビ法による計算

反復の
回数 n

x1()

x2()

x3()

0
1
2
3
4
5
⫶
69
70

0.000000
0.500000
0.750000
0.583333
0.520833
0.392361⫶
0.434782
0.434783

0.000000−0.500000−0.500000−0.541667−0.395833−0.409722⫶−0.391305−0.391304

0.000000
1.000000
1.333333
1.583333
1.402778
1.388889⫶
1.304348
1.304348

表 4.5 ガウス - ザイデル法による計算
反復の
回数 n
0
1
2
3
4
5
⫶
30
31

x1()
0.000000
0.500000
1.000000
0.645833
0.310764
0.303964x2()0.000000−0.750000−0.687500−0.380208−0.275608−0.347258⫶⫶
0.434782 −0.391304
0.434783 −0.391304

x3()0.000000
1.250000
1.770833
1.519097
1.218895
1.188211⫶
1.304347
1.304348

4.7

誤差の評価

101

表 4.6 SOR 法による計算 (ω ＝ 0.8 の場合)
反復の
回数 n

x1()

x2()

x3()

0
1
2
3
4
5
⫶
27
28

0.000000
0.400000
0.763733
0.725222
0.566453
0.447115⫶
0.434782
0.434783

0.000000−0.560000−0.623360−0.526983−0.426492−0.375161

0.000000
0.970667
1.438891
1.527427
1.444917
1.346632⫶−0.391305−0.391304⫶
1.304347
1.304348□

この例題 4.4 が示すように，未知数の少ない方程式に反復法を適用してもそれほど効果的ではない．本来，反復法は，未知数が多く，係数行列が疎行列のときにこそ威力を発揮する方法である．

4.7

誤差の評価

4.7.1 誤差と残差
前節までは，連立 1 次方程式に対し，いかに効率良く解を求めるかについて述べてきた．この節では，数値解法で求めた解の誤差評価について述べる．
x で表すと，そ連立 1 次方程式 Ax ＝ b の真の解を x  で，その近似解を の誤差 e はe＝x − x
である．ただし，真の解は実際にはわからないので，近似解からはその誤差が直接的には評価できない．そこで，
r ＝ Ax−bから誤差を評価することを考える．このベクトル r を残差あるいは残差ベク
トルといい，係数行列 A を近似解 x に掛けて b との差をとることで簡単に

102

第4章

連立 1 次方程式計算できる．もちろん，残差ベクトルが 0 のときは 
x は真の解であるが，そうでないときに残差ベクトルから誤差をどのようにしたら評価できるのだ
ろうか．
まず，最も簡単な 1 次方程式 ax ＝ b について考えてみよう．真の解 x  とその近似解 x ＝ x  ＋ e に対し，
r ＝ ax − b ＝ a(x  ＋ e) − b ＝ ae
である．これより誤差は e ＝ ra となるが，これは  a  が非常に大きいと，たとえ残差が大きくても誤差は小さく，逆に  a  が非常に小さいと残差が小
さくても誤差は大きいということである．連立 1 次方程式に対しても，もし係数の絶対値がすべて小さいと同様のことが起こることは容易に想像できる．しかしながら，係数の絶対値が小さくなくても，誤差が残差に比べて非常に大きくなることがある．
次の例を見てみよう．
例2
連立方程式



2x1 ＋ x2 ＝ 1
1.99x1 ＋ 0.99x2 ＝ 1を考える．行列とベクトルを使って Ax ＝ b と表すと
A＝



2

1

1.99 0.99



，

x＝


 
x1
x2

，

b＝


1
1



である．簡単にわかるように，真の解は x ＝ (1，−1) である．
一方，x
 ＝ (0，1) として誤差と残差を計算してみると，それぞれ
e＝
x − x ＝

 
−1
2

，

r ＝ Ax
−b＝



0
−0.01



となり，残差に比べて誤差が非常に大きい．これは，係数行列が非正則行列

 
2 1
2 1

に近いために，係数行列の逆行列の成分が大きくなり，単独の 1 次方程式において係数が小さい場合と同様のことが生じたためである．

□

116

第5章

行列の固有値問題b2 ＋ b2 ＋ b2 ＋ b2 ＝ a2 ＋ a2 ＋ a2 ＋ a2

(5.6)

および
b2 ＋ b2 ＝ a2 ＋ a2

( j ≠ p，q)を得る．したがって，（ i ) が成り立つ．
また，(5.5) のときには，b ＝ b ＝ 0
である．これと (5.6) より，b2 ＋ b2 ＝ a2 ＋ a2 ＋ a2 ＋ a2である．よって，（ ii ) が成り立つ．

5.2.2

□

ヤコビ法

定理 5.2（ i ) より，変換行列 R(p，q；θ) による相似変換によっては，成分の 2 乗和は変化しない．一方，定理 5.2（ ii ) より，θ を (5.5) のように選べ
ば，非対角成分の 2 乗和は小さくなる．したがって，A よりも B ＝ R  AR の方が対角行列に近くなり，修正された行列 B もまた実対称行列である．そこで，(p，q) をいろいろ変えて行列 R(p，q；θ) を繰り返し A に作用させる
ことにより，A を徐々に対角行列に近づけることができる．このようにして，固有値と固有ベクトルを逐次近似によって求める方法をヤコビ法という．
ヤコビ法による実対称行列の固有値と固有ベクトルの計算
Step 1 ：A0 ＝ A とする．

Step 2 ：1，
2，…) の非対角成分のうち，絶対値が最大の成A ( ＝ 0，q) 成分に対する回転行列を R とし，A＋1
分 a を選ぶ．(p，＝ R A R とおき，これを繰り返す．
Step 3 ：A＋1 の非対角成分の絶対値の最大値が十分小さくなったら反復を打ち切り，A＋1 の対角成分を固有値の近似値とし，U ＝ R0 R1 R2 … Rの各列を固有ベクトルの近似ベクトルとする．

5.2

実対称行列の固有値問題

117

例題
ヤコビ法を用いて次の行列を対角化し，固有値と固有ベクトルを求めよ．



2
1
4
5 −2
A＝ 1
4 −2 −2

【解】 まず，A0 ＝ A とする．A0 の非対角成分のうち，絶対値が最大の成分は (1，3)
あるいは (3，1) 成分であるから，行列 R0 をR0 ＝とする．ここでθ0 ＝

cos θ 0 0 sin θ 0010−sin θ 0 0 cos θ 0



−2a13
1
1
arctan＝ arctan (−2)
22a11 − a33

と選ぶと，



4.472136 −0.200811
0.000000
A1 ＝ R0 A0 R0 ≃ −0.200811
5.000000 −2.227033
0.000000 −2.227033 −4.472136

と計算される．A1 の非対角成分のうち，絶対値が最大の成分は (2，3) あるいは (3，2)
成分なので



1
R1 ＝ 0とする．ここでθ1 ＝

00
cos θ 1 sin θ 1

0 −sin θ 1 cos θ 1



−2a23
1
arctan2
a22 − a33と選び，A2 ＝ R1 A1R1 を計算する．
同様の手順を繰り返すと



4.436151 0.000000
0.000000
A10 ≃ 0.000000 5.533668
0.000000
0.000000 0.000000 −4.969818
および



118

第5章



行列の固有値問題

0.859661 −0.042144 −0.509124

0.959887
U10 ＝ R0 R1 R2 … R10 ≃ 0.176219
0.479511 −0.277202が得られる．

0.218091
0.832603



よって，行列 A の固有値は A10 より
λ1 ≃ 4.436151，

λ2 ≃ 5.533668，

λ3 ≃ −4.969818と計算され，対応する固有ベクトルは U10 より

0.859661
v1 ≃ 0.176219 ，0.479511と計算される．

5.3

v2 ≃





−0.042144
0.959887 ，−0.277202

v3 ≃

−0.509124
0.218091
0.832603

□

一般の行列の固有値問題

この節では，必ずしも実対称ではない行列の固有値問題を数値解析的に解く方法について解説しよう．一般の行列に関する固有値問題では，すべての固有値を求めるというよりは，特別な固有値，例えば絶対値最大の固有値，
絶対値最小の固有値，ある値に最も近い固有値などを求めることを目的とするのが普通である．この節では，一般の行列に対するこのような特別な固有値の計算方法について述べる．

5.3.1

支配的固有値とべき乗法

A を n 次正方行列とし，（説明を簡単にするために）行列 A は n 個の実固有値 λ1，λ2，…，λ をもち，これらは
 λ1  >  λ2  ≥  λ3  ≥ … ≥  λ 

(5.7)

v2，
を満たしていると仮定する．また，各固有値に対応する固有ベクトル v1，…，v ∈ R は 1 次独立であると仮定する．このとき，R の基底として固v2，…，v) をとることができる．
有ベクトルの組 (v1，
絶対値が最大の固有値 λ1 のことを行列 A の支配的固有値といい，これに6.3

中 点 公 式131 p (x) dx ＝   ∑ f (x )L (x) dx ＝ ∑  L (x) dx f (x )










＝0





＝0









が得られる．そこで，I[ f ] の近似公式として
I[ f ] ＝ ∑ α f (x)， α ＝
＝0

 L (x) dx (i ＝ 0，1，…，n) (6.6)


を選ぶと，f(x) が n 次以下の多項式であれば正しい積分値が計算できることになる．したがって，積分公式 (6.6) の次数は少なくとも n 次である．
以上の議論から，分点が与えられているときには重み α を (6.6) のように定めればよいことがわかる．残るのは分点の選び方である．
分点を等間隔に選ぶ公式としては，ニュートン - コーツ型積分公式 (6.3
〜 6.5 節) と呼ばれる一連の公式がある．特に，数表などのように，関数
値が等間隔の分点における数値データとして与えられている場合には，必然的にニュートン - コーツ型積分公式が使われることになる．また，関数f (x) が具体的な数式で与えられていて，任意の点における f(x) の値が計算可能な場合には，必ずしも分点を等間隔にする必要はない．分点が等間隔であるという条件をはずして最適に選び，次数をできるだけ高くしようという方向の公式としてガウス型積分公式 (6.6 節) がある．

6.3

中点公式

6.3.1 中 点 公 式
まず，最も簡単な形のニュートン - コーツ型積分公式である中点公式について説明しよう．中点公式は分点が 1 個の公式であり，分点として積分区間[a，b] の中点を選んだものである．
中点公式
I0[f ] ＝ (b − a)f (x0)，

図 6.2 に中点公式の考え方を示す．

x0 ＝

a＋b
2

132

第6章

定
（b）
f

（x
f 0）

積

分

（x）
f

（a）
f

a

b

x0

x

図 6.2 中点公式

これは，f(x) の変動を無視して f(x) を一定値 f(x0) で近似したものと考b] 上のどのえることができる．この考え方に基づけば，分点 x0 は区間 [a，点でもよいのであるが，後で説明するように，分点として中点を選ぶことによって誤差を小さく（正確には公式の次数を高く）できる．
なお，中点公式は積分区間の両端の値 f(a) および f(b) を用いていない．
このように，両端を分点としない積分公式を開型積分公式と呼ぶ．
中点公式の次数は 1 である．すなわち，被積分関数が 1 次関数 f(x) ＝
cx ＋ d のときには，簡単な計算から
I [f ] ＝
＝







(cx ＋ d)dx ＝




1 2
cx ＋ dx
2






1
a＋b
c(b 2 − a 2) ＋ d(b − a) ＝ (b − a) c⋅
＋d
2
2



＝ (b − a)f (x0) ＝ I0[f ]となり，正しい積分値を与えていることが確かめられる．
分点として中点以外の点を選ぶと，正確な積分値が得られるのは f(x) が定数関数のときに限られ，積分公式の次数は 0 となる．中点を分点として選んだのは，公式の次数を高くできるという理由があったのである．

6.3.2

複合中点公式

次に，複合中点公式について述べよう．複合積分公式を導くには，h ＝

6.3

中 点 公 式133(b − a) n として積分区間 [a，b] を n 個の小区間[a ＋ ih，a ＋ (i ＋ 1)h](i ＝ 0，
1，…，n − 1)に分割し，各小区間に中点公式を適用して足し合わせる．複合中点公式の場
合，分点を



x ＝ a ＋ i ＋



1
h
2

(i ＝ 0，1，…，n − 1)として，I[f ] ＝ h f (x0) ＋ h f (x1) ＋ … ＋ h f (x−1)＝ h { f (x0) ＋ f (x1) ＋ … ＋ f (x−1)}とすればよい．これより，次の積分公式が得られる．
複合中点公式I[f ] ＝ h(f 0 ＋ f 1 ＋ … ＋ f −1)，x ＝ a ＋ i ＋



1
h，
2

f  ＝ f (x)

h＝

b−a
n

(i ＝ 0，1，…，n − 1)

複合中点公式は積分区間の両端の値 f(a) および f(b) を用いていないので，開型公式の一種である．また，分点は各小区間の中点であるが，これは数値データを用いる場合には都合の良いことではない．
例えば，関数 f(x) の値が数値データとして x ＝ 0 から x ＝ 1 まで 0.01刻みで与えられているとしよう．区間 [0，1] 上の定積分を複合中点公式で計
fn − 1

（x）
f

…
f1
f0

a

x0

x1

…
h

xn − 1 bx

図 6.3

複合中点公式

7.6

連立常微分方程式と高階常微分方程式

189

は，n 個の未知関数からなる 1 階の連立常微分方程式に変形できる．
例題 7.5
平らな板の上に重りをおき，左端の壁とバネでつなぐ（図 7.4 参照）．
水平な座標軸を設定し，バネが自然な位置（伸びても縮んでもいない状態）の重りの位置を x ＝ 0 とする．

x
O
図 7.4

x（t）

バネの振動の概念図

重りを右に移動させて手を離すと，重りが運動を始める．
時刻 t での重りの位置を x(t) とするとき，常微分方程式
m x″(t) ＝ − x(t) − c x′(t)が成立する．ここで，m > 0 は重りの質量， > 0 はバネ定数（バネを単位
長さ伸ばすのに必要な力），c ≥ 0 は運動系にはたらく抵抗係数である．バ
ネから受ける力はバネの伸びに比例し（フックの法則），抵抗は重りの速度に比例すると仮定した．ここで，y(t) ＝ x′(t)，a ＝ m > 0，b ＝ cm ≥ 0
とおくと，



x′
(t) ＝ y(t)y′(t) ＝ −a x(t) − b y(t)を得る．これを 4 次のルンゲ - クッタ法で数値計算せよ．

【解】 ステップ幅は h ＝ 0.01 とし，常微分方程式内の定数は a ＝ 1，b ＝ 0，0.1，
2.5 とした計算結果を図 7.5 に示す．b は抵抗を表すので，抵抗がない場合，小さい場合，大きい場合の 3 通りを計算することになる．
この図より，b ＝ 0 のときは周期的に運動し，b ＝ 0.1 では徐々に振幅が小さくなる減衰振動であり，b ＝ 2.5 では振動することなく単調に重りが x ＝ 0 の位置へ戻ることが観察できる．なお，（減衰）振動をするかしないかの境目は，数学的な解析から b ＝  4 a であることが示される．

190

第7章

常微分方程式b＝0

b ＝ 0.1

0.4

0.4

x
（t）

0.8

x（t）

0.8

0

0

−0.4

−0.4

−0.8

−0.8
10

0

20

30

40

50

0

10

20

t
b ＝ 2.5

1

t

30

40

50

x
（t）

0.8
0.6
0.4
0.2
0

0

2

4

6

8

t

10

12

14

16

図 7.5 4 次のルンゲ ‑ クッタ法によるバネの振動のシミュレーション．定
数 a の値はいずれも a ＝ 1 で，ステッ
プ幅は h ＝ 0.01 である．

□

7.7
7.7.1

非自励系
非自励系の常微分方程式

ここまでは，自励系の常微分方程式を扱ったが，本節では非自励系の常微分方程式の初期値問題

x′(t) ＝ f (t，x(t))x(t0) ＝ x0

(7.36)を扱う．
自励系では，初期時刻が変わったとしても，解を初期時刻の差の分だけ時間シフトすれば，解の様相には変化がない．一方，非自励系では初期時刻の値が重要である．初期時刻を例えば t0 ＝ 0 とした場合と t0 ＝ 1 とした場合
では，初期値が同じであっても解の様相が変わりうる．

7.7

非 自 励 系

191

この節では，非自励系の常微分方程式の解法として，オイラー法，2 次のルンゲ - クッタ法，4 次のルンゲ - クッタ法について簡潔に述べる．自励系と非自励系の数値解法の考え方は同様であるが，非自励系の解法の導出のためには偏微分が現れる．そのため，偏微分に関する知識を必要とする読者
は，最低限のことを次章の 8.1.1 項にまとめておいたので，まずは目を通して頂きたい．

7.7.2 非自励系に対する数値解法
オイラー法の考え方は単純である．すでに述べた自励系に対する解法を
x＋1 ＝ x ＋ h f (t，x)(i ＝ 0，1，…，
N − 1)と変更するだけである．なぜならば，オイラー法とは (7.36) の常微分方程式の左辺の x′(t) を差分を用いて近似を行なうというもので，右辺に対して
は f(t，
x) の t における値をそのまま使うからである．
次に，(7.36) に対する 2 次と 4 次のルンゲ - クッタ法を扱うため，陽的な 1 段法の公式 (7.14) をx＋1 ＝ x ＋ h Φ(t，x，h)(i ＝ 0，1，…，
N − 1)と修正する．この公式に対して，2 次のルンゲ - クッタ法は (7.21) と同様に，

Φ(t，x，h) ＝ α1 ＋ β21 ＝ f (t，x)
2 ＝ f (t ＋ ph，x ＋ qh1)とする．ここで，定数 α，β，p，q は，次数が 2 となるように，次のように決める．
まず，テイラー展開より，
2 ＝ f (t ＋ ph，x ＋ qh1)＝ f (t，x) ＋ ph f (t，x) ＋ qh1 f (t，x) ＋ O(h 2)
である．一方，非自励系常微分方程式に対する初期値問題 (7.36) の解 x(t)
に対して，

208

第8章

偏微分方程式I，J の値を 100 〜1000 程度にすれば，係数行列 A は 10000 〜1000000 次程度の大きさになる．このような大規模な連立 1 次方程式をコンピュータで解くときには多くの記憶容量を必要とするため，SOR 法（4.6.4 項参照）などの反復法を使うのが一般的である．
例題 8.2
領域を Ω ＝ (−1，1) × (−1，1) とすると，関数u(x，y) ＝ cos πx2 (1 ＋ y) (1 − y)3はディリクレ境界条件を満たす．そこで，この u がポアソン方程式を満たすように関数 f(x，y) を与える．このとき，問題 B を差分法で解き，そのときの誤差について調べよ．

【解】 ここで，x，y の分割数 I，J を I ＝ J ＝ N として，分割数 N と数値解の相対誤差
e(N )：＝

max  u − u(x，y ) 

0≤，≤max  u(x，y ) 

0≤，≤の関係を調べる．すると，表 8.2 のような計算結果が得られる．
表 8.2 2 次元領域における差分法による誤差
分割数 N

相対誤差 e(N)e(N) Δx 2

10
20
40
80

9.16152 × 10−3

0.229038
0.223557
0.224005
0.224003

2.23557 × 10−35.60013 × 10−4
1.40002 × 10−4

これより，誤差はほぼ N 2 に反比例（Δx 2 ＝ Δy 2 にほぼ比例）することがわかる．
□

例題 8.2 の連立 1 次方程式は SOR 法で解いた．N ＝ 80 での計算において，その加速パラメータ ω の値を ω ＝ 1 とした SOR 法では，解に収束する

8.3

熱方程式に対する差分法

209までの反復回数は 12324 回であったが，ω ＝ 1.94 のときは 338 回であり，
約 36 倍の差がある2)．この差は，ほぼ計算時間の差と等しい．したがって，
SOR 法は ω の値により計算時間が大きく改善されることに注意しよう．

8.3

熱方程式に対する差分法

8.3.1 初期条件と境界条件
未知変数の値が場所 x と時刻 t で決まるものを時間発展問題と呼ぶ．本節では，その中で特に重要なものの 1 つである熱方程式の数値解法について述
べる．
熱方程式は
∇

u ＝ κ

uと表せる．未知関数 u(x，t) が温度を表すものとすると，κ > 0 は熱の伝わりやすさを表す量であり，拡散係数と呼ばれる．以下では簡単のため κ ＝ 1とするが，κ が 1 以外の定数のときには，空間変数 x を変数変換する（言い換えれば，空間の単位を変える）ことにより，いつでも κ ＝ 1 とできることを注意しておく．
以下では，領域が閉区間 [a，b] の場合を扱い，このとき，熱方程式はu ＝ u(8.9)と表せる．境界条件としてディリクレ境界条件を採用すると，次のように定式化できる．
問題 C
T ] → R，u：[0，T ] → R，u0：[a，b] → R と定数 T > 0
関数 u：[0，が与えられたとき，

2)収束判定条件により収束までの反復回数は変わるため，これらの回数は参考値として見てほしい．

210



第8章

偏微分方程式u ＝ u(x ∈ (a，b)，t ∈ (0，T ])u(a，t) ＝ u(t)，u(b，t) ＝ u(t)u(x，0) ＝ u0(x)(t ∈ (0，T ])(x ∈ (a，b))(8.10)

を満たす関数 u：[a，b] × [0，T ] → R を求めよ．

関数 u，u はポアソン方程式の場合と同じく，ディリクレ境界条件における境界値を与えるものである．関数 u0 は時刻 t ＝ 0 における解の値を指定する初期値である．問題 C は初期条件と境界条件から解を求める問題であり，初期値 - 境界値問題と呼ばれる（図 8.4 参照）．
t＝T

u ＝ ua

t＝0
x＝a

8.3.2

vt ＝ Duu ＝ u0u ＝ ubx＝b

図 8.4 ディリクレ境界条件下での1 次元熱方程式

各種の差分法

初期値境界値問題 (8.10) を差分法で数値的に解くために，次のような格子を導入する（図 8.5）
．自然数 I，N ∈ N と実数 Δt > 0 を与え，x ＝ a ＋ i Δx

(i ＝ 0，1，…，I )t ＝ n Δt(n ＝ 0，1，…，N ； t  ≤ T )とおく．ここで，Δx：＝ (b − a) I である．以下では，格子点 (x，t) ＝ (x，t)t) の近似値を u() で表す．
における解 u(x，
この近似解を次のような条件から定める．まず，初期条件より，
u(0)：＝ u0(x)(i ＝ 0，1，…，I )