数値解析についての資料
一橋大学経営管理研究科 小林健太

2023/7/20 12:05 更新

0

この資料の概要

現代においてコンピュータシミュレーションは，自然科学や工学，社会科学から医療に至るまで様々な分野で用いられ，我々の生活を基礎から支えている．講義では，コンピュー
タシミュレーションを行う際，現象を記述する方程式をどのようにして計算機上で計算すればよいのか，また，元の問題の解と計算機で計算した解にはどのような関係があるのか等，数値計算にまつわる諸問題について説明する．
この資料はhttp://kobayashi.hub.hit-u.ac.jp/lecture/からダウンロードできる．随時，更新する予定．

0.1

プログラミング

数値解析においては，理論の理解も大切だが，実際にプログラムを組んで理論を確かめることも重要である．この資料では，プログラムは C 言語（C++）を用いる．C++は C 言語の拡張となっており，C 言語のプログラムは C++としても解釈できるが，C++は C 言
語としては実行できない．この資料では，プログラムは基本的には C++の文法で記述する．ただ，慣れの問題もあると思うので，C 言語の関数も適宜用いている（printf 関数など）．パソコン上で C++を使うにはコンパイラをインストールする必要があるが，フリーのものでも各種の選択肢がある．
初心者にとっては Embarcadero 社の Free C++ Compiler がインストールも簡単で良いだろう．

https://www.embarcadero.com/jp/free-tools/ccompiler/free-download
こちらからダウンロードできる．ダウンロードには登録が必要だが，登録しても宣伝メールが送られてくるとかの不都合は特に生じない．インストールの方法について

http://kobayashi.hub.hit-u.ac.jp/topics/compiler.html

1

に簡単に解説しておいたので，参考にしてもらいたい（ただし，情報は古い可能性がある）．

Visual Studio は開発環境などが色々と揃っていて，Windows アプリケーションなど、そこそこ本格的なプログラミングも可能である．フリーだが登録が必要．
Cygwin は Windows 上の仮想 Linux 環境だが，インストールすると C++コンパイラが付いてくる．サーバー管理等，Linux に興味のある人はこれも良いだろう．
Cygwin からコンパイラのみを抜き出したものとして MinGW がある．非常に軽くてお勧
めだが，若干，インストールに苦労することがあるかもしれない．
他にも，Open Watcom C/C++，Digital Mars C/C++ Compilers，なんてのもある．

Mac については，私は詳しくは知らないが，色々とやり方はあるようだ．Windows でな
く，最初から Linux 機でプログラムするのも手かもしれない．

0.2

線形演算ライブラリ

C 言語には，行列やベクトルを扱う機能が基本機能としては付属していないので，それを実現するライブラリを用意した．
http://kobayashi.hub.hit-u.ac.jp/topics/linear.html
ここから Linear.h というファイルをダウンロードしてプログラムと同じフォルダに入れておくこと．その上で，プログラムの最初に#include "Linear.h"と書いておくと行列やベクトルを扱うことができる．

0.3

グラフ表示

数値計算結果をグラフ化して表示するために，gnuplot をインストールする必要がある．

http://www.gnuplot.info/
こちらから Download from SourceForge のリンクを辿るとダウンロードサイトに行けるが，バナー広告で「Download」とか書かれているのが出たりするので間違ってクリック
しないように．

32bit ウィンドウズ機なら「gp***-win32-setup.exe」を，64 ビットウィンドウズ機なら「gp***-win64-setup.exe」を選ぶと良いだろう（ただし，***はバージョン番号で数字が
入る）．

2

0.4

この資料のプログラミング例について参考のためにこの資料にもプログラムの例を載せてあるが，あくまで参考という事で，すべての数値解析手法についてプログラムが載っているわけではない．プログラムが載っていないものについては，各自で頑張ってプログラミングしてもらいたい．

0.5

参考文献

C 言語の参考書としては
柴田望洋「新版 明解 C 言語 入門編」（ソフトバンククリエイティブ）
皆本晃弥「やさしく学べる C 言語入門−基礎から数値計算入門まで」（サイエンス社）
を勧める．
数値解析について，参考にしたい，もしくはもっと詳しく勉強したいという人には
森正武「数値解析」（共立出版）
齊藤宣一「数値解析入門」（東京大学出版会）
皆本晃弥「C 言語による数値計算入門」（サイエンス社）などの書籍を勧める．理解もさらに深まることだろう．

0.6

質問・連絡先など

この資料について，何か質問や連絡したいことがある場合は，私のメールアドレス

kenta.k@r.hit-u.ac.jp
まで連絡すること．

3

1

連立一次方程式

何か現象を表すモデル方程式があったとして，そのモデル方程式を離散化して計算機で解こうとすると，連立一次方程式が現れることが多い．そこで，まずは連立一次方程式をいかにして解くか，について説明しよう．以下では，i 行 j 列の要素が aij であるような行列を (aij ) と書く．

1.1

ガウスの消去法

まず連立一次方程式の解法として思いつくのは，学部一年生の線形代数で習った掃き出し法である．これは，数値解析ではガウスの消去法として知られている．
最初に具体例を考えてみよう．連立一次方程式
   −1 −2 4x1
4   
7 −2 x2  = −5 2
−3 −8 6x3
14
を解きたいとする．まず，1 行目の定数倍を他の行から引いて，2 行目以降の 1 列目を零にする．
   4−1 −2 4x1
   
3
6  x2  =  3
 02
x3
0 −2 −6
右辺のベクトルも変化していることに注意．線形代数で，右辺ベクトルを行列にくっつけて拡大係数行列を作ったことを思い出そう．次に，2 行目の定数倍を他の行から引いて，3行目以降の 2 列目を零にする．
   −1 −2 4x1
4   
3
6  x2  =  3
 0
0
0 −2
x3
4
このようにして，まずは上三角行列が得られた．この操作を前進消去という．
前進消去が終わったら，x3 から x1 まで順に解を求めていく．

1· 4 = −2−21x2 = (3 − 6x3 ) = 5
31x1 =(4 − (−2)x2 − 4x3 ) = −22
−1
x3 =

4

これで解が x1 = −22, x2 = 5, x3 = −2 と求まる．このように，上三角行列から代入によって解を求める操作を後退代入という．
これを一般的に書くと以下のようになる．解きたい連立一次方程式を



(1)

(1)

(1)

a11
a12
a13
(1)
(1)
 a(1)
a22
a23
 21
 (1)
(1)
(1)
 a31
a32
a33
 .
..
..
 .
.
.
 .
 (1)
(1)
(1)
an−1 1 an−1 2 an−1 3
(1)
(1)
(1)
an1
an2
an3

(1)

···
···
···
...

a1 n−1(1)
a2 n−1(1)
a3 n−1
..
.

···
···

an−1 n−1(1)
an n−1(1)  (1) (1)  a1nb1
x1(1)  (1)
a2n   x2   b2   (1) (1)  a3n   x3   b3 

 

.. 
 .  =  . .   ..   ..   (1) (1)an−1 n  xn−1  bn−1 
(1)
(1)
xn
bn
ann

とし，k = 1, 2, 3, · · · , n − 1 について
(k)

aik

(k)
a
(k) kjakk(k)aik (k)(k)(k+1)= bi − (k) bk
bi
akk
(k+1)

aij

(k)

= aij −(i = k + 1, · · · , n,

j = k + 1, · · · , n)(i = k + 1, · · · , n)と計算すると，連立一次方程式は

 (1) (1) (1)a11 a12 a13
 0 a(2) a(2)

22
23

(3)
 0
0 a33
 .
..
..
 .
.
.
 .

 0
0
0
0
0
0

(1)

···
···
···
..
.

a1 n−1(2)
a2 n−1(3)
a3 n−1
..
.

···
···

an−1 n−10(n−1)

  (1) (1)  
b1
a1n
x1
(2)
(2)  
a2n   x2   b2   (3) (3)  a3n   x3   b3 


 
.. 
 .  =  . .   ..   ..   (n−1) (n−1)  an−1 n  xn−1  bn−1 
(n)
(n)
xn
bn
ann

と変形され，前進消去が実行される．この連立一次方程式を改めて

 


b̃1
x1
ã11 ã12 ã13 · · · ã1 n−1
ã1n

 
 0 ã
ã2 n−1ã2n   x2   b̃2 
22 ã23 · · ·

 


0
0 ã33 · · · ã3 n−1ã3n   x3   b̃3 

 

 .
..
.. . .
..
..  .
 .  =  . .
.
.
.
.   ..   .. 
 .

 


0
0
0 · · · ãn−1 n−1 ãn−1 n  xn−1  b̃n−1 
b̃n
xn
0
0
0 ···
0
ãnn
と書けば，

xn =

1
b̃n
ãnn

5

xn−1 =

1ãn−1 n−1b̃n−1 − ãn−1 n xn

···
1
x2 =
ã22

b̃2 −

n
X

!
ã2j xj

j=3

1
x1 =
ã11

b̃1 −

n
X



!
ã1j xj

j=2

によって後退代入が実行され，解が求まる．
以上の説明では，行列やベクトルの要素が変化するときに，変化前と変化後を右上の添え字で区別したが，実際のプログラムでは単に代入して値を更新していけばよい．

1.2

ピボット選択(k)(k)

ガウスの消去法においては，途中で akk が零になると計算が破綻してしまう．また，akk
が零でなかったとしても，非常に小さな値のときには，計算精度が悪くなることがある．
まずは実例でその様子を見てみよう．
以下の連立一次方程式を考える．解は x1 = x2 = x3 = 1 である．
  
4.501
0.001 2.000 2.500
x1  
−1.000 1.648 4.745 x2  = 5.393
9.703
x3−2.000 4.438 7.265
この方程式をガウスの消去法で，4 桁の精度（丸めは四捨五入）で計算してみよう．1 段
目の前進消去が終わると以下のようになる．
  
4.501
0.001 2.000 2.500
x1  
2002 2505  x2  =  4506  004004 5007
x39012

2 段目の前進消去として，2 行目を 2 倍して 3 行目から引くと
  

4.501
x1
0.001 2.000 2.500
  

2002 2505  x2  =  4506 
 0
0
x3
0
0
−3.000となる．引き続いて後退代入を行うと，解が x1 = −1.000, x2 = 2.251, x3 = 0.000 と得られるが，これは本当の解とは大きく異なっている．

6このような現象が起きる原因は，今回の場合，前進消去の第 1 段目にある．対角成分が0.001 と非常に小さいので，1 行目の定数倍を 2 行目と 3 行目から引く際に，絶対値の大
きな数と小さな数の演算が生じ，小さい方の数の情報が失われるからである．このような現象を情報落ちという．
(k)

ガウスの消去法において，akk をピボットという．情報落ちが出来るだけ生じないようにするには，計算の過程でピボットの絶対値が大きい方がよい．そこで，前進消去の k 段目
(k)
(k)
(k)
を行う前に，|akk |, |ak+1 k |, · · · , |ank | を比べて一番大きなものを探し，その行と第 k 行を入れ換えると精度が良くなることが期待できる．このような操作を，ピボット選択という．
先ほどの例では，まず第 1 行と第 3 行を入れ換えて  

−2.000 4.438 7.265
x1
9.703

  

−1.000 1.648 4.745 x2  = 5.393
0.001 2.000 2.500
x34.501
としてから計算を開始することになる．このようにして解くと，解が x1 = 1.002, x2 =1.001, x3 = 0.9998 となり，4 桁の精度で計算していることを考えれば妥当な結果が得ら
れる．

1.3

ガウスの消去法のプログラム

ガウスの消去法のプログラムを以下に示そう．プログラムは C++の機能も用いているので，ファイル名は”〜.c”ではなく”〜.cpp”とする必要がある．
線形演算ライブラリ Liner.h では，既にガウスの消去法の計算が実装されており，行列 Aに大して A.Gauss(b) とすることで，Ax = b の解を求めることができる．そこで，既に実装されているガウスの消去法の結果と比較して答え合わせを行っている．

#include <stdio.h>#include <math.h>#include "Linear.h"
int main(){
printf("***** ガウスの消去法 *****\n");次ページへ続く…

7

int n = 3;
Matrix A(n, n); // n × n 次元行列を宣言
Vector b(n);// n 次元ベクトルを宣言
A[1][1] = 2; A[1][2] = 4; A[1][3] = 5; b[1] = 1;
A[2][1] = 1; A[2][2] = 6; A[2][3] = 1; b[2] = 2;
A[3][1] = 6; A[3][2] = 2; A[3][3] = 7; b[3] = 3;
Matrix A2 = A; // 答え合わせ用にコピーしておく
Vector x = b;// 前進消去
for(int k=1; k<=n-1; k++){// ピボット選択
int p = k;double pv = fabs(A[k][k]);for(int i=k+1; i<=n; i++){double a = fabs(A[i][k]);if(a>pv){
p = i; pv = a;}}if(p>k){for(int i=k; i<=n; i++){double a = A[k][i]; A[k][i] = A[p][i]; A[p][i] = a;}double a = x[k]; x[k] = x[p]; x[p] = a;}for(int i=k+1; i<=n; i++){double a = A[i][k] / A[k][k];for(int j=k+1; j<=n; j++){A[i][j] -= a*A[k][j];}x[i] -= a*x[k];}}次ページへ続く…

8

// 後退代入
for(int k=n; k>=1; k--){for(int i=k+1; i<=n; i++){x[k] -= A[k][i]*x[i];}x[k] /= A[k][k];}printf("解は以下のようになりました\n");x.Show(16); // 有効数字 16 桁で表示
printf("答え合わせ\n");x = A2.Gauss(b);x.Show(16);return 0;}
プログラムが正しければ，以上のように出力される筈である．

***** ガウスの消去法 *****
解は以下のようになりました
0.7872340425531914
0.2553191489361702-0.3191489361702127
答え合わせ0.7872340425531914
0.2553191489361702-0.3191489361702128

1.4

定常反復法

零でない要素（非零要素）の個数が行列全体の要素の個数（n 次行列なら n2 個）に比べて非常に少ない行列を疎行列という．方程式を離散化すると，疎行列を係数に持つ連立一次方程式になることが多い．次数の非常に大きな連立一次方程式を計算機で解くとき，行列の成分を全てメモリに保持するのは現実的ではない．例えば，100 万次元の問題を解く
場合，100 万× 100 万の変数をメモリに確保しようとすると，8T（テラ）バイトのメモリが必要になってしまう．非零要素が対角要素周辺に集まっている場合は，ガウスの消去法

9を用いても新たな非零成分は対角要素の周辺にしか現れないので，対角成分の周辺のみメモリに確保すればよいが，疎行列であっても非零成分が行列全体に散らばっている場合には，ガウスの消去法を用いると新たな非零成分が行列全体に多数発生してしまい，メモリ消費の点で問題が生じる．そこで用いられるのが反復法である．ここでは，反復法の中でも定常反復法という種類に属する，ヤコビ法，ガウス・ザイデル法，および SOR 法について紹介する．
連立一次方程式の i 行目は

ai1 x1 + ai2 x2 + · · · + ain xn = bi
であるが，これは aii が零でない場合には

1
xi =
aii

bi −

X

!
aij xjj̸=iと変形できる．そこで，

1
(k+1)
=
xi
aii

bi −

X

!
(k)

,

aij xji = 1, 2, · · · , nj̸=iなる反復を考えてみる．このような反復によって解を求める方法をヤコビ法という．
少々天下りに思えるかもしれないが，例えば高校でやった漸化式

an+1 = pan + q
を思い出すと，この数列は |p| < 1 のときに収束し，収束先を a とすると

a = pa + qを満たすので，それとの類似で考えれば，ヤコビ法の反復は，対角成分 aii が非対角成分
aij (i ̸= j) に比べて大きいときに収束することが期待できる．実際，対角成分が大きいよ
うな例
!
!
!
3 1
x1
4
=
−1 5
x2
4
(0)

(0)

について，初期値 x1 = 0, x2 = 0 に対してヤコビ法を適用すると

1
(1)
(0)
x1 = (4 − x2 ) = 1.33333 · · ·
3
1
(2)
(1)
x1 = (4 − x2 ) = 1.06666 · · ·
3

1
(1)
(0)
x2 = (4 + x1 ) = 0.8
5
1
(2)
(1)
x2 = (4 + x1 ) = 1.06666 · · ·
5

10

1
(3)
(2)
x1 = (4 − x2 ) = 0.97777 · · ·
3
1
(4)
(3)
x1 = (4 − x2 ) = 0.99555 · · ·
3
1
(5)
(4)
x1 = (4 − x2 ) = 1.00148 · · ·
3

1
(3)
(2)
x2 = (4 + x1 ) = 1.01333 · · ·
5
1
(3)
(4)
x2 = (4 + x1 ) = 0.99555 · · ·
5
1
(5)
(4)
x2 = (4 + x1 ) = 0.99911 · · ·5と，解 x1 = 1, x2 = 1 に近づいていくことがわかる．
方程式を離散化すると，比較的対角成分の大きな行列が得られることが多いので，ヤコビ法は有効に機能することが多い．ただ，いつも収束するとは限らない．収束のための厳密な十分条件については後ほど説明する．
(k)

(k)

(k)

(k+1)

(k+1)

(k+1)

ヤコビ法では，x1 , x2 , · · · , xn を全て求めてから，それらを元に x1, x2, · · · , xn
(k+1)
(k+1)
(k+1)
を計算している．しかし，xiを計算しようとする時点では既に x1, · · · , xi−1 が求
(k)
(k)
(k+1)
(k+1)
, · · · , xi−1 を使おうというまっているので，x1 , · · · , xi−1 よりも最新の結果である x1のが，ガウス・ザイデル法である．具体的に書くと
!
i−1
n
X
X
1
(k+1)
(k)
(k+1)
bi −
aij xj−aij xjxi=,
i = 1, 2, · · · , n
aii
j=1
j=i+1
なる反復を行うことになる．ガウス・ザイデル法は xi について常に最新の結果だけを保持していればいいので，計算上は同じベクトルに上書きしていけばよく，プログラムは簡単になる．
ガウス・ザイデル法は

1
(k+1)
yi
=
aii

bi −

i−1
X

(k+1)

aij xj

−

j=1

n
X

!
(k)

aij xj

,

j=i+1

(k)
(k+1)
(k)
(k+1)
− xi ),xi= xi + (yi
(k)

(k+1)

i = 1, 2, · · · , n
(k)

(k+1)

と書き換えることができる．すなわち，xi に (yi− xi ) を足したものが xiている．このとき，足す量を緩和係数 ω で調整して
!
i−1
n
X
X
1
(k+1)
(k+1)
(k)
yi
=
bi −
aij xj−aij xj
,
aii
j=1
j=i+1
(k+1)

xi

(k)

(k+1)

= xi + ω(yi(k)− xi ),になっ

i = 1, 2, · · · , nとしたものが SOR 法である．緩和係数を調整することにより，ガウス・ザイデル法よりも高速に収束させることができる．緩和係数は 0 < ω < 2 の範囲に，もっとも高速に収束する値があるということが知られているが，問題に応じて色々と試して決定する必要がある．ω = 1 のとき，SOR 法はガウス・ザイデル法に一致する．

11

ヤコビ法および 0 < ω ≤ 1 についての SOR 法はX
A = (aij ),|aij | < |aii |,i = 1, 2, · · · , nj̸=iを満たすような行列を係数に持つ連立一次方程式に対して収束することが知られている（ガウス・ザイデル法は SOR 法に含まれる）．このような行列を狭義優対角行列という．
狭義優対角行列以外にもヤコビ法や SOR 法が収束するケースは色々とあるのだが，一般的な収束証明は難しいので，ここでは狭義優対角行列の場合にヤコビ法が収束することを証明してみよう．

1.5

ベクトルのノルムヤコビ法の収束に関する証明のためには，ベクトルや行列の大きさを測る基準が必要になる．収束というのは，近づいていくことであるが，近づいていくというのは差が小さくなっていくことである．ということは，小さいとか大きいとかを判断する基準が必要になる．そこで，ノルムというものを導入する．なお，ベクトルや行列の要素などは全て実数の範囲で考える．
任意のベクトル x, y に対して以下の 3 条件を満たすような実数値関数 ∥ · ∥ をノルムという．

(1)
(2)
(3)

正値性
∥x∥ ≥ 0 であり， ∥x∥ = 0 ⇐⇒ x = 0
同次性
任意の実数 α に対して ∥αx∥ = |α|∥x∥
三角不等式 ∥x + y∥ ≤ ∥x∥ + ∥y∥さらに，1 ≤ p ≤ ∞ とし，
x1
 x2  
x= .  .. 
xnに対して∥x∥p =


!1/p
n

X



|xi |p
1≤p<∞
i=1



 max |xi |
1≤i≤n

p=∞

と定義すると，∥ · ∥p はノルムとなる．
ここで，p = ∞ のときだけ式の形が違っていて，おやっ，と思う人がいるかもしれない
が，p → ∞ のときの極限を考えれば p = ∞ のときの形になる．実際，M = max |xi | と
1≤i≤n

12

すると，
n
X

M = (M p )1/p ≤!1/p|xi |p≤ (nM p )1/p = n1/p M −→ M (p → ∞)i=1が成り立つ．よって挟み打ちの原理より
n
X

!1/p
|xi |p→ max |xi |(p → ∞)

1≤i≤n

i=1

となる．
さて，∥ · ∥p がノルムになることを p = 1, 2, ∞ の場合に証明しておこう（ここでは証明しないが，一般の p ≥ 1 場合もそれほど難しくない）．ノルムの条件のうち，(1) と (2) は難しくないので，(3) の三角不等式のみ示す．

 y1
 y2  y=.
 .. 


x1
 x2  
x =  . , .. 

ynxnとおくと，p = 1 のときは

∥x + y∥1 =nX|xi + yi | ≤

i=1

n
X

n
X

(|xi | + |yi |) =i=1

|xi | +

i=1nX|yi | = ∥x∥1 + ∥y∥1

i=1

p = ∞ のときには∥x + y∥∞ = max |xi + yi | ≤ max (|xi | + |yi |) ≤ max |xi | + max |yi | = ∥x∥∞ + ∥y∥∞
1≤i≤n

1≤i≤n

1≤i≤n

1≤i≤n

より成り立つ．p = 2 のときは少し難しいが，コーシー・シュワルツの不等式を用いると∥x + y∥22 =nX|xi + yi | ≤
2

i=1

≤

n
X
i=1

n
X
i=1

|xi |2 +nXi=1= (∥x∥2 + ∥y∥2 )|xi | +
2nX

|yi | + 2
2

i=1

n
X

|xi ||yi |

i=1

v
v
v
2
u n
u n
u n
n
X
X
X
X
u
u
u
|yi |2 + 2t
|xi |2|yi |2 = t|xi |2 + t|yi |2 
i=1

i=1

2

と証明できる．

13

i=1

i=1

1.6

行列のノルム行列に対してもノルムを定義する．いま，∥ · ∥ をベクトルのノルムとするとき，行列 A の
ノルムを
∥Ax∥
∥A∥ = maxx̸=0 ∥x∥と定義する．とくに，1 ≤ p ≤ ∞ に対して

∥Ax∥px̸=0 ∥x∥p∥A∥p = maxと定義する．

行列ノルムの意味は，∥Ax∥ が ∥x∥ に比べてどれくらい大きくなるか考え，その比率が最大になるときの比率を行列の大きさとする，ということである．
行列ノルムは以下の性質を満たす．ただし，O は零行列，I は単位行列である．

(1)
(2)
(3)
(4)
(5)
(6)

∥A∥ ≥ 0 であり， ∥A∥ = 0 ⇐⇒ A = O
任意の実数 α に対して ∥αA∥ = |α|∥A∥∥A + B∥ ≤ ∥A∥ + ∥B∥∥Ax∥ ≤ ∥A∥∥x∥∥AB∥ ≤ ∥A∥∥B∥∥I∥ = 1(1),(2),(3),(6) は難しくない．(4) は x ̸= 0 のとき
∥Ax∥
∥Ay∥
≤ max
= ∥A∥
y̸=0 ∥y∥∥x∥となることからすぐにわかり，(5) は∥ABx∥∥A∥∥Bx∥∥Bx∥≤ max= ∥A∥ max
= ∥A∥∥B∥x̸=0
x̸=0x̸=0 ∥x∥∥x∥∥x∥∥AB∥ = max
よりわかる．

p = 1 のときと p = ∞ のときは，行列ノルムは行列の要素を用いて具体的に書くことができる．すなわち，n 次正方行列 A = (aij ) について次が成り立つ．
(1)(2)∥A∥1 = max

n
X

1≤j≤n

∥A∥∞ = max

1≤i≤n|aij |

i=1nX|aij |

j=1

14

以下，




x1
 x2  
x= .  .. 
xnとする．まず (1) を示そう．αj =nX|aij | とおくと

i=1

∥Ax∥1 ==nnXX
i=1 j=1nX

aij xj ≤nn XXαj |xj | ≤ max αj
1≤j≤n

j=1

|aij | |xj | =i=1 j=1
n
X

n
n
X
X
j=1

!
|aij | |xj |

i=1|xj | = max αj ∥x∥1
1≤j≤n

j=1

よって

max αj ∥x∥1
∥Ax∥1
1≤j≤n
≤ max
= max αj∥A∥1 = max
1≤j≤n
x̸=0
x̸=0 ∥x∥1∥x∥1
一方で，α1 , α2 , · · · , αn のうちで最大のもの（最大のものが複数ある場合はそのうち一つ）を αk とし，
 
0
 
0
.
 ..  
v= 1 ← k 番目
.
.
.
0
とすると，

n
X

|aik |
∥Ax∥1
∥Av∥1
i=1
∥A∥1 = max≥== αk = max αj1≤j≤nx̸=0 ∥x∥1
∥v∥1
1
よって

∥A∥1 = max αj = max
1≤j≤n

1≤j≤n

n
X

|aij |

i=1が言えた．

(2) については，αi =nX|aij | とおくと

j=1

∥Ax∥∞ = max

1≤i≤n

n
X
j=1

aij xj ≤ max

1≤i≤n

n
X
j=1

15

|aij | |xj | ≤ max

1≤i≤n

n
X
j=1

|aij | · max |xj |1≤j≤n

= max αi max |xj | = max αi ∥x∥∞
1≤i≤n

1≤j≤n

よって

1≤i≤n

max αi ∥x∥∞
∥Ax∥∞
1≤i≤n
≤ max
= max αi1≤i≤nx̸=0 ∥x∥∞
x̸=0
∥x∥∞

∥A∥∞ = max

一方で，α1 , α2 , · · · , αn のうちで最大のもの（最大のものが複数ある場合はそのうち一つ）を αk とし，
 v1(
 v2 
1(akj ≥ 0) 
v =  . ,
vj =
 .. 
− 1 (akj < 0)

vn
とすると，
n
X

max
aij vj
1≤i≤n
∥Ax∥∞
∥Av∥∞
j=1
∥A∥∞ = max≥=x̸=0 ∥x∥∞
∥v∥∞
1
n
n
X
X
≥
akj vj =|akj | = αk = max αi
j=1

1≤i≤n

j=1

よって

∥A∥∞ = max αi = max
1≤i≤n

1≤i≤n

n
X

|aij |

j=1

が言えた．

1.7

ヤコビ法の収束証明

ヤコビ法の反復は



0

 a21

a
 22
 a31

M = −a
 33
 .
 ..


 an1
ann

x(k+1) = M x(k) + c,a1n 
a12 a13···a11 a11
a11 
a23
a2n 

0
···
a22
a22 

a32
a3n 

0 ···
,
a33
a33 

..
..
.. 
...
.
.
. 



an2 an3···0
ann ann

16




b1
 a11 
 b2 a  22  b c= 3  a33  .. 
 . 




 bn 
annと書ける．ここで，元の連立一次方程式の解を x∗ とすると

x∗ = M x ∗ + c
より

x(k+1) − x∗ = M (x(k) − x∗ )が成り立つ．したがって

x(k) − x∗ = M k (x(0) − x∗ )より

∥x(k) − x∗ ∥∞ = ∥M k (x(0) − x∗ )∥∞ ≤ ∥M ∥k∞ ∥x(0) − x∗ ∥∞が成り立つので，∥M ∥∞ < 1 ならばヤコビ法が収束するということがわかる．
実際に計算すると∥M ∥∞ = max

1≤i≤n

X aij
j̸=i

aii

1 X
|aij |
1≤i≤n |aii |
j̸=i

= max

となり，A の狭義優対角性から ∥M ∥∞ < 1 がいえる．よって，狭義優対角行列を係数に持つような連立一次方程式に対してヤコビ法が収束するということが証明できた．

1.8

ヤコビ法のプログラム

参考までにヤコビ法のプログラムを紹介しておく．行列は 10 次元の 3 重対角行列とした．
定常反復法においては，いつ反復を終了するかが重要になるが，ここでは相対残差ノルム∥b − Ax∥∞∥b∥∞
が予め決めたある小さな値より小さくなったときに反復を終了するものとする．以下のプログラムでは，相対残差ノルムが 10−12 以下になったときに反復を終了している．

#include <stdio.h>#include "Linear.h"
int main(){
printf("***** ヤコビ法 *****\n");int n = 10;
次ページへ続く…

17

Matrix A(n, n);
Vector b(n);//行列と右辺ベクトルを設定
for(int i=1; i<=n; i++){for(int j=1; j<=n; j++){if(i==j){A[i][j] = 4;
} else if(i==j+1){
A[i][j] = 1;
} else if(i==j-1){
A[i][j] = -2;} else {A[i][j] = 0;}}b[i] = i;}
Vector x(n);Vector x2(n);// 初期値は 0 ベクトル
for(int i=1; i<=n; i++){x[i] = 0;}//反復開始int f=0; //フラグ: 1=計算成功, 0=計算失敗
int k;
for(k=1; k<=1000; k++){for(int i=1; i<=n; i++){x2[i] = b[i];for(int j=1; j<=n; j++){
次ページへ続く…

18

if(j!=i){x2[i] -= A[i][j] * x[j];}}x2[i] /= A[i][i];}x = x2;// 収束判定
Vector z = b - A*x; // 残差ベクトル
double r = z.NormInfty() / b.NormInfty(); // 相対残差ノルム
if(r < 1.0E-12){
f = 1; // 成功
break;}
if(r > 1.0E100){
break; // 失敗}}
if(f==0){
printf("収束しませんでした\n");} else {printf("反復%d 回目で収束\n", k);printf("解は以下のようになりました\n");x.Show(16);printf("答え合わせ\n");x = A.Gauss(b);x.Show(16);}return 0;}
相対残差ノルムが 10100 を超えたり，1000 回の反復でも終了しなかった場合には，「収束しませんでした」と表示するようにしている．
ガウス・ザイデル法や SOR 法についてもプログラムして，収束の速さを比べてみると良い．

19

1.9

補足

時間の都合で触れなかった事項についてここで述べておく．
一度きりの計算ではガウスの消去法はそこそこ効率的だが，ある決まった行列 A について，複数の右辺ベクトル b(1) , b(2) , · · · , b(l) に対してそれぞれ Ax = b(j) の解を求めたい場
合には，ガウスの消去法より効率の良い手法として LU 分解というアルゴリズムが知られている．
数値解析においては，行列の固有値や固有ベクトルを計算したいことも多い．そのようなときには，コレスキー分解や QR 分解が用いられることが多い．
本章では，連立一次方程式の反復解法としてヤコビ法，ガウス・ザイデル法，SOR 法について述べたが，これらの方法は歴史的には重要だが，実際にはあまり使われない（全く使われないというわけではないが）．そのかわり，反復法としては CG 法という方法が現在よく用いられている．ただ，CG 法については説明が煩雑になるので説明は省略した．

1.10

課題

課題１ ガウスの消去法の前進消去では，行列を上三角行列に変換していくことがポイン
トである．しかし，実際のプログラムでは行列は上三角行列には変形されていない（添え字の動く範囲を確認すること）．それでもちゃんと解が計算できる理由を説明せよ．
課題２ ガウスの消去法のプログラムにおいて，ピボット選択を用いない場合に精度が悪くなるような例を考え，計算結果を比較せよ．
課題３ 10 次元以上の同じ連立一次方程式を，ヤコビ法，ガウス・ザイデル法，パラメー
タを色々と変えた SOR 法で解き，計算結果を比較せよ．
課題４

ヤコビ法は

A = (aij ),X|aij | < |ajj |,j = 1, 2, · · · , ni̸=jを満たす行列に関する連立一次方程式についても収束することを示せ．

202

非線形方程式

本章では，非線形方程式の解法について考える．

2.1

二分法

1 変数の非線形方程式 F (x) = 0 の解を求めたいとする．一番素朴な方法は，中間値の定理（の変形バージョン）を利用する方法だろう．すなわち，F (x) が連続のとき，「a < b について，F (a) と F (b) が異符号もしくはどちらかが 0 ならば，[a, b] に F (x) = 0の解が存在する」という事実を用いる方法である．具体的には，まず，x(0) < x(0) かつ F (x(0) )F (x(0) ) ≤ 0となるような初期値 x(0) , x(0) を求め，k = 0, 1, 2, · · · について以下のような反復を行う．

x(k) + x(k)2(k)
F (x )F (m(k) ) ≤ 0 のとき

m(k) =(1)x(k+1) = x(k) ,(2)

x(k+1) = m(k) ,F (m(k) )F (x(k) ) ≤ 0 のとき
x(k+1) = m(k) ,x(k+1) = x(k) .

(1),(2) の両方が成立するときには，どちらを選んでもよい．
二分法については，方程式の真の解を x̂ とすると，
1  (0)1  (k)
(k)
(0)
(k)
x −x
= k+1 x − x|x̂ − m | ≤22が成り立つ．すなわち，m(k) を解の近似値とすると，反復するごとに誤差は半分になる
ことがわかる．
二分法の利点は，F (x) が微分可能でなくとも，連続でありさえすれば適用可能であること，さらに，解の範囲が明確に得られることである．
二分法を用いて e−x = x の解を求めるプログラムを以下に示す．

21

#include <stdio.h>#include <math.h>double f(double x){
return exp(-x)-x;}
int main(){
printf("***** 二分法 *****\n");// 初期値
double xl = 0;
double xh = 1;
if( f(xl)*f(xh)>0 ){printf("関数値が異符号となるように 2 つの初期値を取って下さい\n");return 0;}
for(int k=1; k<=60; k++){double mid = (xl+xh)/2;
if( f(xl)*f(mid)<=0 ){xh = mid;} else {xl = mid;}
printf("%d 回目の反復の結果、解の範囲は [%.15f,%.15f]\n",k,xl,xh);}
return 0;}

2.2

一次元ニュートン法

次に，二分法より速く収束する計算方法であるニュートン法を紹介する．まずは 1 変数の場合を考えてみよう．いま，非線形方程式 F (x) = 0 の近似解 x̃ があったとして，この近22

似解の精度を改善することを試みる．F の微分可能性を仮定すれば，F を一次の項までテ
イラー展開することによりF (x) ≈ F (x̃) + (x − x̃)F ′ (x̃)
と近似できる．よって，右辺 = 0 を x について解き，x = x̃ −F (x̃)F ′ (x̃)として x を求めれば，x̃ よりも良い近似解になると期待できる．これは要するに，関数F (x) を，x̃ における接線で近似し，解を求めていることになる．よって，適当な初期値
x(0) から出発し，
F (x(k) )x(k+1) = x(k) − ′ (k)F (x )
なる反復を行えば，{x(k) } は解に収束することが期待できる．これが一次元ニュートン法である．ニュートン法は，適当な滑らかさの仮定の下に，初期値を真の解の近くに取れば，解へ収束することが保証されている．初期値が真の解から遠い場合は，収束しないこともある．
それでは，実際にニュートン法を用いて非線形方程式を解いてみよう．例として，二分法のときと同じ e−x = x の数値解を考える．
まず，方程式を F (x) = 0 の形に直す．この例の場合は F (x) = e−x − x とすればよい．よって，ニュートン法の反復は

e−x − x(k)−−e−x(k) − 1
(k)

x

(k+1)

=x

(k)

F (x)

x(k+1) x(k)
図 1: ニュートン法23

となる．実際のプログラムは以下のようになる．

#include <stdio.h>#include <math.h>double f(double x){
return exp(-x)-x;}
double df(double x){
return -exp(-x)-1;}
int main(){
printf("***** 一次元ニュートン法 *****\n");// 初期値
double x = 0;
for(int k=1; k<=10; k++){double diff = f(x)/df(x); // 変化量
x -= diff;printf("x^(%d)=%.15f, 変化量=%.15f\n", k, x, diff);}
return 0;}
これを実行すると

x(1) = 0.537882842739990,

変化量 = 0.462117157260010

x(2)= 0.566986991405413,

変化量 = −0.029104148665423

x(3)= 0.567143285989123,

変化量 = −0.000156294583710

x(4)= 0.567143290409784,

変化量 = −0.000000004420661

x(5) = 0.567143290409784,

変化量 = −0.000000000000000

x(6) = 0.567143290409784,

変化量 = −0.000000000000000

24x(7) = 0.567143290409784,

変化量 = −0.000000000000000

x(8) = 0.567143290409784,

変化量 = −0.000000000000000

x(9) = 0.567143290409784,

変化量 = −0.000000000000000

x(10) = 0.567143290409784,

変化量 = −0.000000000000000となり，誤差が急激に減少することがわかる．

2.3

一次元ニュートン法の収束速度

ここでは一次元ニュートン法がどれくらいの速度で収束するのか調べてみよう．

F (x) = 0 の真の解を x̂ とし，F (x) は x̂ を含むある閉区間 I で C 2 級であるものとする．また，区間 I で F ′ (x) は 0 にならないものとし，
M = max |F ′′ (x)|,m = min |F ′ (x)| > 0

x∈I

x∈I

とする．このとき，x(k) ∈ I とし，F (x) を x(k) のまわりで 2 次までテイラー展開し，x = x̂を代入するとF (x̂) = F (x(k) ) + (x̂ − x(k) )F ′ (x(k) ) +(x̂ − x(k) )2 ′′F (ξ)
2となる．ここで，2 次の項はラグランジュの剰余項であり，ξ は x(k) と x̂ の間の実数である．さて，F (x̂) = 0 より

−F (x(k) ) = (x̂ − x(k) )F ′ (x(k) ) +(x̂ − x(k) )2 ′′F (ξ)
2が成り立つので，これを用いてニュートン法の反復の式を変形すると

x

(k+1)

F (x(k) )− x̂ = x − x̂ − ′ (k)F (x )(x̂ − x(k) )F ′ (x(k) ) + (x̂ − x(k) )2 F ′′ (ξ)/2= x(k) − x̂ +F ′ (x(k) )(x̂ − x(k) )2 F ′′ (ξ)=
2F ′ (x(k) )(k)となる．よって|x(k+1) − x̂| ≤M (k)|x − x̂|22mが成り立つ．これを変形するとM (k+1)|x− x̂| ≤
2m



25

2
M (k)
|x − x̂|2mとなるので，M (k)|x − x̂| ≤
2m



2k
M (0)
|x − x̂|
2m

なる誤差評価が得られる．つまり，|x(0) − x̂| <2mMを満たすように初期値を取れば，ニュートン法は収束することがわかる．さらに，M (k+1)
M (k)
log10
|x
− x̂| ≤ 2 log10|x − x̂|
2m
2m
より，ニュートン法の数値解が真の解に近いときには，一回反復を行うごとに，解の有効桁数が倍になることがわかる．
一般に，解が|x(k+1) − x̂| ≤ C|x(k) − x̂|pなる形で収束するとき，数列 {x(k) } もしくはその反復法は p 次収束するという．ニュートン法は x̂ のまわりで F ′ (x) ̸= 0 のとき，2 次収束することがわかる．

2.4

多次元ニュートン法

次に，一次元ニュートン法の考え方を多次元に拡張する．F (x) を n 変数 n 次元関数とする．ここで， 
f1 (x)x1
 f2 (x)  x2  
F (x) =  . 
x =  . , .. 
 .. 

xnとする．さて，F = 0 の近似解を

fn (x)
x̃1
 x̃2  
x̃ =  . 
 .. 
x̃nとし，この近似解のまわりで F を一次の項までテイラー展開するとF (x) ≈ F (x̃) + [F ′ (x̃)](x − x̃)26と近似できる．ここで，[F ′ (x)] はヤコビ行列であり，∂f1 ∂f1∂f1···
 ∂x1 ∂x2
∂xn 
 ∂f2 ∂f2
∂f2 


···

[F ′ (x)] = 
∂x
∂x
∂x
1
2
n
 .

.
.
.
 .
..
..
.. 
 .

 ∂f

∂f
∂f
n
n
n
···
∂x1 ∂x2∂xnで定義される．よって，一次元ニュートン法の場合と同様に考えると，

x(k+1) = x(k) − [F ′ (x(k) )]−1 F (x(k) )なる反復を行うことで数値解を計算することができる．
多次元ニュートン法を別の視点から考えてみよう．一変数の場合，ニュートン法の反復は，関数を近似解における接線で近似することにより，より良い近似解を求めることであった．
同様に，多変数の場合，ニュートン法は，関数を接平面で近似して近似解を計算していると考えることができる．それを具体的に見ていこう．

fk (x) に x = x̃ で接する接平面を
y = fk (x̃) + ak1 (x1 − x̃1 ) + ak2 (x2 − x̃2 ) + · · · + akn (xn − x̃n )とする．ここで，実際にこれが接平面となっていることから，l = 1, 2, · · · , n について∂y∂fk(x̃) =
= akl
∂xl
∂xl
が成り立つ．よって，この接平面は
∂fk
∂fk
∂fk
y = fk (x̃) +(x̃) · (x1 − x̃1 ) +(x̃) · (x2 − x̃2 ) + · · · +(x̃) · (xn − x̃n )
∂x1
∂x2
∂xn
と書ける．よって，これを fk (x) の近似として用いれば，F (x) はF (x) ≈ F (x̃) + [F ′ (x̃)](x − x̃)と近似できることがわかる．
ニュートン法の反復には

x(k+1) = x(k) − [F ′ (x(k) )]−1 F (x(k) )のようにヤコビ行列の逆行列が出てくるが，実際の数値計算においては，([F ′ (x(k) )]d(k) = F (x(k) )x(k+1) = x(k) − d(k)と考え，第一式の連立一次方程式を解いて x(k) から d(k) を求め，それを元に x(k+1) を求める．逆行列を求めるよりも連立一次方程式を解く方が楽なので，逆行列は計算しないことに注意する．ただし，n = 2 の場合は逆行列を用いてもよい．

27

2.5

多次元ニュートン法のプログラム

3 次元ニュートン法のプログラムは以下のようになる．
#include <stdio.h>#include <math.h>#include "Linear.h"// x^2+y^2+z^2=1// y=sin(x)// z=x+y// を解く。プログラムでは x を x[1]、y を x[2]、z を x[3] とする。
// F(x) = ( x[1]^2+x[2]^2+x[3]^2-1, x[2]-sin(x[1]), x[3]-x[1]-x[2] )// 関数値
Vector f(Vector x){
Vector y(3);y[1] = x[1]*x[1] + x[2]*x[2] + x[3]*x[3] - 1;y[2] = x[2] - sin(x[1]);y[3] = x[3] - x[1] - x[2];return y;}// ヤコビ行列
Matrix df(Vector x){
Matrix y(3, 3);
y[1][1] = 2*x[1];y[1][2] = 2*x[2]; y[1][3] = 2*x[3];y[2][1] = -cos(x[1]); y[2][2] = 1;y[2][3] = 0;y[3][1] = -1;
y[3][2] = -1;
y[3][3] = 1;
return y;}
次ページへ続く…

28

int main(){
printf("***** 多次元ニュートン法 *****\n");int n = 3;
Matrix J(n, n);
Vector x(n);// 初期値x[1]=0.5; x[2]=0.5; x[3] = 0.5;
for(int k=1; k<=10; k++){J = df(x);x = x - J.Gauss(f(x));
printf("x^(%d)=(%.15f,%.15f,%.15f)\n",k,x[1],x[2],x[3]);}return 0;}
計算結果が 2 次収束しているようなら，プログラムは正しいと考えられる．

2.6

補足

本稿では触れなかったが，真の解で一回微分が 0 になったり，ヤコビ行列が非正則になる
ような場合，ニュートン法は二次収束ではなく一次収束することが知られている．
一般的にニュートン法は，初期値を近似解の近くに取れば急速に真の解に近づくが，真の解から離れた初期値では上手く収束しないことが多い．そのため，まずは二分法などの他の方法で，ある程度，真の解に近い近似解を計算し，最後の仕上げにニュートン法を用いることが多い．
解を求めたい方程式が代数方程式，つまりxn + a1 xn−1 + a2 xn−2 + · · · + an−1 x + an = 0なる形の方程式の場合には，ニュートン法ではなく，DKA 法など，代数方程式に特化した方法を用いることが多い．DKA 法を用いると，n 次元方程式の n 個の解を同時に求めることができる．

29

2.7

課題

課題１ sin x = 0 をニュートン法で解くとき，どの範囲に初期値を取れば解 x = 0 に収束するか，数値実験により調べなさい（-1.57〜1.57 の範囲で 0.01 刻みに初期値を取って調べよ）．
課題２ sin x = 0 をニュートン法で解くとき，必ず解 x = 0 に収束するような初期値の範囲を本文を参考に求めなさい（まず区間 I を定める必要がある）．
課題３ 4 変数以上の 4 変数以上の非線形方程式を一つ考え，多次元ニュートン法を用いて解け．その際，収束の様子を観察し，2 次収束していることを確認すること．

30

3

数値積分

本章では，関数 f (x) の積分

Z bQ=

f (x)dxaの計算について考える．積分は，曲線で囲まれた部分の面積を求める場合などに必要になってくる．
関数 f (x) と a, b の値によってはZ 1

1x dx = ,
3
0

Z 2

2

1

dx
= log 2,
x

Z 1

dx
π
=
2
2
−1 1 + xなどのように，公式を用いて数学的に積分値を計算することができるが，公式が無いような場合には，数学的に積分値を厳密に計算することができない場合が多い．そこで，数値的な近似計算を考える．

3.1

積分則

以下，c は a と b の中点，すなわち c =

a+b
とする．
2

積分を近似するため，関数 f (x) について次の 3 種類の近似を考える．

(a) f (x) を，a と b の中点での値 f (c) を取る定数関数 f˜0 (x) = f (c) で近似する．
(b) f (x) を，(a, f (a)) と (b, f (b)) を通る一次関数 f˜1 (x) で近似する．ここでb−xx−af˜1 (x) =f (a) +
f (b)
b−a
b−a
である．実際，x = a, b を代入すれば f˜1 (a) = f (a) かつ f˜1 (b) = f (b) となることが確かめられる．

(c) f (x) を，(a, f (a)) と (c, f (c)) と (b, f (b)) を通る二次関数 f˜2 (x) で近似する．ここで(b − x)(c − x)(b − x)(x − a)(x − c)(x − a)f˜2 (x) =f (a) +
f (c) +f (b)(b − a)(c − a)(b − c)(c − a)(b − c)(b − a)である．実際，x = a, b, c を代入すれば f˜2 (a) = f (a), f˜2 (c) = f (c), f˜2 (b) = f (b) が成り立つことが確かめられる．
これらの近似関数を用いて，積分値 Q について以下の近似を考える．

31

中点則：



Z b
f˜0 (x)dx = (b − a)f (c) = (b − a)f

Ra,b (f ) =
a

a+b
2



台形則：

Z b
f˜1 (x)dxaZZf (a) bf (b) b=(b − x)dx +(x − a)dx
b−a a
b−a a

b − a
f (a) + f (b)=2Ta,b (f ) =シンプソン則：Z b
f˜2 (x)dx

Sa,b (f ) =
a

Z b
f (c)
(b − x)(c − x)dx +(b − x)(x − a)dx
(b − c)(c − a) a
a
Z b
f (b)
+
(x − c)(x − a)dx
(b − c)(b − a) a

Z b
Z b
2
=
(b − x)(x − a)dx
(b − x)(c − x)dx + 2f (c)f (a)(b − a)2aa
Z b(x − c)(x − a)dx
+f (b)
a




b−a
a+b
=
f (a) + 4f
+ f (b)62f (a)=(b − a)(c − a)Z b

これらの近似積分公式を一般的に積分則という．実際にこれらの積分則を使用する際には，積分区間を複数の区間に分割し，個々の区間ごとに積分則を適用する．

3.2

積分則の誤差ここでは，積分則による近似の誤差を解析しよう．簡単のため，f (x) を平行移動して，中点が x = 0 になるようにしたものを φ(x) とする．すなわち


a+b
+x
φ(x) = f
2

32

b−a
とおくと
2


Z b
a+b
Q − Ra,b (f ) =f (x)dx − (b − a)f
2
a
Z k
=
φ(x)dx − 2kφ(0)である．このとき，k =

−k

Z b


b − a
f (a) + f (b)
2
a
Z k


=
φ(x)dx − k φ(k) + φ(−k)
−k




Z b
b−a
a+b
Q − Sa,b (f ) =f (x)dx −f (a) + 4f
+ f (b)
6
2
a
Z k

k
φ(k) + φ(−k) + 4φ(0)=φ(x)dx −3−k
Q − Ta,b (f ) =f (x)dx −となる．そこでZ x
g0 (x) =g1 (x) =g2 (x) =

Z−x
x
Z−x
x
−x

φ(t)dt − 2xφ(0)
φ(t)dt − x φ(x) + φ(−x)φ(t)dt −x
φ(x) + φ(−x) + 4φ(0)3と置いて，|g0 (k)|, |g1 (k)|, |g2 (k)| の大きさを考えることにする．
いま，f (x) ∈ C 2 [a, b] とし，M2 = max |f ′′ (x)| = max |φ′′ (x)| とする．このとき，−k≤x≤ka≤x≤bg0′ (x) = φ(x) + φ(−x) − 2φ(0),
Z x
′′
′
′
φ′′ (t)dt
g0 (x) = φ (x) − φ (−x) =−xが成り立つので，0 ≤ x ≤ k について
Z x
Z x
′′
′′
|g0 (x)| =φ (t)dt ≤M2 dt = 2M2 x
−x
−x
Z x
Z x
′
′
′′
|g0 (x)| = g0 (0) +g0 (t)dt ≤2M2 t dt = M2 x2
Z 0x
Z0 x
1
M2 t2 dt = M2 x3
g0′ (t)dt ≤|g0 (x)| = g0 (0) +
3
0
0

33

が成り立つ．よって

1(b − a)3|Q − Ra,b (f )| = |g0 (k)| ≤ M2 k 3 =max |f ′′ (x)|3
24 a≤x≤bとなり，中点則の誤差評価ができた．
次に





g1′ (x) = φ(x) + φ(−x) − φ(x) + φ(−x) − x φ′ (x) − φ′ (−x)
Z x


′
′
= −x φ (x) − φ (−x) = −xφ′′ (t)dt−xが成り立つので，0 ≤ x ≤ k について
Z x
Z x
′
′′
|g1 (x)| = x
φ (t)dt ≤ x
M2 dt = 2M2 x2
−x
−x
Z x
Z x
2
′
|g1 (x)| = g1 (0) +g1 (t)dt ≤2M2 t2 dt = M2 x3300が成り立つ．よって

2(b − a)3|Q − Ta,b (f )| = |g1 (k)| ≤ M2 k 3 =max |f ′′ (x)|3
12 a≤x≤bとなり，台形則の誤差評価ができた．台形則より中点則の方が精度が良いというのは意外かもしれない．
シンプソン則については，f (x) ∈ C 4 [a, b] とし，M4 = max |f ′′′′ (x)| = max |φ′′′′ (x)| と
a≤x≤b

すると，

 x

′
′
φ(x) + φ(−x) + 4φ(0) −φ (x) − φ (−x) ,
3
3




2
x ′
=
φ(x) + φ(−x) − 2φ(0) −φ (x) − φ′ (−x) ,
3
3
 1
 x

2 ′
′
′
′
′′
′′
′′
g2 (x) =φ (x) − φ (−x) −φ (x) − φ (−x) −φ (x) + φ (−x) ,
3
3
3




1 ′
x ′′
=
φ (x) − φ′ (−x) −φ (x) + φ′′ (−x) ,
3
3
 1
 x

1  ′′
′′
′′
′′
′′′
′′′
′′′
φ (x) + φ (−x) −φ (x) + φ (−x) −φ (x) − φ (−x) ,g2 (x) =
3
3
3
Z x


x
x
= − φ′′′ (x) − φ′′′ (−x) = −φ′′′′ (t)dt
33 −x
1g2′ (x) = φ(x) + φ(−x) −−k≤x≤kが成り立つので，0 ≤ x ≤ k について
Z x
Z
x x
x
2
′′′′
′′′
φ (t)dt ≤|g2 (x)| =M4 dt = M4 x2
3 −x
3 −x
3

34

Z x

Z x

2
2
M4 t2 dt = M4 x3
3
9
Z 0x
Z 0x
2
1
|g2′ (x)| = g2′ (0) +g2′′ (t)dt ≤M4 t3 dt = M4 x4
9
18
Z0 x
Z 0x
1
1
|g2 (x)| = g2 (0) +g2′ (t)dt ≤M4 t4 dt = M4 x5
90
0
0 18

|g2′′ (x)| =g2′′ (0) +g2′′′ (t)dt≤が成り立つ．よって|Q − Sa,b (f )| = |g2 (k)| ≤1(b − a)5M4 k 5 =max |f ′′′′ (x)|902880 a≤x≤bとなり，シンプソン則の誤差評価ができた．

3.3

複合積分則

実際に積分の近似値を計算する際には，積分区間を複数の区間に分割し，個々の区間ごとに積分則を適用する．簡単のため，区間は N 等分するものとし，区間幅を h とする．具
体的には，
b−a
,
xk = a + khh=Nとする．そのうえで，各区間で積分則を適用した結果は次のようになる．
複合中点則：Rh (f ) = h

N−1Xf (xk+1/2 )= h f (x1/2 ) + f (x3/2 ) + f (x5/2 ) + · · · + f (xN −1/2 )k=0

複合台形則：N −1h X
Th (f ) =f (xk ) + f (xk+1 )
2 k=0

h
f (x0 ) + 2f (x1 ) + 2f (x2 ) + · · · + 2f (xN −1 ) + f (xN )=2

複合シンプソン則：
N −1

h X
f (xk ) + 4f (xk+1/2 ) + f (xk+1 )Sh (f ) =
6 k=0

35

=


h
f (x0 ) + 4f (x1/2 ) + 2f (x1 ) + 4f (x3/2 ) + 2f (x2 ) + · · · + 4f (xN −1/2 ) + f (xN )6となる．
複合積分則の誤差は各区間の誤差の合計で評価することができるので

h3(b − a)h2
max |f ′′ (x)| =max |f ′′ (x)|a≤x≤b
24 a≤x≤b24h3(b − a)h2|Q − Th (f )| ≤ N ·max |f ′′ (x)| =max |f ′′ (x)|a≤x≤b
12 a≤x≤b12(b − a)h4
h5
max |f ′′′′ (x)| =max |f ′′′′ (x)||Q − Sh (f )| ≤ N ·2880 a≤x≤b2880 a≤x≤b|Q − Rh (f )| ≤ N ·と誤差評価が得られる．すなわち，刻みを細かくしていった場合の誤差限界は，複合中点
則と複合台形則は刻み幅の 2 乗，複合シンプソン則は刻み幅の 4 乗のオーダーで誤差が減少することがわかる．

3.4

複合積分則のプログラム

複合中点則のプログラムは以下のようになる．

#include <stdio.h>#include <math.h>// 積分したい関数
double f(double x){
return exp(x);}// 原始関数
double F(double x){
return exp(x);}次ページへ続く…

36

// [a,b] における 2 階微分の絶対値の最大値
double M2(double a, double b){
return exp(b);}int main(){
printf("***** 複合中点則 *****\n");double a = 1.0, b = 2.0; // 積分区間
for(int N=5; N<=160; N*=2){ // 分割数
double h = (b-a)/N;// 区間幅double S = 0;
for(int k=0; k<=N-1; k++){double x = a + (k+0.5)*h;S += f(x);}S *= h;double Q = F(b) - F(a); // 積分の理論値
double E = (b-a)*h*h/24 * M2(a,b); // 理論的な誤差限界
printf("N=%d、誤差は%.15f、理論誤差限界は%.15f\n", N, S-Q, E);}
return 0;}
積分する関数を変えるときは，最初の 3 つのユーザー定義関数を変更する．
プログラムが間違っていなければ，刻みが倍になるごとに誤差が 4 分の一になる．
複合台形則のプログラムは以下のようになる．

37

#include <stdio.h>#include <math.h>// 積分したい関数
double f(double x){
return exp(x);}// 原始関数
double F(double x){
return exp(x);}// [a,b] における 2 階微分の絶対値の最大値
double M2(double a, double b){
return exp(b);}int main(){
printf("***** 複合台形則 *****\n");double a = 1.0, b = 2.0; // 積分区間
for(int N=5; N<=160; N*=2){ // 分割数
double h = (b-a)/N;// 区間幅double S = 0;
for(int k=0; k<=N; k++){double x = a + k*h;if(k==0 || k==N){
S += f(x);} else {次ページへ続く…

38

S += 2*f(x);}}S *= h/2;
double Q = F(b) - F(a); // 積分の理論値
double E = (b-a)*h*h/12 * M2(a,b); // 理論的な誤差限界
printf("N=%d、誤差は%.15f、理論誤差限界は%.15f\n", N, S-Q, E);}
return 0;}
N が同じ場合，中点則と比較して誤差がおよそ倍になっていればよい．
複合シンプソン則のプログラムは以下のようになる．

#include <stdio.h>#include <math.h>// 積分したい関数
double f(double x){
return exp(x);}// 原始関数
double F(double x){
return exp(x);}// [a,b] における 4 階微分の絶対値の最大値
double M4(double a, double b){
return exp(b);}次ページへ続く…

39

int main(){
printf("***** 複合シンプソン則 *****\n");double a = 1.0, b = 2.0; // 積分区間
for(int N=5; N<=160; N*=2){ // 分割数
double h = (b-a)/N;// 区間幅double S = 0;
for(int k=0; k<=2*N; k++){double x = a + k*0.5*h;if(k==0 || k==2*N){
S += f(x);} else if(k%2==0){
S += 2*f(x);} else {S += 4*f(x);}}S *= h/6;
double Q = F(b) - F(a); // 積分の理論値
double E = (b-a)*h*h*h*h/2880 * M4(a,b); // 理論的な誤差限界
printf("N=%d、誤差は%.15f、理論誤差限界は%.15f\n", N, S-Q, E);}
return 0;}
プログラムが間違っていなければ，刻みが倍になるごとに誤差が 16 分の一になる．

3.5

周期関数に対する複合台形則

滑らかな周期関数に対しては，複合台形則の精度が非常に良くなることが知られている．
この事実は複素関数論の理論により示すことができるのだが，ここでは概略のみ説明する．

f (x) を，周期 2π の滑らかな周期関数とする．ここで，滑らかというのは「解析的」という意味で用いている．
「解析的」という概念については，ここでは詳細は述べないので，詳40

しく知りたい人は複素関数論を勉強してもらいたい．ただ，[0, 2π] の全ての点について定義され，絶対値などを含まない普通の関数で定義されるような関数は，大抵は解析的な周期関数となる．例えば，

f1 (x) = sin x,

f3 (x) = ecos x ,f2 (x) = cos x,

f4 (x) = log

2 − sin x
1 + cos2 xなどは解析的な周期関数である．一方，

f5 (x) = tan x,(
f8 (x) = x (−π ≤ x < π),1f7 (x) =,1 − cos x

f6 (x) = | sin x|,

f8 (x + 2π) = f8 (x)などは解析的ではない．
さて，一定の条件を満たす周期関数は

f (x) =a0
+ a1 cos x + a2 cos 2x + a3 cos 3x + a4 cos 4x + · · ·2+ b1 sin x + b2 sin 2x + b3 sin 3x + b4 sin 4x + · · ·とフーリエ級数展開できることが知られている．特に，f (x) が解析的な周期関数の場合には，ある正の定数 C1 , C2 が存在して，係数 an , bn について

|an | + |bn | ≤ C1 e−C2 nが成り立つことが知られている．
それでは，上のようにフーリエ級数展開されている関数 f (x) の一周期積分を考えてみる
ことにしよう．三角関数の一周期積分は 0 なので，積分すると a0 だけが残り
Z 2π
Q=
f (x)dx = πa00となる．一方で，複合台形則の区間分割数を N とすると，整数 m が N の倍数ではないと
きには
N
−1
X

cos mxk +i

N
−1
X

k=0

k=0

sin mxk =N−1X

eimxk =

k=0
2πmi/N N=

1 − (e)2πmi/N
1−e

N
−1
X

N
−1
X

=

N
−1
X
k=0
2πmi

1−e
=0
1 − e2πmi/N

より

k=0

cos mxk =k=041

e2πmki/N

sin mxk = 0

が成り立ち，m が N の倍数のときにはN−1X

cos mxk + i

k=0

N
−1
X

sin mxk =

k=0

N
−1
X

e

imxk

=

N
−1
X

e

2πmki/N

k=0

k=0

=

N
−1
X

1=N

k=0

より
N
−1
X

N
−1
X

cos mxk = N,

k=0

sin mxk = 0

k=0が成り立つことから，複合台形則による値は

2π X
Th (f ) =f (xk ) = πa0 + 2π(aN + a2N + a3N + · · · )
N k=0
N −1

となる．ここで，f (x) が解析的であって，フーリエ級数の係数について|an | + |bn | ≤ C1 e−C2 nなる評価が成り立つとすれば，|Q − Th (f )| = 2π|aN + a2N + a3N + · · · |≤ 2πC1 (e−C2 N + e−2C2 N + e−3C2 N + · · · )=2πC1 e−C2 N
1 − e−C2 Nとなり，区間分割 N に対して指数的に誤差が減少することがわかる．
実際に，

Z 2π
0

2π
dx
=√
4 + 2 sin x + cos x11などの例を台形則で計算してみると，誤差の指数的減少が確かめられる．
通常の閉区間上の積分も，両端をつないで周期関数とみなしたときの滑らかさが上がるほど，積分則の精度が向上することが知られている．例えば，関数f (x) =

1 (13 − 3x)x
+
x
8

9 ′
1
, f (1) = f ′ (2) = − を満たすので，区間 [1, 2] の両端をつないで周期48関数とみなしたものは，繋いだ点で値と一回微分が連続になる．この関数の積分
Z 2
25
Q=
f (x)dx = log 2 +16
1

は f (1) = f (2) =42

を複合台形公式で計算すると，誤差は h4 に比例するスピードで減少し，理論的に求めた
複合台形則の誤差限界よりも良い精度が実現している．
滑らかな関数 f (x) に対する無限積分

Z ∞Q=

f (x)dx−∞についても，複合台形則は良い精度を与えることが知られている．以下ではそれを説明しよう．
適当な L > 0 を取って∞Xf˜(x) =f (x + Lk)k=−∞とすれば f˜(x) は周期 L の滑らかな周期関数となるので，複合台形則

h=

L
,
N

は

h

N
−1
X

f˜(kh)

k=0

Z L
f˜(x)dx0を高精度に近似する．ここで，

h

N
−1
X

f˜(kh) = h

k=0

Z L

∞
X

f˜(x)dx =f (kh),

が成り立つので，無限積分

f (x)dx
−∞

0

k=−∞

Z ∞

Z ∞
Q=

f (x)dx−∞

は複合台形則h∞Xf (kh),

k=−∞を用いて高精度に計算することができる．
無限積分の場合，複合台形則は無限和になるが，f (x) が e−x のように |x| → ∞ で急激に減少するような関数の場合は，無限和を有限和で打ち切ることで近似値を精度良く計算することができる．その際，打ち切りによる誤差が十分小さくなるように配慮する必要が
ある．
2

43

3.6

補足

本章では，複合積分則により，等間隔の分点における関数値を用いて積分を計算したが，分点を等間隔ではない特別な位置に取ることで，さらに精度を上げることができる．そのような積分公式で有名なものにガウス・ルジャンドル積分公式があり，実際にもよく用いられているが，時間の都合で説明は省略した．
既に見たように，積分する関数の両端をつないで周期関数とみなしたとき，つないだ点における滑らかさが上がると複合台形則の精度が上昇する．これはフーリエ級数論やオイ
ラー・マクローリンの総和公式などを用いて証明することができるので，興味のある人は調べてみれば面白いと思う．解析的な周期関数の場合の複合台形則の誤差の指数的減少については，既に書いたように，複素関数論の知識を用いて証明できる．これも興味のある
人は勉強してみるとよいと思う．
置換積分を利用して積分の積分変数を変換すると，被積分関数の形を変えることができる．その事実を利用し，置換積分により被積分関数を |x| → ∞ で急激に減少する関数に変換して，有限項で打ち切った複合台形則で積分を計算すると，精度良く計算できる．このような考えに基づく数値積分は，変数変換型数値積分と呼ばれ，有名なものに二重指数型変数変換積分公式（DE 公式）がある．DE 公式の誤差は，多くの積分において，分点数に対して指数的減少に近いスピードで減少する．

3.7

課題

課題１ 原始関数のわかっている関数を自由に選び，複合中点則，複合台形則，複合シンプソン則で計算を行い，理論的な誤差限界と比較しなさい．
Z 2 −xe
課題２
dx の値を小数点以下 8 桁まで正確に計算しなさい．小数点以下 8 桁まで
x1
正確な値であるということは，積分則の誤差評価を用いて確かめなさい．積分則の誤差評
価を行う上で，2 階微分や 4 階微分の最大値が必要だが，これは正確な最大値でなくとも，最大値を上から押さえる値であれば誤差評価できることに注意すること．

44

4

常微分方程式

一変数関数とその導関数からなる方程式を常微分方程式という．本章では，常微分方程式の初期値問題の数値解法について説明する．

4.1

色々な常微分方程式の初期値問題

微分方程式とは，関数そのものや，その微分などの間に成り立つ関係式であるが，それだ
けでは解が一つに定まらない．そこで，ある地点での関数値や微分の値を与えて解を考えることが多い．このような条件を初期条件という．解が決まるには，通常は，微分方程式に出てくる微分の階数分だけ初期条件が必要である．
微分方程式の解を y(x) としたとき，y を従属変数，x を独立変数という．これらの用語は，x に応じて y が決まることから名づけられている．以下，微分の記号 ′ や ′′ や (n) などは独立変数による微分を表すものとする．
それでは，微分方程式にどのようなものがあるか見ていくことにしよう．まず，世の中には数学的に（紙と鉛筆で）解ける微分方程式がある．例えば(
y ′ = y,y(0) = 1
は数学的に解くことができて，解は

y = exとなる（独立変数は式に現れていないが，x としている）．ここで y(0) = 1 が初期条件である．初期条件という用語は，時間変数を独立変数に取ることが多いという理由による．
微分方程式の右辺には独立変数が入ってもよい．例えば(
y ′ = xy,y(0) = 1の解は2

y = ex /2となる．
初めにも述べたが，微分の階数が高いと，それに応じて初期条件が必要である．例えば(y ′′ + 2y ′ + 2y = 0,y(0) = 1, y ′ (0) = 0

45

は数学的に解けて，解はy(x) = (cos x + sin x)e−xとなる．以上のような微分方程式は数学的に解が求まるが，世の中には式変形では解けない方程式も多い．
以下は，ロトカ・ボルテラ方程式という，捕食者と被食者（例えばライオンとシマウマ）の増減関係をモデル化した方程式である．
(
x′ = (a − by)x,y ′ = (−c + dx)y
ここで，x と y は時間 t の関数（つまり，t が独立変数，x と y が従属変数）であり，x が被食者の，y が捕食者の個体数を表す．a, b, c, d は正の定数である．初期条件としては，適当な時刻（例えば t = 0）における x と y の値を与える．
以下は，ローレンツ方程式という，大気変動モデルを元にした方程式である．
′x = −px + py,

y ′ = −xz + rx − y,z ′ = xy − bz

ここで，t が独立変数，x, y z が従属変数であり，p, r, b は定数である．ローレンツの原論文では，p = 10, r = 28, b = 8/3 と取られており，現在でもこの値を用いることが多い．初期条件としては，適当な時刻における x, y, z の値を与える．ローレンツ方程式
は，カオスの発見につながった方程式として有名である．
ロトカ・ボルテラ方程式やローレンツ方程式の微分の階数は一階であるが，従属変数の数
は複数である．しかし，これらの方程式は，変数をベクトルとして考えれば，従属変数が一つの場合と同様，(
y ′ = f (x, y),y(x0 ) = y0の形で書くことができる．例えば，ロトカ・ボルテラ方程式では!
!
xx(a − by)y=,f (t, y) =yy(−c + dx)と置けば

y′ = f (t, y)の形になり，t と y を改めて x と y と書き直せば

y ′ = f (x, y)46

の形になる．
以下は，ファン・デル・ポール方程式という，非線形振動を表す微分方程式である．

x′′ − µ(1 − x2 )x′ + x = 0
ここで，x は時間 t の関数，µ は定数である．初期条件としては，ある時刻における x お
よび x′ の値を与える．この微分方程式は，電気回路の振動現象をモデル化するために用いられた．
一般に，x を独立変数とする n 階の微分方程式

y (n) = f (x, y, y ′ , y ′′ , · · · , y (n−1) )
について



y
y′
..
.





y′
y ′′
..
.




 


 





d 

=

 
dx 

 (n−2)  y (n−1)

 
y
′ ′′
(n−1)
(n−1)
f (x, y, y , y , · · · , y)yが成り立つので，各階数の微分をそれぞれ別の関数と考えて yk = y (k) とし，




y0
y1




y2


 y1 


 . 
.

. 
..
F (x, y) = y= . ,




yn−1


yn−2 yn−1
f (x, y0 , y1 , y2 , · · · , yn−1 )と置くことにより，元の微分方程式はy′ = F (x, y)の形に書くことができる．例えば，ファン・デル・ポール方程式は
!
!
d x
x′
=
dt x′
µ(1 − x2 )x′ − xと書けるので，x1 = x, x2 = x′ と置くと!
x1x′ = F (t, x),
x=
,
x2

!
F (t, x) =47

x2µ(1 − x21 )x2 − x1と書ける．
以下は，n 個の天体の運動を表した方程式である．
nX
d 2 ximj=−
G(xi − xj ),
2
3
dt
|x
−
x
|
i
j
j=1

i = 1, 2, · · · , nj̸=i

ここで，i 番目の天体の質量を mi ，位置を xi としている．G は万有引力定数である．この方程式も，微分の階数を減らすことで

x′ = F (t, x)の形に書くことができる．この方程式は，万有引力の法則とニュートンの運動法則から導くことができる．
 
x1
例えば，平面上を動く 2 つの天体の運動は，一つ目の天体の位置を x1 =
，二つ目
y1
 
x2
の天体の位置を x2 =，さらに ui = x′i , vi = yi′ (i = 1, 2) と置くと
y2


 
u
1
x1


 
v1


 y1 


 


u2
 x2 


 


v2
 y2 


 


 


m2 (x1 − x2 )
 
−G



u
2
2
3/2
1

{(x1 − x2 ) + (y1 − y2 ) } 
,
(t,
x)
F
=
x′ = F (t, x),
x=


 


m
(y
−
y
)
 
2
1
2
−G

 v1 

2 + (y − y )2 }3/2 
 
{(x
−
x
)


1
2
1
2
 


 
m
(x
−
x
)


1
2
1
 
−G

u2 
2
2
3/2

{(x1 − x2 ) + (y1 − y2 ) } 
 


 


m1 (y2 − y1 )−G
v2{(x1 − x2 )2 + (y1 − y2 )2 }3/2と書ける．
以上で見たように，多くの微分方程式の初期値問題は(
y ′ = f (x, y),y(x0 ) = y0の形に書くことができる．そこで，以降はこの形の微分方程式の数値解法について考えていくことにする．

48

4.2

ランダウの記号

今後よく使うので，ランダウの記号について説明しておこう．|x| が十分小さいときに，ある定数 M > 0 が存在して |f (x)| ≤ M |g(x)| が成り立つことをf (x) = O(g(x))と書くことにする．この O(·) をランダウの記号と言う．例えば

cos x = O(1)sin x = O(x)ex − 1 − x = O(x2 )xm = O(xn ),m≥n≥01= 1 − x + x2 + O(x3 )1+xなどが成り立つ．簡単に言うと，O(xs ) は，x の s 乗以上のオーダーで減少する項を表し
ている．

4.3

オイラー法

微分方程式

(
y ′ = f (x, y),y(x0 ) = y0の数値解法を考えるため，独立変数 x を分割幅 h で等間隔に離散化したものとして

x0 < x 1 < x 2 < · · · < x N ,xn+1 − xn = h (n = 0, 1, · · · , N − 1)を考え，y(xn ) の近似を yn で表すものとする．
いま，解 y(x) が十分に滑らかだとすると，テイラー展開により

y(xn+1 ) = y(xn + h) = y(xn ) + hy ′ (xn ) + O(h2 )が成り立つ．
ここで y ′ = f (x, y) であるので，上式はy(xn+1 ) = y(xn ) + hf (xn , y(xn )) + O(h2 )と変形できる．そこで，近似スキームとして

yn+1 = yn + hf (xn , yn )49

が考えられる．具体的には，(

y0 = y(x0 ),yn+1 = yn + hf (xn , yn )により，y(x1 ), y(x2 ), · · · , y(xN ) の近似である y1 , y2 , · · · , yN を求める．これをオイラー法という．
オイラー法の精度を調べるため，yn+1 と y(xn+1 ) がどの程度近くなるか調べてみることにしよう．
(
y(xn+1 ) = y(xn ) + hf (xn , y(xn )) + O(h2 ),yn+1 = yn + hf (xn , yn )より，辺々引くと


y(xn+1 ) − yn+1 = y(xn ) − yn + h f (xn , y(xn )) − f (xn , yn ) + O(h2 )となる．すなわち，仮に yn = y(xn ) となったとすると，yn+1 と y(xn+1 ) は h の 1 次の項まで一致する．そのため，オイラー法は 1 次精度であると言われる．
実際には，yn と y(xn ) の間にも誤差があり，さらに誤差は累積していくので，n = 1 から
n = N まで計算して，トータルでどの程度の誤差が発生するかという点については，もう少し詳しく分析する必要がある．これは後ほどまた検証する．

4.4

修正オイラー法

オイラー法の考え方は，y ′ (xn ) = f (xn , y(xn )) を用いて，y(xn+1 ) をテイラー展開の O(h)の項まで近似するというものであった．
そこで，さらに精度を上げるため，y ′ (xn ) = f (xn , y(xn )) だけではなく y ′ (xn+1 ) = f (xn+1 , y(xn+1 ))も用いることで，今度は y(xn+1 ) を O(h2 ) の項まで近似してみよう．
そのためには，仮に yn = y(xn ) となったとしたときに
2y(xn+1 ) = y(xn ) + hy ′ (xn ) + h y ′′ (xn ) + O(h3 ),2 yn+1 = yn + h αf (xn , y(xn )) + βf (xn+1 , y(xn+1 ))が O(h2 ) の項まで一致するように定数 α, β の値を決めたい．ここで，yn + h αf (xn , y(xn )) + βf (xn+1 , y(xn+1 )) = yn + h αy ′ (xn ) + βy ′ (xn+1 )

50



′
′
= yn + h αy (xn ) + βy (xn + h)
= yn + h αy ′ (xn ) + βy ′ (xn ) + βhy ′′ (xn ) + O(h2 ))1と取れば目標が達成されることがわかる．よって y ′ (xn ) と y ′ (xn+1 ) を
2
それぞれ yn ，yn+1 で置き換えたh
yn+1 = yn +f (xn , yn ) + f (xn+1 , yn+1 )2より，α = β =なる数値計算スキームが考えられるが，この形では右辺に yn+1 が入っているので使いにくい．そこで，右辺の yn+1 をオイラー法を用いて yn + hf (xn , yn ) で近似し，
hyn+1 = yn +f (xn , yn ) + f (xn+1 , yn + hf (xn , yn ))2というスキームを考える．このままでは括弧が多くて見にくいので書き直し，初期条件も加えると
y0 = y(x0 ),





k1 = hf (xn , yn ),k2 = hf (xn + h, yn + k1 ),





1
y
n+1 = yn + (k1 + k2 )2という形になる（ここで，k1 , k2 は n ごとに計算する）．この数値計算スキームを修正オ
イラー法（もしくはホイン法）という．
修正オイラー法を導出する際には，右辺に出てきた y(xn+1 ) を yn+1 で置き換え，さらにそれをオイラー法で近似したので，それでもちゃんと O(h2 ) の精度が保たれているか調べる必要がある．そのため，ここで二変数関数についての連鎖公式について復習しておこう．
いま，f (x, y) を C 2 級関数とし，x および y を変数 t の関数とすると

df (x, y)∂f dx ∂f dy
=
·
+
·
dt
∂x dt
∂y dt
が成り立つ．特に x = t のときにはdf (x, y)
∂f
∂f dy
=
+
·
dx
∂x ∂y dxとなる．
それでは，修正オイラー法の精度を調べよう．仮に yn = y(xn ) となったとし，f (x, y) が十分に滑らかだとすると，2 変数関数のテイラー展開を用いてh
f (xn , yn ) + f (xn+1 , yn + hf (xn , yn ))yn+1 = yn +
2

51


h
f (xn , yn ) + f (xn + h, yn + hy ′ (xn ))
2

h
∂
∂
′
2
= yn +f (xn , yn ) + f (xn , yn ) + h f (xn , yn ) + hy (xn ) f (xn , yn ) + O(h )
2
∂x
∂y

2
h ∂
∂
= yn + hf (xn , yn ) +f (xn , yn ) + y ′ (xn ) f (xn , yn ) + O(h3 )2 ∂x∂y= yn +

ここで，上で確認した連鎖公式を用いると

h2 d· f (xn , yn ) + O(h3 )2 dx
h2 d ′= yn + hf (xn , yn ) +· y (xn ) + O(h3 )2 dx
h2= yn + hy ′ (xn ) + y ′′ (xn ) + O(h3 )2yn+1 = yn + hf (xn , yn ) +が成り立つ．一方，h2 ′′y (xn ) + O(h3 )2なので，仮に yn と y(xn ) が一致したときには，yn+1 と y(xn+1 ) は h の 2 次の項まで一致する．以上により，修正オイラー法が 2 次精度であることがわかった．
y(xn+1 ) = y(xn + h) = y(xn ) + hy ′ (xn ) +4.5

4 次のルンゲ・クッタ法

実際の応用上は，オイラー法や修正オイラー法は精度の点からあまり使われない．よく使われるのは以下の 4 次のルンゲ・クッタ法である．

y0 = y(x0 ),





k1 = hf (xn , yn ),





k1
h



k2 = hf (xn + 2 , yn + 2 ),hk2k3 = hf (xn + , yn + ),



2
2




k4 = hf (xn + h, yn + k3 ),





1

yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )6
ここで，k1 , k2 , k3 , k4 は n ごとに計算する．4 次のルンゲ・クッタ法はその名称の通り 4 次
精度であることが知られている．4 次精度であることは，オイラー法や修正オイラー法と同様に示すことができるが，かなりの式変形を必要とするので容易ではない．ただ，腕力のある人は挑戦してみてもいいと思う．
数値計算スキームとしては 5 次以上のスキームも考えられるのだが，精度が 5 次以上になると手順が途端に複雑になるので，通常はあまり用いられない．

52

4.6

リプシッツ条件

次節以降の解析ではリプシッツ条件という概念を用いるので，ここで説明しておく．
関数 f (x) が領域 D でリプシッツ条件を満たすとは，ある定数 L > 0 が存在して，任意のx1 , x2 ∈ D について|f (x1 ) − f (x2 )| ≤ L|x1 − x2 |が成り立つことを言う．f (x) や x がベクトル値を取る場合には，絶対値のかわりに適当なノルムを用いる．f (x) が D でリプシッツ条件を満たすとき，f (x) は D でリプシッツ連
続であるという．簡単に言うと，リプシッツ条件というのは，D 内の任意の 2 点間の f (x)の平均変化率の絶対値が L 以下になるということである．
例えば，f (x) = x2 は D = [−1, 1] でリプシッツ条件を満たす．実際，L = 2 と取れば

|f (x1 ) − f (x2 )| = |x21 − x22 | = |x1 + x2 | |x1 − x2 | ≤ L|x1 − x2 |が成り立つ．一般に，f (x) が微分可能で，|f ′ (x)| が D で有界，つまり |f ′ (x)| ≤ M とな
るような M が存在するならばリプシッツ条件を満たす．実際
Z x2
Z x2
′
|f (x1 ) − f (x2 )| =f (x)dx ≤M dx = M |x1 − x2 |x1x1が成り立つ．もしくは，平均値の定理より，x1 と x2 の間の数 c が存在して|f (x1 ) − f (x2 )| = |f ′ (c)(x1 − x2 )| ≤ M |x1 − x2 |が成り立つ．
リプシッツ連続という概念は微分可能と連続の中間の概念と考えることができる．つまり，微分不可能だがリプシッツ条件を満たす関数や，連続だがリプシッツ条件を満たさない関数が存在する．例えば，f (x) = |x| は (−1, 1) で微分可能ではないが，リプシッツ条
件は満たす．実際，三角不等式から|x1 | ≤ |x1 − x2 | + |x2 |,が成り立つので，|x2 | ≤ |x2 − x1 | + |x1 ||f (x1 ) − f (x2 )| = ||x1 | − |x2 || = max |x1 | − |x2 |, |x2 | − |x1 | ≤ |x1 − x2 |pが成り立つ．一方，f (x) = |x| は [−1, 1] で連続だがリプシッツ条件は満たさない．実際，どんなに大きく L > 0 を取っても，0 < x1 < min(1/L2 , 1), x2 = 0 と取れば√
x1
x1
|f (x1 ) − f (x2 )| = x1 = √ > p= Lx1 = L|x1 − x2 |x1
1/L2となるので，リプシッツ条件を満たさないことがわかる．

53

4.7

オイラー法の大域誤差

前項までで，オイラー法と修正オイラー法の誤差について考察したが，その誤差は，y(xn )と yn が一致したときに次のステップで生じる誤差であった．実際には，誤差は累積するものであるし，また，ある一定区間を分割して数値スキームを適用するケースを考えると，分割の数が増えるに従い，各ステップの誤差は小さくなっても，ステップの数は増えるので，総合的にみて誤差が小さくなるかどうかはすぐにはわからない．そこで，本節ではオイラー法について，総合的に誤差がどうなるかを考えてみることにする．なお，既にみたような 1 ステップで生じる誤差を局所誤差，トータルとして生じる誤差を大域誤差と
いう．
微分方程式

(
y ′ = f (x, y),y(a) = y0

を区間 [a, b] 上で解く．そのため，分割数を N とし，

xn = a + nh,と分点を取って

(

h=

b−a
N

y0 = y(x0 ),yn+1 = yn + hf (xn , yn )により近似解を計算する．
誤差解析を行う上で，解 y(x) および f (x, y) について，ある程度の仮定を置く必要がある．
まず，解 y(x) は [a, b] で C 2 級であるものとし，

M = max ∥y ′′ (x)∥a≤x≤bとする．ここで ∥ · ∥ は任意のノルムである．f (x, y) や y がベクトルではなくスカラー関数である場合には，単なる絶対値を考えればよい．
次に，全ての a ≤ x ≤ b について y(x) ∈ D となるような適当な領域 D を取り，f (x, y) はa ≤ x ≤ b, y ∈ D で，y に関するリプシッツ条件を満たすものとする．すなわち，定数
L > 0 が存在して，すべての a ≤ x ≤ b および z1 , z2 ∈ D について∥f (x, z1 ) − f (x, z2 )∥ ≤ L∥z1 − z2 ∥が成り立つものとする．
さて，解 y(x) は [a, b] で C 2 級であるから

y(xn+1 ) = y(xn ) + hy ′ (xn ) +

h2
h2 ′′
y (ξ) = y(xn ) + hf (xn , y(xn )) + y ′′ (ξ)
2254が成り立つ．ここで h2 の項はラグランジュの剰余項であり，ξ は xn と xn+1 の間の数である．一方で，オイラー法のスキームは

yn+1 = yn + hf (xn , yn )なので，辺々引くと

 h2
yn+1 − y(xn+1 ) = yn − y(xn ) + h f (xn , yn ) − f (xn , y(xn )) − y ′′ (ξ)
2となる．これより∥yn+1 − y(xn+1 )∥ ≤ ∥yn − y(xn )∥ + h f (xn , yn ) − f (xn , y(xn )) +M h22が成り立つ．ここで，yn ∈ D を仮定すると，リプシッツ条件から∥yn+1 − y(xn+1 )∥ ≤ ∥yn − y(xn )∥ + Lh∥yn − y(xn )∥ += (1 + Lh)∥yn − y(xn )∥ +

M h2
2

M h2
2

が成り立つ．
一般的に，数列 {en } と実数 A, B （ただし A > 1 とする）がen+1 ≤ Aen + B(n ≥ 0)を満たすとき，An − 1B
A−1が成り立つ．証明は簡単で，n = 0 のときには明らかに成り立ち，n での成立を仮定すると
An − 1
An+1 − 1nen+1 ≤ Aen + B ≤ A A e0 +B + B = An+1 e0 +
B
A−1
A−1
en ≤ An e0 +より，n + 1 での成立も言えることから帰納的に示される．
すべての n について yn ∈ D を仮定するとき，上の不等式を用いると∥yn − y(xn )∥ ≤ (1 + Lh)n ∥y0 − y(x0 )∥ +=((1 + Lh)n − 1)M h
2L(1 + Lh)n − 1 M h2·
Lh2が成り立つ．さらに 1 + Lh ≤ eLh を用いると∥yn − y(xn )∥ ≤(eLN h − 1)M h(eL(b−a) − 1)M h(enLh − 1)M h
≤
=
2L
2L
2L

55

が成り立つことがわかる．ここで，a, b, L, M は全て定数なので，オイラー法の大域誤差
は O(h) となることがわかる．
上の解析では yn ∈ D を仮定していたが，最初に D を解 y(x) の動く範囲より大きく取っておけば，h を小さく取ることで yn がすべて D に入るようにできるので，特に問題は生
じない．

4.8

オイラー法のプログラム

まずは，答えの分かっている問題を解いて精度を確認してみよう．N が 10 倍になるごとに誤差が 10 分の一になっていればよい．

#include <stdio.h>#include <math.h>#include "Linear.h"// y’=y, y(0)=1 を解く。解は y=e^x
double f(double x, double y){return y;}// 解double sol(double x){
return exp(x);}int main(){
printf("***** オイラー法 *****\n");double a = 0.0, b = 1.0; // 計算範囲
for(int N=10; N<=10000; N*=10){
Vector x(N);Vector y(N);
次ページへ続く…

56

double h = (b-a)/N; // 区間幅// 分点を作成
for(int n=0; n<=N; n++){x[n] = a + n*h;}// オイラー法
y[0] = 1; // 初期条件
for(int n=0; n<N; n++){y[n+1] = y[n] + h*f(x[n],y[n]);}// 誤差を調査
double maxerr = 0;
for(int n=0; n<=N; n++){double err = fabs( y[n] - sol(x[n]) );if(err>maxerr){ maxerr = err;}}
printf("N=%d のとき、誤差は%.15f\n", N, maxerr);}
return 0;}

4.9

ルンゲ・クッタ法のプログラム

次に，同じ問題を用いて 4 次のルンゲ・クッタ法の精度を確認してみよう．N が 10 倍になるごとに誤差が 1 万分の一になっていればよい．

57

#include <stdio.h>#include <math.h>#include "Linear.h"// y’=y, y(0)=1 を解く。解は y=e^x
double f(double x, double y){return y;}// 解double sol(double x){
return exp(x);}int main(){
printf("***** 4 次のルンゲ・クッタ法 *****\n");double a = 0.0, b = 1.0; // 計算範囲
for(int N=10; N<=10000; N*=10){
Vector x(N);Vector y(N);double h = (b-a)/N; // 区間幅// 分点を作成
for(int n=0; n<=N; n++){x[n] = a + n*h;}// ルンゲ・クッタ法
y[0] = 1; // 初期条件
for(int n=0; n<N; n++){double k1 = h * f(x[n],y[n]);
double k2 = h * f(x[n]+h/2,y[n]+k1/2);
double k3 = h * f(x[n]+h/2,y[n]+k2/2);
double k4 = h * f(x[n]+h ,y[n]+k3 );y[n+1] = y[n] + (k1 + 2*k2 + 2*k3 + k4)/6;}
次ページへ続く…

58

// 誤差を調査
double maxerr = 0;
for(int n=0; n<=N; n++){double err = fabs( y[n] - sol(x[n]) );if(err>maxerr){ maxerr = err;}}
printf("N=%d のとき、誤差は%.15f\n", N, maxerr);}
return 0;}
次に，ロトカ・ボルテラ方程式を解いて，結果をファイルに出力してみよう．ファイルへの出力には printf の代わりに fprintf を用いる．出力した結果は，Excel や gnuplot などのアプリケーションを用いてグラフ化することができる．

#include <stdio.h>#include "Linear.h"// x_1’=(a-b*x_2)*x_1 を解く。ここで、x1, x2 は時間 t の関数// x_2’=(-c+d*x_1)*x_2
Vector f(double t, Vector x){
Vector y(2);y[1] = (1.0 - 1.0*x[2])*x[1];y[2] = (-1.0 + 0.5*x[1])*x[2];return y;}
int main(){
printf("***** ロトカボルテラ方程式 *****\n");double a = 0.0, b = 30.0; // 計算範囲
int N=1000; // 分割数
次ページへ続く…

59

int d = 2; // 従属変数の次元
Matrix S(d, N);
Vector t(N);Vector x(d), k1(d), k2(d), k3(d), k4(d);double h = (b-a)/N; // 区間幅x[1] = 0.5; x[2] = 0.5; // 初期条件
S.SetColVector(0, x);for(int n=0; n<=N; n++){t[n] = a + n*h;}for(int n=0; n<N; n++){k1 = h * f(t[n],x);k2 = h * f(t[n] + h/2, x + k1/2);
k3 = h * f(t[n] + h/2, x + k2/2);
k4 = h * f(t[n] + h , x + k3 );x += (k1 + 2*k2 + 2*k3 + k4)/6;
S.SetColVector(n+1, x);}
FILE *fp;
fp = fopen("l-v.dat", "w");for(int n=0; n<=N; n++){fprintf(fp, "%f", t[n]);for(int i=1; i<=d; i++){fprintf(fp, " %f", S[i][n]);}fprintf(fp, "\n");}fclose(fp);return 0;}
このプログラムを実行すると l-v.dat というファイルができる．gnuplot でplot "l-v.dat" using 1:2 with lines, "l-v.dat" using 1:3 with linesと入力するとグラフを表示することができる．

60

4.10

補足

オイラー法，ルンゲ・クッタ法などは，正確には，陽的オイラー法，陽的ルンゲ・クッタ法と呼ばれることもある．常微分方程式の数値解法には，他にも，陰的オイラー法，陰的ルンゲ・クッタ法，ピカールの逐次近似法，予測子・修正子法など，様々な方法があり，一長一短があるが，おそらく最もよく用いられているのは，今回説明した 4 次の陽的ルン
ゲ・クッタ法である．
本稿では，常微分方程式の初期値問題の数値解法について説明したが，境界値問題，例えば
y(a) = ya ,y(b) = yb ,
 y ′′ = f (x, y, y ′ )のような形の方程式も重要である．境界値問題の数値解法としては，差分法などが用いられることが多いが，時間の関係で今回は省略した．

4.11

課題

課題１ 解のわかっている微分方程式を自由に選び，オイラー法，修正オイラー法，4 次のルンゲ・クッタ法で計算を行い，誤差が刻み幅に対してどのように減少していくか調べなさい
課題２ 2 階の微分方程式を自由に選び，4 次のルンゲ・クッタ法を用いて解を計算し，グラフを描きなさい．横軸に独立変数，縦軸に y と y ′ を取ったグラフと，横軸に y ，縦軸にy ′ を取ったグラフ（相図）を作成すること．

61

5

偏微分方程式の差分解法

本章では，偏微分方程式の数値解法について説明する．偏微分方程式とは，関数の偏微分についての関係を表した方程式であり，最も基本的なものとして，熱方程式が知られている．熱方程式は，その名前からもわかる通り，元々は物理の方程式であるが，自然現象や社会現象など，様々な所に現れる極めて重要な方程式である．また，熱方程式に対する数値解法は，他の様々な偏微分方程式の数値解法へ応用が可能であるので，熱方程式の数値解法を知ることは大いに意味がある．そこで，まずは熱方程式とその数値解法を見ていくことにしよう．

5.1

熱方程式

以下の偏微分方程式が熱方程式である（導出については後ほど説明する）．

∂u
∂ 2u



=
λ
(t > 0, a < x < b),


∂t
∂x2


u(0, x) = u0 (x)(a ≤ x ≤ b), u(t, a) = α(t)(t > 0) u(t, b) = β(t)(t > 0)
ここで，x は空間方向の位置を表し，t は時間を表す変数である．u(t, x) が求めたい未知関数で，u0 (x), α(t), β(t) は与えられた関数である．また，λ は正の定数である．
熱方程式は，区間 [a, b] で表される針金の熱分布の時間変化を表す方程式とみることができる．ここで，条件 u(0, x) = u0 (x) は，時刻 t = 0 における温度分布を表しているので，初期条件と呼ばれ，条件 u(t, a) = α(t) と u(t, b) = β(t) は，端点 x = a および x = b での温度を α(t) および β(t) に固定することを意味しているので，境界条件と呼ばれる．特にこのような，関数値を指定する条件をディリクレ境界条件という．場合によっては，端点における微分の値を指定することもあり，そのような条件をノイマン境界条件と言うが，ここでは簡単のためディリクレ境界条件のみを考える．
ここで，熱方程式の導出について説明しておこう．一本の針金を考え，時刻 t，位置 x ∈ [a, b]の温度を u(t, x) とする．熱は時間とともに針金上を移動するが，その際の熱の流量は，熱分布の傾き（温度勾配）の大きさに比例することが知られている．つまり，温度変化が急な部分ほど熱が多く移動するということである．具体的には，各点における単位時間かつ単位断面積あたりの熱流量 Q は∂u
Q = −K∂xと表される．ここで K は針金の熱伝導率である．ここで，マイナスの符号が付いているのは，熱は温度勾配と逆方向に流れるからである．

62ax x+∆xb

図 2: 針金上の熱分布
さて，針金上の位置 x および x + ∆x における断面で囲まれる部分を V とすると，単位時間あたりに V に流入する熱量は，断面積を S として


∂
∂
∂
∂
−SK u(t, x) + SK u(t, x + ∆x) = SKu(t, x + ∆x) −u(t, x)∂x∂x∂x∂xとなる．さらに，V の体積は S∆x であるので，針金の比熱を c，密度を ρ とすると，単位時間あたりの V の温度変化は，平均すると




SK
∂
K
∂
∂
∂
u(t, x + ∆x) −u(t, x) =u(t, x + ∆x) −u(t, x)cρS∆x ∂x∂xcρ∆x ∂x∂xとなる．ここで，∆x → 0 とすると


1
∂
∂2
∂
u(t, x + ∆x) −u(t, x)→u(t, x)
∆x ∂x
∂x
∂x2
より，

∂u
K ∂ 2u=·∂tcρ ∂x2が成り立つことがわかる．λ =5.2

Kとおけば，最初に紹介した熱方程式の形になる．
cρ

差分

微分方程式の数値計算では，微分を関数の差で近似することが多い．微分を近似する関数の差を差分といい，差分を用いる数値計算法を差分法という．それでは，いくつかの差分を紹介しよう．

f (x) を C 2 級関数とすると，テイラー展開より
f (x + h) = f (x) + hf ′ (x) + O(h2 )63

が成り立つ．よって

f (x + h) − f (x)+ O(h)hf (x + h) − f (x)が成り立つ．すなわち，
で f ′ (x) を近似することができる．これを（一h階の）前進差分という．
f ′ (x) =同様に，

f (x − h) = f (x) − hf ′ (x) + O(h2 )より

f ′ (x) =も成り立つ．

f (x) − f (x − h)+ O(h)
hf (x) − f (x − h)を（一階の）後退差分という．
hまた，f (x) が C 3 級のとき，

h2 ′′f (x) + O(h3 )2h2f (x − h) = f (x) − hf ′ (x) + f ′′ (x) + O(h3 )2f (x + h) = f (x) + hf ′ (x) +より

f (x + h) − f (x − h) = 2hf ′ (x) + O(h3 )が成り立つので，

f ′ (x) =f (x + h) − f (x − h)+ O(h2 )2hf (x + h) − f (x − h)が成り立つ．
を（一階の）中心差分という．f (x) が滑らかな場合，中2h
心差分は前進差分や後退差分より精度が良いということがわかる．前進差分と後退差分
は，まとめて片側差分と呼ばれることもある．
次に，2 回微分の近似を考えてみよう．f (x) が C 4 級のとき，

h2 ′′
h3
f (x) + f ′′′ (x) + O(h4 )
2
6
2
3
h
h
f (x − h) = f (x) − hf ′ (x) + f ′′ (x) − f ′′′ (x) + O(h4 )26f (x + h) = f (x) + hf ′ (x) +が成り立つので，f (x + h) − 2f (x) + f (x − h)+ O(h2 )h2
f (x + h) − 2f (x) + f (x − h)を（二階の）中心差分という．
が成り立つ．
h2
f ′′ (x) =

64

5.3

陽解法

それでは，熱方程式の数値解法に入ろう．時間方向 t を刻み幅 ∆t で分割し，上付きの添え字を用いてtm = m∆tb−aとする（m 乗ではなく添え字なので注意）．また，空間方向 x を刻み幅 ∆x =で分N割し，
xn = a + n∆xとする．さらに，u(tm , xn ) の近似を umn とし，微分の近似として前進差分および中心差分∂uu(tm+1 , xn ) − u(tm , xn )um+1 − um
n
≈
≈ n
,
∂t
∆t
∆t
m
m
um
u(tm , xn+1 ) − 2u(tm , xn ) + u(tm , xn−1 )∂ 2u
n+1 − 2un + un−1
≈
≈
∂x2
∆x2
∆x2
を導入し（≈ は近似を意味する），熱方程式に代入すると
m
m
um
um+1
− um
n+1 − 2un + un−1
n
n
=λ
∆t
∆x2

となる（差分による近似と umn による近似が二重に入っているので，本当に近似になっているかどうかは厳密にはわからないが，今のところはよしとしておく）．上の式を見ると，
m
m
0
um+1
を um
n
n−1 , un , un+1 で書くことができるので，t = t から順に計算していき，熱方程式の近似解を求めることができる．
初期条件や境界条件も入れてまとめて書くと，r = λu0n = u0 (xn )um+1 = rum + (1 − 2r)um + rum
n
n+1
n
n−1
m+1
m+1

u0 = α(t
)



 m+1
m+1
uN = β(t
)

∆t
と置いて
∆x2
(0 ≤ n ≤ N ),(m ≥ 0, 0 < n < N ),(m ≥ 0),(m ≥ 0)となる．この計算スキームは，次の時間ステップの近似解を前の時間ステップの近似解を用いて書き下せる（これを陽に書き下せる，という）ので，差分法の中でも陽解法と呼ば
れている．

5.4

解の収束性

1
を満たしながら ∆t と ∆x を小さくしていけば，真の解に収束すること2が知られている．ここではこれを示そう．まず，陽解法は，r ≤mmem
n = u(t , xn ) − un65

とする．つまり，emn は真の解と数値解の差とする．このとき，解が十分に滑らかだとすm
ると，t = t , x = xn において∂uu(tm+1 , xn ) − u(tm , xn )=+ O(∆t),∂t∆t∂ 2uu(tm , xn+1 ) − 2u(tm , xn ) + u(tm , xn−1 )=+ O(∆x2 )22∂x∆xが成り立つので，u(tm+1 , xn ) − u(tm , xn )u(tm , xn+1 ) − 2u(tm , xn ) + u(tm , xn−1 )= O(∆t) + O(∆x2 )
−λ
∆t
∆x2
となる．一方で，陽解法のスキームよりmum − 2um
um+1− um
n + un−1nn=0− λ n+1∆t∆x2

なので，辺々引くとmem − 2em
em+1− em
n + en−1nn− λ n+1= O(∆t) + O(∆x2 )∆t∆x2となる．さらに両辺に ∆t をかけて整理すると



m
m
m
2
em+1
=
re
+
(1
−
2r)e
+
re
+
∆t
O(∆t)
+
O(∆x
)
n
n−1
n
n+1

1
のとき，n について両辺の絶対値の最大値を取ると
2


m
m
m
2
max |em+1
|
≤
r
max
|e
|
+
(1
−
2r)
max
|e
|
+
r
max
|e
|
+
∆t
O(∆t)
+
O(∆x
)
n
n
n
n
0<n<N
0<n<N
0<n<N
0<n<N


m
2
= (r + (1 − 2r) + r) max |en | + ∆t O(∆t) + O(∆x )
0<n<N


m
= max |en | + ∆t O(∆t) + O(∆x2 )となる．ここで r ≤0<n<Nが成り立つ．ここで，時間方向の計算範囲を 0 ≤ t ≤ T とし，∆t =m ≤ M についてTとすると，0 <
M





0
2
2
max |em
|
≤
max
|e
|
+
m∆t
O(∆t)
+
O(∆x
)
=
m∆t
O(∆t)
+
O(∆x
)
n
n
0<n<N
0<n<N




2
2
≤ M ∆t O(∆t) + O(∆x ) = T O(∆t) + O(∆x )1のとき，∆t と ∆x を小さくしていけば，陽解法による数値解は2真の解に収束することがわかった．

となる．これより，r ≤1のときには収束性は必ずしも言えない．実際，このときには陽解法の数値解
2が収束しないような解の例を作ることができる．
逆に，r >

66

5.5

安定性

実際の数値計算においては，計算過程で丸め誤差などの誤差が混入する．よって，安定して数値計算を行うには，混入した誤差が増幅されないようにしなければならない．
m
m+1
いま，um
に
n に εn の誤差が混入したとしよう．このとき，この誤差の影響によって un
m+1
加わる誤差 εn を求めてみると
m
m
m
m
m
um+1
+ εm+1= r(umnnn−1 + εn−1 ) + (1 − 2r)(un + εn ) + r(un+1 + εn+1 )

 

m
m
m
m
m
= rum
+
(1
−
2r)u
+
ru
+
rε
+
(1
−
2r)ε
+
rε
n−1
n
n+1
n−1
n
n+1

より
m
m
εm+1
= rεm
n
n−1 + (1 − 2r)εn + rεn+1と表されることがわかる．r ≤

1
のときには
2

m
m
m
max |εm+1| ≤ r max |εmnn | + (1 − 2r) max |εn | + r max |εn | = max |εn |

0<n<N

0<n<N

0<n<N

0<n<N

0<n<N

より，時間ステップを経ても混入した誤差は増大しないということがわかる．
一方で r >

1nmのとき，少々天下り的だが，umn に εn = η(−1) の誤差が混入したとすると，
2
m
m
εm+1
= rεm
n
n−1 + (1 − 2r)εn + rεn+1= rη(−1)n−1 + (1 − 2r)η(−1)n + rη(−1)n+1= (1 − 4r)η(−1)n1のとき
2|1 − 4r| > 1 となるので，誤差は時間ステップを経るごとに増幅され，振幅が増大していくことになる．
より，um+1には εm+1= (1 − 4r)η(−1)n の誤差が混入することになる．r >nn

数値計算においては，通常はありとあらゆる形の誤差が混入するので，いま述べたような1nεm（ごくわずかかもしれないが）混入してしまう．r > のn = η(−1) という形の誤差も，2ときには，その誤差は指数関数的に増大するので，数値解を安定して計算することが困難になってしまう．

1のときには混入した誤差が増加しない．このような，時2間ステップが進んでも誤差が増幅されないための条件を安定性条件という．陽解法の安定
1性条件は r ≤ である．
2
一方で，上で見たように，r ≤

67

5.6

陰解法

陽解法において安定性条件を満たすためには，時間刻みを空間刻みの 2 乗のオーダーに取らなければならない．もし，空間方向の構造を細かく見るために ∆x を細かく取ったとすると，∆t を非常に小さく取らねばならず，計算量という点で難がある．例えば，λ = 1, ∆x = 0.001のときには，陽解法の安定性条件を満たすために ∆t を 0.0000005 以下に取らねばならず，
時間ステップを進めるのに非常に多くの計算が必要になってしまう．そこで用いられるのが陰解法である．陰解法では，時刻 tm+1 において，差分化した関係式が成り立つように
um· を決める．具体的には
m+1+ um+1
um+1
um+1
− um
n−1
n+1 − 2un
n
n
=λ
2
∆t
∆x

と差分化することになる．証明には行列の固有値解析が必要なので詳しくは述べないが，陰解法は無条件安定，つまりどんな ∆t や ∆x に対しても安定に計算できることが知られている．
初期条件や境界条件を加味すると，陰解法は
u0n = u0 (xn )





um+1
= α(tm+1 )
0
m+1
m

−rum+1− rum+1
n−1 + (1 + 2r)un
n+1 = un
um+1
= β(tm+1 )N(0 ≤ n ≤ N ),(m ≥ 0),(m ≥ 0, 1 ≤ n ≤ N − 1),(m ≥ 0)となる．解は毎回，連立一次方程式を解いて求めることになる．ここで，端点 um+1と um+10
Nは境界条件から決まるのだが，この端点も未知数として連立一次方程式を立てると，プログラムが簡単になる．具体的には，N + 1 次の連立一次方程式
  m+1  

1
0
0
···
0
0
0
u0
α(tm+1 )  m+1  m
−r 1 + 2r −r · · ·

 

0
0
0

 u1   u1   m+1  −r 1 + 2r · · ·000  u2   um
0

2


 

 ..

..
..
..
..
..   ..  = 
..
..
 .





.
.
.
.
.
.  .  .
  m+1   m
00· · · 1 + 2r −r
0  uN −2   uN −2 
0


 

0
 m+1   m

0
0
···
−r 1 + 2r −r uN −1   uN −1 um+1β(tm+1 )
0
0
0
···
0
0
1
N
m+1
を解くことで，umを計算する．連立一次方程式の係数行列は 3 重対角行列とな· から u·り，最初と最後の行以外は，対角成分が 1 + 2r，その両脇の成分が −r なる．

68

5.7

クランク・ニコルソン法

今まで述べた陽解法や陰解法では，時間に関する片側差分の誤差が O(∆t)，空間に関する
中心差分の誤差が O(∆x2 ) となっており，時間方向の精度があまり良くない．ところで，
時間方向の差分u(tm+1 , xn ) − u(tm , xn )
∆t
′ m
′ m+1
は，u (t , xn ) あるいは u (t, xn ) の近似として見れば片側差分となり，誤差は O(∆t)′ m+1/2
であるが，u (t, xn ) の近似と見れば中心差分と考えることができ，誤差は O(∆t2 ) となる．よって，t = tm+1/2 で方程式を離散化すれば，時間に関しても O(∆t2 ) の精度を実現できると考えられる．しかし，t = tm+1/2 は時間分割の分点ではないというのが問題である．
一般に，f (x) が C 2 級のとき，f (x + h) = f (x) + hf ′ (x) + O(h2 )f (x − h) = f (x) − hf ′ (x) + O(h2 )より，

f (x + h) + f (x − h)+ O(h2 )2となることから，t = tm+1/2 における熱方程式の左辺を t = tm のときと t = tm+1 のときの平均で近似すると


∂
1 ∂
∂
m+1/2
m+1
m
u(t
, x) =u(t, x) + u(t , x) + O(∆t2 )
∂t
2 ∂t
∂t
f (x) =

が成り立つ．よって，右辺の 2 項をそれぞれ 2 階の中心差分で近似することにより，t =tm+1/2 での離散化を実現することができる．具体的には


m
m
m+1
um
+ um+1− um
um+1
λ um+1n+1 − 2unn−1
n+1 − 2un + un−1
n
n
=
+
∆t
2
∆x2
∆x2
となる．r を用いて書き換えると
m+1
m
m
m
−rum+1− rum+1
n−1 + 2(1 + r)unn+1 = run−1 + 2(1 − r)un + run+1となる．これをクランク・ニコルソン法という．クランク・ニコルソン法は O(∆t2 )+O(∆x2 )の精度を持ち，さらに無条件安定であることが知られている．クランク・ニコルソン法も陰解法の一種なので，解は毎回，連立一次方程式を解いて求めることになる．具体的には，

69  m+1 u0  m+1 −r 2(1 + r) −r


···
0
0
0

  u1 −r 2(1 + r) · · ·000  um+1
0
 2 

 ..

..
..
.
.
.
.
...
..
..
..   .. 
 .

.
.




  m+1 
00· · · 2(1 + r) −r
0  uN −2 
0



0
 um+1 
0
0
·
·
·
−r
2(1
+
r)
−r

  N −1 
0
0
0
···
0
0
1
um+1
N


α(tm+1 )


m
m
 rum

+
2(1
−
r)u
+
ru
0
1
2




m
m
m
 ru1 + 2(1 − r)u2 + ru3





.
.
=

.


 m
m 
+
ru
ruN −3 + 2(1 − r)um
N −2
N −1  rum + 2(1 − r)um + rum 
N −1
N 
 N −2
m+1
β(t
)
1

0

0

···

0

0

0

を解けばよい．

5.8

陽解法のプログラム

陽解法のプログラムを以下に示す．熱が変化していく様子をグラフにするには，Excel や
gnuplot を用いて 3 次元グラフを表示すればよい．全ての空間ステップ，時間ステップでの u の値を出力するとデータが多くなり過ぎてうまくグラフが描けないので，MM や NN
でファイルに出力する間隔を指定している．初期条件 u0 は u0 (x) = min(x, 1 − x) という山型の関数を与えている．
全ての時間刻みと空間刻みについて数値解を計算するとデータ量が膨大になるので，時間方向は 1 ステップごとに前回の結果を上書きするようにしている．

70

#include <stdio.h>#include "Linear.h"#define M 1000 // 時間分割
#define MM 50 // 時間のファイル出力間隔#define N 20 // 空間分割
#define NN 1 // 空間のファイル出力間隔// 境界条件
double alpha(double t){return 0;}
double beta(double t){return 0;}// 初期条件
double u0(double x){
if(x<0.5){
return x;} else {return 1-x;}}
int main(){double lambda = 1;
double a = 0; // 左端
double b = 1; // 右端
double T = 0.5; // 時間の計算範囲
double dt = T/M;
double dx = (b-a)/N;
次ページへ続く…

71

Vector x(N);Vector u(N);Vector u2(N);printf("***** 熱方程式（陽解法） *****\n");double r = lambda*dt/dx/dx;
printf("r = %.15f\n", r);// 初期条件および x_n を設定
for(int n=0; n<=N; n++){x[n] = a + n*dx;
u[n] = u0(x[n]);}FILE *fp;
fp = fopen("heat_explicit.dat","w");for(int m=0; m<=M; m++){double t = m*dt;
if(m%MM==0){
for(int n=0; n<=N; n+=NN){
fprintf(fp, "%f %f %f\n",t, x[n], u[n]);}fprintf(fp, "\n");}t = (m+1)*dt;
u2[0] = alpha(t);u2[N] = beta(t);for(int n=1; n<=N-1; n++){u2[n] = r*u[n+1] + (1-2*r)*u[n] + r*u[n-1];}
u = u2;}
fclose(fp);return 0;}72

グラフを表示するには，gnuplot でsplot "heat_explicit.dat" with linesと入力する．

Excel でグラフを表示することもできる．この場合は，プログラムの一部を以下のように書き換え，出力された csv ファイルを Excel で読み込み，グラフの種類として「ワイヤーフレーム 3D 等高線」を描けばよい．
--------- ここまで同じ --------FILE *fp;
fp = fopen("heat_explicit.csv","w");for(int n=0; n<=N; n+=NN){
fprintf(fp, ",%f",x[n]);}fprintf(fp, "\n");for(int m=0; m<=M; m++){double t = m*dt;
if(m%MM==0){
fprintf(fp, "%f",t);for(int n=0; n<=N; n+=NN){
fprintf(fp, ",%f",u[n]);}fprintf(fp, "\n");}t = (m+1)*dt;--------- これ以降同じ --------陽解法では，安定性条件 r ≤ 0.5 が満たされないときには計算が上手く行かないことも確
認しておこう．

5.9

陰解法のプログラム

陰解法の場合は，毎時間ステップごとに連立一次方程式を解く必要があるが，ここではガウスの消去法を用いよう．プログラムは以下のようになる．

73

#include <stdio.h>#include "Linear.h"#define M 100 // 時間分割
#define MM 5 // 時間のファイル出力間隔#define N 20 // 空間分割
#define NN 1 // 空間のファイル出力間隔// 境界条件
double alpha(double t){return 0;}
double beta(double t){return 0;}// 初期条件
double u0(double x){
if(x<0.5){
return x;} else {return 1-x;}}
int main(){double lambda = 1;
double a = 0; // 左端
double b = 1; // 右端
double T = 0.5; // 時間の計算範囲
double dt = T/M;
double dx = (b-a)/N;
次ページへ続く…

74

Matrix A(N+1, N+1);
Vector x(N);Vector u(N+1);
printf("***** 熱方程式（陰解法） *****\n");double r = lambda*dt/dx/dx;
printf("r = %.15f\n", r);// 行列を作成
for(int i=1; i<=N+1; i++){for(int j=1; j<=N+1; j++){if(i==j){A[i][j] = 1+2*r;} else if(i-j==1 || i-j==-1){
A[i][j] = -r;} else {A[i][j] = 0;}}}
A[1][1]= 1; A[1][2]= 0;
A[N+1][N] = 0; A[N+1][N+1] = 1;// 初期条件および x_n を設定
for(int n=0; n<=N; n++){x[n] = a + n*dx;
u[n] = u0(x[n]);}次ページへ続く…

75

FILE *fp;
fp = fopen("heat_implicit.dat","w");for(int m=0; m<=M; m++){double t = m*dt;
if(m%MM==0){
for(int n=0; n<=N; n+=NN){
fprintf(fp, "%f %f %f\n",t, x[n], u[n]);}fprintf(fp, "\n");}t = (m+1)*dt;
u[0] = alpha(t);u[N] = beta(t);u = u.Shift(1); // 添え字を 1 始まりにするために 1 だけシフトする
u = A.Gauss(u);u = u.Shift(-1); // 添え字を 0 始まりに戻す}
fclose(fp);return 0;}
Excel でグラフを描く場合は，ファイルに出力する部分を陽解法のときと同様に変更すれ
ばよい．

5.10

差分解法の応用

熱方程式以外の偏微分方程式についても，差分解法を用いて数値解を求められることが多い．ここでは，ファイナンスの分野に現れる Black-Scholes 方程式
∂Y
σ2S 2 ∂ 2Y
∂Y

=−
·
+ rY
(0 ≤ t ≤ T ),− rS ·
2
∂t
2
∂S
∂S

Y (T, S) = F (S)の数値解法を考えてみよう．ここで，r は無リスク金利（正の定数），σ はボラティリティ―（株価の変動の大きさを表す正の定数）である．Black-Scholes 方程式は，満期 t = T
において，その時点での株価 S に応じた払い戻し F (S) が受けられるような金融商品の，時刻 t < T かつ株価が S のときの価格 Y (t, S) が満たす偏微分方程式である．

76

K を正の定数とし，F (S) = max(S − K, 0) とする場合，すなわち，満期における株価 S
が K を上回っていた場合には K と株価との差額の払い戻しが受けられ，株価 S が K を下回っていた場合には払い戻しが無いような金融商品をコール・オプションという．また，F (S) = max(K − S, 0) であるような金融商品をプット・オプションという．いずれの場
合でも，K を権利行使価格という．
株価が 2 倍になれば日々の値動きの大きさも 2 倍になる，というように，一般的に株価は比率で動くので，Black-Scholes 方程式は価格そのもので見るよりも，価格の対数を見る方が自然である．そこで，x = log S とおいて，Y を x の関数として求めよう．

∂Y
∂Y dx
1 ∂Y
=
·
= ·
∂S
∂x  dS
S ∂x
2
∂ Y
∂
1 ∂Y
1 ∂Y
1 ∂ 2 Y dx
1 ∂Y
1 ∂ 2Y
=
·
=
−
·
+
·
·
=
−
·
+
·
∂S 2
∂S S ∂x
S 2 ∂xS ∂x2 dSS 2 ∂xS 2 ∂x2
より，代入すると

∂Y
σ2 ∂ 2Yσ 2 ∂Y
=− ·
− r−
+ rY
∂t
2 ∂x2
2 ∂x

が得られる．これを解くために，時間刻み ∆t，価格刻み ∆x を決め，等間隔に取った離散点を (tm , xn ) とし，Y (tm , xn ) の近似を Ynm とする．
まずは熱方程式の場合と同様，陽解法を考えよう．差分による近似として

Y ≈ Ynm+1 ,∂YY m+1 − Ynm
≈ n
,
∂t
∆t
m+1
m+1
Yn+1
− Yn−1
∂Y
≈
,
∂x
2∆x
m+1
m+1
Yn+1
− 2Ynm+1 + Yn−1
∂ 2Y
≈
,
∂x2
∆x2
を導入し，Black-Scholes 方程式に代入すると
m+1
m+1− 2Ynm+1 + Yn−1
σ 2 Yn+1
Ynm+1 − Ynm
=− ·
∆t
2
∆x2


m+1
m+1
2
Yn+1
− Yn−1σ− r−+ rYnm+1
2
2∆x
m+1
m+1
となる．上の式を見ると，Ynm が Yn−1
, Ynm+1 , Yn+1
で書けるので，満期から時間を逆に計算して行くことで時刻 0 での値を求めることができる．時間を逆向きに計算するところが熱方程式との違いである．

Black-Scholes 方程式において x の範囲は全ての実数となるが，数値計算の際には，適当な区間に限定して計算を行う．計算する金融商品がコール・オプションやプット・オプショ

77

ンの場合は，権利行使価格を中心とし，適当な R > 1 を取って log K/R ≤ x ≤ log KR とすることが多い．さらに，詳しい説明は省くが，境界条件は，コール・オプションの場合

Ykm = er(T −t ) F (exk )mプット・オプションの場合

Ykm = e−r(T −t ) F (exk )mと取ればよいことがわかっている．ただし，k は境界の添え字である．
熱方程式の場合と同様に考えると，Black-Scholes 方程式の陽解法の安定性条件は

1∆t≤ 2 + O(∆t)
2∆xσであることがわかる．

Black-Scholes 方程式においても熱方程式と同様，陰解法やクランク・ニコルソン法を考えることができる．陰解法が無条件安定になること，クランク・ニコルソン法の誤差が時間刻みに対して二次精度である点も熱方程式の場合と同じである．

5.11

Black-Scholes 方程式の陽解法のプログラム

Black-Scholes 方程式を陽解法で解くプログラムを以下に示す．
#include <stdio.h>#include <math.h>#include "Linear.h"#define M 1000 // 時間分割
#define MM 50 // 時間のファイル出力間隔#define N 20// 空間分割#define NN 1// 空間のファイル出力間隔
次ページへ続く…

78

// ペイオフ関数
double F(double S){
double K = 100;if(S>K){
return S-K;
} else {return 0;}}
int main(){// パラメータの設定
double SL= 50;// 下限価格
double SH= 200;// 上限価格
double r= 0.1;// 無リスク金利
double sig = 0.2;// ボラティリティ
double T= 2.0;// 満期
double dt= T / M;// 時間分割幅
double ds= log(SH/SL)/N; // log 価格分割幅
printf("安定性条件チェック: dt=%f, ds=%f, sigma^2*dt/ds^2=%f\n",dt, ds, sig*sig*dt/ds/ds);// 原資産価格
Vector S(N);
for(int n=0; n<=N; n++){S[n] = SL*exp(n*ds);}
Vector Y(N);
Vector Y2(N);
次ページへ続く…

79

FILE *fp;
fp = fopen("b-s.dat","w");// 満期のデータを出力
for(int n=0; n<=N; n+=NN){
Y[n] = F(S[n]);fprintf(fp,"%f %f %f\n", T, S[n], Y[n]);}fprintf(fp,"\n");// 係数は一定なので、あらかじめ求めておく
double a = dt*sig*sig/ds/ds/2 - dt*(r-sig*sig/2)/2/ds;
double b = 1 - dt*sig*sig/ds/ds - dt*r;
double c = dt*sig*sig/ds/ds/2 + dt*(r-sig*sig/2)/2/ds;
for(int m=M-1; m>=0; m--){double t = m*dt;
double er = exp(r*(T-t));Y2[0] = F(er*S[0]);Y2[N] = F(er*S[N]);for(int n=1; n<N; n++){Y2[n] = a*Y[n-1] + b*Y[n] + c*Y[n+1];}
Y = Y2;if(m%MM==0){
for(int n=0; n<=N; n+=NN){
fprintf(fp,"%f %f %f\n", t, S[n], Y[n]);}fprintf(fp,"\n");}}fclose(fp);return 0;}
グラフは gnuplot で表示できる．Excel の「ワイヤーフレーム 3D 等高線」で表示する場
合は，途中を80

--------- ここまで同じ --------FILE *fp;
fp = fopen("b-s.csv","w");for(int n=0; n<=N; n+=NN){
fprintf(fp,",%f", S[n]);}fprintf(fp,"\n");// 満期のデータを出力
fprintf(fp,"%f", T);for(int n=0; n<=N; n+=NN){
Y[n] = F(S[n]);fprintf(fp,",%f", Y[n]);}fprintf(fp,"\n");--------- 中略 --------if(m%MM==0){
fprintf(fp,"%f", t);for(int n=0; n<=N; n+=NN){
fprintf(fp,",%f", Y[n]);}fprintf(fp,"\n");}--------- これ以降同じ --------と変更すればよい．

5.12

補足

差分法は偏微分方程式を解く上で最も基本的な解法である．今回は線形の方程式のみを取り扱ったが，例えば
∂u
∂ 2u
=
+ u2
2
∂t
∂x
のような非線形方程式も差分法で解くことができる．その場合，陰解法を適用すると非線形方程式になるので，毎時間ステップごとに多次元ニュートン法を用いることになる．
問題が線形で，右辺に時間変数が入らない場合，陰解法に現れる連立一次方程式の係数行
列は時間ステップによらず毎回同じなので，そのような場合には，LU 分解などを用いる

81

ことで計算を効率化することができる．また，係数行列は一般的に疎行列となるので，その行列の特性に応じた連立一次方程式の解法を用いることで効率化を実現できる．特に今回述べたような空間 1 次元の問題の場合は，係数行列が 3 重対角行列となるので，3 重対角行列専用の方法を用いることができる．興味のある人は調べてもらいたい．
空間 2 次元以上の問題も差分法で解くことができるが，方程式を考える領域の形状が長方形でない場合には少々面倒である．領域が長方形ではない場合には，次章で紹介する有限要素法を用いることが多い．

5.13

課題

課題１ 初期条件を自由に定めて熱方程式の解を陽解法で計算し，解のグラフを表示しな
さい．
課題２ 時間方向の計算区間を [0, 1]，空間方向の計算区間を [0, 1] とし，初期条件を u0 (x) =2
sin πx，境界条件を 0 とすると，λ = 1 のときの熱方程式の解は u = e−π t sin πx となる．
この解を陰解法およびクランクニコルソン法で解いたときの最大誤差を，(∆t, ∆x) =(0.05, 0.0005), (0.1, 0.0005), (0.0005, 0.05), (0.0005, 0.1) の 4 通りについて調べ，誤差がおおよそ理論通りになっていることを確かめなさい（∆t による誤差と ∆x による誤差を完全に分離することはできないので，おおよそで確かめられればよい）．

82

6

有限要素法

本章では，偏微分方程式の解法として，差分法とならんでよく用いられる有限要素法について紹介する．有限要素法は，空間の 2 次元以上の，長方形でないような領域で方程式を考える際に威力を発揮する．有限要素法は時間発展する方程式にも適用できるが，今回は簡単のため 2 次元の定常問題を考えることにする．
有限要素法の定式化や誤差評価には，関数解析の理論が大きな役割を果たしている．この
資料では関数解析の知識は仮定していないので，適当な滑らかさがあれば，とか，関数の範囲を適当に制限すれば，等の曖昧な言い回しが多く出てくるが，これらは関数解析をベースに有限要素法を勉強すればクリアになる．興味のある人は是非勉強してもらい
たい．

6.1

準備

2 次元領域上の偏微分方程式を扱うため，記号をいくつか導入する．まず，Ω は R2 の有界領域を表すとする．具体的には平面上の多角形や円の内部をイメージするとよい．ちな
みに，数学で領域と言えば普通は開集合を指すので，Ω は境界を含まないものとする．また，Ω の境界を ∂Ω と表記する．
本資料にける関数は，特に断りがないかぎり実数の範囲で考えるものとする．また，下付
きの添え字で偏微分を表すものとする．例えば

ux =

∂u
,
∂x

uxy =

∂ 2u
∂x∂y

等である．
スカラー値関数 u に対し

 
ux
∇u =
,
uy

∆u = uxx + uyy

なる表記を用いる．ここで，∇ はナブラもしくはグラディエント，∆ はラプラシアンと読 
u1む．また，ベクトル値関数 u =に対し
u2

div u = (u1 )x + (u2 )yなる表記を用いる．div はダイバージェントと読む．
ベクトル値関数同士の積は，通常の意味でのベクトルの内積を意味するものとする．具体
 
 
u1
v1
的には，ベクトル値関数 u =
とv =
について
u2
v2

uv = u1 v1 + u2 v283と定める．また，ベクトル値関数 u に対して，|u| は通常の意味でのベクトルの長さを意 u1
味するものとする．具体的には，ベクトル値関数 u =
について
u2
q
|u| = u21 + u22
と定める．

6.2

ポアソン方程式与えられた関数 f, g に対し(−∆u = f in Ω
u = g on ∂Ωを満たす u を求める問題を考える．この方程式はポアソン方程式と呼ばれる．特に f ≡ 0のときはラプラス方程式という．この方程式は，熱伝導，薄膜の微小変形，流体の運動など，様々な現象の定式化に現れる．まずは，熱伝導および薄膜の微小変形について，どのようにポアソン方程式が現れるのかを説明する．

6.3

熱伝導

薄い板があり，板は発熱しているものとする．発熱量は単位時間，単位面積あたり f であるとする．f は位置の関数なので，発熱の仕方は場所によって異なっていても構わない．
板の形状を Ω で表し，Ω の境界では温度が g になるように冷却されているものとする．g
も場所によって異なっても構わない．状況としては，コンピュータの CPU が発熱し，周囲で冷却しているような状況をイメージするとよい．
板上の温度分布を u とする．仮定より，境界条件は u = g on ∂Ω となる．ここで，板上の微小領域 s = (x, x + h) × (y, y + h) を考える．熱の移動量は熱伝導率および温度勾配に比例するので，単位時間あたりに線分 x × (y, y + h) を横切って s から流出する熱量は
Z y+h
λux (x, t)dtyとなる．ここで，λ は熱伝導率である．同じことを s の 4 辺で考えると，単位時間あたりに s から流出する熱量はZ y+h
Z y+h
Z x+h
Z x+h
ux (x, t)dt −ux (x + h, t)dt +uy (t, y)dt −uy (t, y + h)dt
λ
y

= −λ

Z y+h yxx

Z x+h  ux (x + h, t) − ux (x, t) dt +uy (t, y + h) − uy (t, y) dt

y

x

84

ux

uy

y+h
y

x+h

x

図 3: 温度分布と熱の流入の様子となる．さて，長い時間が経過し，熱分布が定常状態に落ち着いた場合を考える．このとき，s から流出する熱量と s における発熱量が釣り合うので，
 Z y+hZ
ux (x + h, t) − ux (x, t)
1 x+h uy (t, y + h) − uy (t, y)
1
dt +
dt
−λ
h y
h
h x
h
Z x+hZ y+h1= 2
f (t1 , t2 ) dt1 dt2
h xyが成り立つ．ここで，u および f が十分に滑らかならば，テイラー展開

ux (x + h, t) − ux (x, t) = huxx (x, y) + O(h2 )uy (t, y + h) − uy (t, y) = huyy (x, y) + O(h2 )f (t1 , t2 ) = f (x, y) + O(h)を用いて h → 0 とすることで

−∆u =

f
λ

が得られる．

6.4

薄膜の微小変形

枠に固定された薄い膜（例えばゴム膜）があり，上向きに f の力がかかっているとする．ここでは，薄膜の縦方向の変形がごく微小な場合を考える．薄膜の縦方向の変位を u で表し，
枠の上では関数 g の値で固定されているものとする．よって，境界条件は u = g on ∂Ω
となる．

85

f
g
Ω
図 4: 薄膜の変形（変形の大きさは誇張して描いてある）
さて，薄膜の密度を m，厚さを d，ヤング率を Y とし，薄膜を，長さ h，バネ定数 Y d のバネが格子状に繋がり，バネの節点には重さ h2 dm の重りがついているような構造で近似する（この近似が良いかどうかは議論がある，というか本当は良くないのだが，ここでは細かいことは気にしない）．ここで，簡単のため重力は零であり，バネの自然長は零であるとする（つまり，膜は大きく引き伸ばした状態で枠に固定されているとしている）．
水平方向の変位は無視できるとすると，(x, y) の位置にある重りと (x + h, y) の位置にある重りとの間に働く張力 t は
r



t=Yd

h2 + u(x, y) − u(x + h, y)2となる．よって，この張力 t によって (x, y) の位置にある重りに下向きにかかる力を s と
すると
r

2
u(x, y) − u(x + h, y)s = Y d h2 + u(x, y) − u(x + h, y) r2h2 + u(x, y) − u(x + h, y)= Y d u(x, y) − u(x + h, y)sh
図 5: バネと重りによる近似86tとなる．同様にして，他の重りによって引っ張られる力も考え，それが上向きにかかっている力と釣り合うことから
Y d 4u(x, y) − u(x + h, y) − u(x − h, y) − u(x, y + h) − u(x, y − h) = f h2
すなわちu(x + h, y) + u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y) = −

f 2
h
Yd

が成り立つ．ここで，u が十分に滑らかならば，テイラー展開

u(x + h, y) = u(x, y) + hux (x, y) +等を用いて h → 0 とすることで

−∆u =

h2
uxx (x, y) + O(h3 )2fYdが得られる．

6.5

変分原理

物理現象において，釣り合いの状態は，何らかの物理量の極小状態と対応していることが多い．これを変分原理という（正確な定式化等は解析力学の教科書を参照すること）．例
えば，バネの一方が天井に固定され，もう一方に重りが吊るされているようなモデルを考える．
重りの重さを m，重力加速度を g ，バネ定数を k ，バネの伸びを x とすると，力の釣り合
いから
kx = mgmg

図 6: 釣り合いの位置

87

が成り立つ．また，バネの弾性エネルギーと重りの位置エネルギーの和を E とすると

1
E = kx2 − mgx2と表される．このとき，力の釣り合いからも，エネルギーの最小化からも，
mgx=kが求められる．
前節で考えた，薄膜をバネと重りの集合で近似するモデルについても全エネルギーを求めてみよう．そのため，E(x, y) を，(x, y) の位置にある重りの位置エネルギーに，その重りの右と上につながるバネの弾性エネルギーを足したものとする．この E(x, y) を足し合わせれば全エネルギーが求められる．ニュートンの運動法則から，重りに単位面積あたり fの力がかかっているということは，f /(dm) の加速度がかかっていることと同じなので，


2 
1
2
E(x, y) = Y d h + u(x, y) − u(x + h, y)
2


2 
1
f (x, y)2
+ Y d h + u(x, y) − u(x, y + h)u(x, y)− h2 dm ·
2
dm
2

1
u(x + h, y) − u(x, y)22= Y dh + Y dh
2
h

2
1
u(x, y + h) − u(x, y)2+ Y dh− h2 f (x, y)u(x, y)2hとなり，これまでと同様，テイラー展開を用いると

1
E(x, y) = Y dh2 + Y dh2 |∇u(x, y)|2 − h2 f (x, y)u(x, y) + O(h3 )2が成り立つ．よって，これを足し合わせて h → 0 として積分の定義を用いると，全エネ
ルギーは
Z
Z
1
2
Y d|Ω| + Y d |∇u(x, y)| dxdy − f (x, y)u(x, y) dxdy2ΩΩとなることがわかる．よって，変分原理によれば，境界上で g と一致する関数の中で

−∆u =

f
Yd

を満たす u を求める問題と，汎関数
Z
Z
1
fu
2
|∇u| dxdy −
dxdy
2 Ω
Ω Yd
の最小値を与える u を求める問題は同値であると期待できる．ここで，汎関数とは，関数に対してスカラーを返すような写像を意味する（普通の関数は，スカラーもしくはベクト
ルに対してスカラーを返す）．

88

さて，

f
を改めて f と置くと，
Yd−∆u = fを満たす u を求める問題と

1
E(u) =
2

Z

Z|∇u| dxdy −2Ω

f u dxdyΩの最小値を与える u を求める問題が同値であると期待できる．
ポアソン方程式の解が汎関数 E(u) の最小値を与える関数に一致するということは，19 世
紀前半頃には知られていたようだが，厳密な証明が得られたのは 19 世紀末になってからである．現代では，関数解析の理論を用いて示すことができる．そもそも，汎関数においては，最小値の存在すら自明ではない．汎関数が下に有界なら最小値があるのは当たり前だと思うかもしれないが，これはそれほど明らかな話ではないのである．
例えば C を，u(0) = u(1) = 0 を満たすような一変数連続関数の集合とすると，汎関数

Z 1
G(u) =

2
u(x) − 1 dx

0

は C の中では最小値を取らない（証明は考えてみよう）．
詳しい説明は省くが，領域 Ω や f や g にある程度の滑らかさを仮定し，関数 u の範囲を適当に定めれば，ポアソン方程式の解は汎関数 E(u) の最小値を与える関数に一致することが示される．
さてここで，汎関数 E(u) を最小にする u がどのような性質を持つかを考えてみよう．そこで，適当な滑らかさを持ち，∂Ω で零となる関数の集合を H0 とし，φ ∈ H0 についてE(u + λφ) を考えてみると，
ZZ
12|∇(u + λφ)| dxdy − f (u + λφ)dxdy
E(u + λφ) =
2 Ω
Ω
Z
Z
Z
1
λ2
2
=
|∇u| + λ ∇u∇φ dxdy +|∇φ|2 dxdy
2 Ω
2 Ω
Ω
Z
Z
− f u dxdy − λ f φ dxdyΩΩとなる．ここで上式を λ の二次式と見ると，E(u + λφ) は λ = 0 で最小にならねばならないので，λ の一次の項は零でなければならない．すなわちZZ∇u∇φ dxdy =f φ dxdyfor ∀φ ∈ H0
Ω

Ω

89

が成り立つ．このような条件を満たす u をポアソン方程式の弱解という．また，弱解に出てくる φ ∈ H0 を試験関数という．
もともとポアソン方程式は−∆u = fという形で 2 階微分が入っている．それに対し，弱解は 1 階微分が存在すればいいので，弱解の方が条件が緩い．そういう意味で弱解という名称が使われている．
方程式に現れる微分の階数だけ連続微分可能で，かつ境界までこめて連続な解を古典解という．ポアソン方程式の場合は C 2 級かつ境界までこめて連続な解が古典解となる．一般
的に，方程式の古典解は弱解になることが多いが，逆に弱解が古典解になるかどうかは状況による．方程式によっては，弱解しか存在しないこともある．しかし，古典的な意味では元の方程式を満たさないとしても，弱解には物理的に重要な意味があることが多いし，古典解を考えるよりも弱解を考える方が自然な場合も多い．
例えば，Ω = (0, 1)2 において波動方程式(
uxx − uyy = 0 in Ω

u = 0 on ∂Ω
を考える．このとき，h(t) を周期 1 の周期関数とすると

u(x, y) = h(x + y) − h(x − y) − h(−x + y) + h(−x − y)が波動方程式を満たすことは容易に確かめられる．このとき，h が滑らかならば u は古典解になるが，必ずしも h が滑らかでない場合にも u を解として扱いたいというのは自然な要求である．そのような場合には，弱解を考えれば，必ずしも 2 回微分できないような関数も解として考えることができる（h が不連続点を持つような関数の場合は，さらに解の範囲を広げなければならないが）．

6.6

ガウスの発散定理とグリーンの公式

今後，多用するので，ここでガウスの発散定理とグリーンの公式について説明しておく．
関数は十分に滑らかであるとする．
ガウスの発散定理

ベクトル値関数 u についてZZ
div u dxdy =Ω

u · n ds∂Ωが成り立つ．ここで n は外向き単位法線ベクトル，ds は ∂Ω 上の線要素である．つまり，右辺は境界上における u と n の内積の積分を表す．

90

n= cosθ
sinθ

(

Ω

)

θ
n= 0
-1

( )

図 7: 外向き単位法線方向ベクトル
例えば，(x, y) = (r cos θ, r sin θ) とし，onΩ = (x, y) r < 2, y > 0 ,
とすると
Z

Z
div u dxdy =

Ω

 2
x
u=
y2

Z πZ 2(2x + 2y) dxdy = 2Ω(r cos θ + r sin θ)r drdθ =
0

0

32
3

となり，一方で，∂Ω のうち円弧部分を γ1 ，直線部分を γ2 とすると
ZZ  2 Z  2  
x
cos θ
x
0
u · n ds =
·
ds +
·
ds
2
2
y
sin θ
y
−1
∂Ω
γ1
γ2
 

Z π 2Z 2  2  
2 cos2 θ
cos θ
x
0
=
·
2 dθ +·dx
2
2
2 sin θ
sin θ
0
−1
0
−2
Z π
32
=8
(cos3 θ + sin3 θ) dθ + 0 =30となり，確かに成り立っていることがわかる．ただし，線要素 ds を通常の積分変数で置き換える際には注意を要する．s を境界上の移動距離に比例して増加する変数とすると，γ1 上では s = 2θ となるので，ds = 2dθ となる．
ガウスの発散定理の証明についてはベクトル解析の教科書で勉強してもらいたい．ここでは，ガウスの発散定理の直感的な意味について説明する．
 
u1
平面上を水が流れているものとし，(x, y) における流速を u =で表す．また，平面
u2
からは水が湧き出しているものとする（底から地下水が滲み出してくるイメージ）．さて，

91

u•n

u
θ

n
γ

図 8: 境界を横切る流量
境界 γ を横切って流出する流量は，外向き単位法線ベクトル n を用いて，単位長さ単位時
間あたり|u| cos θ = |u||n| cos θ = u · nと表される．よって，Ω から流出する流量の合計は単位時間あたり
Z
u · n ds∂Ωとなり，これは Ω 全体における単位時間あたりの湧き出し量と等しい．一方，Ω 内の微小領域 (x, x + h) × (y, y + h) から流出する流量は

Z y+h

Z y+h
u1 (x + h, t)dt −y

Z y+h

=h
y
2

Z x+h
u2 (t, y + h)dt −u1 (x, t)dt +yu1 (x + h, t) − u1 (x, t)
dt +
h

Z x+h

x
Z x+h
x

u2 (t, y)dtxu2 (t, y + h) − u2 (t, y)dth3= h div u + O(h )となるので，これはこの微小領域における湧き出し量と等しい．これを Ω 全体で足し合わせると全体の湧き出し量になるので，h → 0 の極限を考えると
ZZ
div u dxdy =u · n dsΩ∂Ωが成り立つ．
ガウスの発散定理から，次のグリーンの公式が成り立つ．
グリーンの公式

スカラー値関数 u およびベクトル値関数 v について
ZZ
Z∇u · v dxdy =uv · n ds − u div v dxdy
Ω

∂Ω

Ω

が成り立つ．

92 v1
グリーンの公式は，v =と置くと
v2
div(uv) = (uv1 )x + (uv2 )y = ux v1 + u (v1 )x + uy v2 + u (v2 )y = ∇u · v + u div vとなるので，左辺にガウスの発散定理を適用することにより示すことができる．
グリーンの公式は一変数関数の部分積分の拡張になっている．
実際，Ω を長方形 (a, b)×(0, 1) wとし，u と v =が x だけの関数とすると，
0Z
Z 1Z bZ b∇u · v dxdy =ux w dxdy =ux w dx
Ω
0
a
a
Z
Z 1
Z 1
uv · n ds =u(b)v1 (b) dy −u(a)v1 (a) dy = [uw]ba∂Ω00Z
Z 1Z bZ b
u div v dxdy =uwx dxdy =uwx dxΩ0aaより，グリーンの公式は一変数関数の部分積分の公式と一致する．また，グリーンの公式自体も部分積分と呼ぶことがある．

6.7

弱定式化

グリーンの公式の応用として，ポアソン方程式の古典解が弱解になることを示そう．u を
ポアソン方程式(−∆u = f in Ω

u = g on ∂Ωの古典解とする．ここで，第一式に任意の φ ∈ H0 をかけて Ω で積分すると
ZZ− ∆uφ dxdy =f φ dxdyΩΩとなる．ここで，部分積分（グリーンの公式）を適用すると
ZZZ− ∆uφ dxdy = −φ∇u · n ds + ∇u∇φ dxdy
Ω

∂Ω

Ω

が成り立つが，φ ∈ H0 より φ が ∂Ω で零となることから
ZZ∇u∇φ dxdy =f φ dxdyforΩ∀φ ∈ H0Ωが得られる．よって u は弱解になる．
このように，部分積分（グリーンの公式）を利用して弱解が満たす条件式を求めることを，弱定式化という．

93

6.8

リッツ法とガレルキン法

本節以降では，簡単のため境界条件は零とする．すなわち，有界領域 Ω と関数 f について，次のポアソン方程式を考える．
(−∆u = f in Ω

u = 0 on ∂Ω
特殊な場合を除くと，この方程式は数学的には解けないので，解の様子を知りたい場合には数値計算が用いられることが多い．
さて，この方程式の近似解を求めたい場合に，近似解を有限個の関数 φk ∈ H0 , k =1, 2, · · · , n の一次結合
n
X
uh =
ck φk
k=1

で近似するという発想は自然である．このような方法を有限要素法といい，得られた近似解を有限要素解という．また，関数 φk (k = 1, 2, · · · n) を有限要素基底という．有限要素法においては，係数 {ck } をどのように決めるべきかが問題となる．係数を決める代表的な方法としては，リッツ法とガレルキン法がある．これらについて順に説明する．
まず，ポアソン方程式の解 u は汎関数
ZZ
1
2
E(u) =|∇u| dxdy − f u dxdy
2 ΩΩの最小値を与える関数として特徴づけられたことを思い出そう．よって，{φk } の一次結合全体で表される関数の集合を Sh とし，Sh の中で E(u) を最小にするものを近似解とす
ることにする．そこで，任意の φ ∈ Sh を取り，以前と同様に E(uh + λφ) を考えてみると，λ = 0 で最小値を取ることからZZ∇uh ∇φ dxdy =f φ dxdyfor∀φ ∈ ShΩΩが得られる．上式は，φ = φj , (j = 1, 2, · · · n) で成り立てば任意の φ ∈ Sh で成り立つので，近似解を求めるための条件はZZ∇uh ∇φj dxdy =f φj dxdy
j = 1, 2, · · · , nΩΩとなる．さらに

uh =

n
X
k=1

94

ck φ kを代入すると，条件は
Z
Z
n
X
ck
∇φk ∇φj dxdy =f φj dxdy
k=1Ω

j = 1, 2, · · · , nΩとなる．これは n 次元連立一次方程式になり，
a11 a12 · · · a1nZ a21 a22 · · · a2n 


A= .
,
ajk =
∇φk ∇φj dxdy,
..
.. 
...
 ..
.
. 
Ω

an1 an2 · · · ann
 
 
c1
b1
 c2 
 b2   c =  . ,b =  . ,.
.
 .. 
cn
bn

Z
bj =

f φj dxdyΩと置くとAc = bと書ける．このとき，行列 A を剛性行列，ベクトル b を外力ベクトルという．このように，方程式の解がある汎関数を最小化するときに，有限次元空間の中からその汎関数を最小化するものを探して近似解とする方法をリッツ法という．
さて，ここで，近似解を求めるもう一つのアプローチについて説明する．近似解を有限個の関数 φj ∈ H0 , j = 1, 2, · · · , n の一次結合

uh =

n
X

ck φk

k=1

で表すのはリッツ法と同様だが，今度は，近似解 uh を，厳密解との誤差を最小にするものとして求める．誤差の測り方には色々あるが，ここではZ|∇(u − uh )|2 dxdyΩ

を最小にするように uh を定める．このとき，任意の φ ∈ Sh を取ると
Z
Z
Z
2
2
|∇(u − uh + λφ)| dxdy =|∇(u − uh )| dxdy + 2λ ∇(u − uh )∇φ dxdy
Ω
Ω
Ω
Z
2
|∇φ|2 dxdy+λΩ

が λ = 0 で最小になることからZ∇(u − uh )∇φ dxdy = 0
Ω

95

for

∀φ ∈ Sh

が成り立つ．ここで，弱解が満たす式を用いると，φ ∈ Sh のとき φ ∈ H0 よりZ
ZZ∇uh ∇φ dxdy =∇u∇φ dxdy =f φ dxdyfor∀φ ∈ ShΩΩΩとなり，リッツ法の手順と一致する．このように，誤差を最小化するものとして近似解を求める方法をガレルキン法という．
今回の問題については，リッツ法とガレルキン法は一致するので，合わせてリッツ・ガレ
ルキン法と言う．しかし一般的には，解が何らかの汎関数を最小にするとは限らないので，リッツ法は常に適用可能ではないし，ガレルキン法においても，どのように誤差を測るかによって求まる近似解は異なるので，一般的にはリッツ法とガレルキン法が一致するとは限らない．

6.9

有限要素基底前節において，有限要素基底 φj (j = 1, 2, · · · n) については，境界で零になるという以外に特に制限は無かった．そういう意味では何でもよいのだが，実際には計算を簡便にするため，領域 Ω を多数の三角形に分割し，その上で連続な区分一次関数を考えることが多い．
具体的には，まず領域 Ω を多数の三角形に分割し，Ω の内部にある節点を p1 , p2 , · · · , pn とする．また，境界上の節点を pn+1 , pn+2 , · · · とする．そして，pj で 1 になり，他の節点では零になり，各々の三角形上では一次関数となるような連続関数を φj とする．ひとまずは有限要素基底には境界上で零になるという制限があるので，φj (j = 1, 2, · · · n) を用いるが，境界条件によっては φn+1 以降も基底として用いることがある．また，φj (j = 1, 2, · · · n)の一次結合全体で表される関数の集合を Sh とする．
分割されたそれぞれの三角形を有限要素，もしくは単に要素という．有限要素分割には三角形分割の他にも長方形分割などもあるが，三角形による分割であることを強調したい場
合には三角形要素などともいう．また，境界がカーブしている場合には，まず領域を多角形で近似してから三角形分割を行うものとする．
有限要素基底を用いて有限要素解

uh =

n
X

ck φ kk=1が得られたとすると，uh (pk ) = ck となるので，係数 ck そのものが u(pk ) の近似になっている．

96

以下は三角形による領域分割と有限要素基底の例である．

p3
p4
p2

Ω

p1

p5

図 9: 三角形による分割
1

1

ϕ1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

1
0.8

ϕ3

1
0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

1
0.8

ϕ5

0.6
0.4
0.2
0

図 10: 有限要素基底

97

ϕ2

ϕ4

6.10

節点一覧表および要素節点対応表

有限要素法では，三角形分割のデータを扱うために，節点一覧表と要素節点対応表を用いる．例えば，Ω = (0, 1) × (0, 1/2) とし，Ω を 16 個の三角形に分割し，図 11，図 12 のよ
うに三角形分割して節点および要素に番号を振ったとする．
節点番号は，内点から番号を振り，その後に境界上の節点に番号を振るとプログラムが簡単になる．ただし実際は，使い易さの点から，内点も境界上の節点も区別せずに番号を振っていくことが多い．その場合のアルゴリズムについては後で述べる．
節点ごとに，その座標と節点の属性（0:内点，1:境界上の点）を記したものが節点一覧表
である．また，三角形要素ごとに，頂点の節点番号を記したものが要素節点対応表である．図 11，図 12 で例に挙げた三角形分割の場合，節点一覧表と要素節点対応表はそれぞれ，表 1，表 2 のようになる．要素節点対応表に記す頂点は，反時計回りに記しておくも
のとする．

11

12

13

14

15

9

1

2

3

10

4

5

6

7

8

図 11: 節点番号

13

14
9

5

15
10

6
1

16
11

7
2

8
3

図 12: 要素番号

98

12

4

表 1: 節点一覧表
節点番号 x 座標 y 座標
1
0.25
0.25
2
0.50
0.25
3
0.75
0.25
4
0.00
0.00
5
0.25
0.00
6
0.50
0.00
7
0.75
0.00
8
1.00
0.00
9
0.00
0.25
10
1.00
0.25
11
0.00
0.50
12
0.25
0.50
13
0.50
0.50
14
0.75
0.50
15
1.00
0.50

属性
0
0
0
1
1
1
1
1
1
1
1
1
1
1
1

表 2: 要素節点対応表
要素番号 頂点 1 頂点 2 頂点 3
1
4
5
1
2
5
6
2
3
6
7
3
4
7
8
10
5
1
9
4
6
2
1
5
7
3
2
6
8
10
3
7
9
9
1
12
10
1
2
13
11
2
3
14
12
3
10
15
13
12
11
9
14
13
12
1
15
14
13
2
16
15
14
3

99

6.11

剛性行列と外力ベクトルの計算

ここでは，剛性行列
a11 a12 · · · a1n
 a21 a22 · · · a2n 


A= .
,
..
.. 
...
 ..
.
. 
an1 an2 · · · ann
および外力ベクトル b1
 b2  
b =  . , .. 

Z∇φk ∇φj dxdy

ajk =
Ω

Z
bj =

f φj dxdyΩ

bnを節点一覧表と要素節点対応表から計算する方法について述べる．本節以降，分割された三角形要素を τl (l = 1, 2, · · · , m) とする．
剛性行列 A や外力ベクトル b を求める際に，ajk や bj を順に求めるというのは実はあまり効率的ではない．例えば，ajk を求めるためには pj と pk を両方含む三角形要素を検索しなければならず，手間がかかる．それよりも，各要素ごとの足し合わせで A や b を求める方法がよく用いられている．具体的には，
 (l) (l)(l)
a11 a12 · · · a1n
 (l) (l)
Z
(l) 

a
a
·
·
·
a
21
22
2n 
(l)
(l)


ajk =
∇φk ∇φj dxdy
A = .
..
..  ,
..
.
.
. 
τl
 ..
(l)
(l)
(l)
an1 an2 · · · ann (l) b1
 (l) 
Z
b2 
(l)
(l)


bj =
f φj dxdyb =  . ,τl
 .. 
(l)
bn
として

A=

m
X

(l)

A ,

l=1

b=

m
X

b(l)

l=1

により A と b を求める．かえって複雑になったように思うかもしれないが，実はこうすることで効率的に計算を行うことができる．実際の手順としては次のようになる．

1. A および b を零行列および零ベクトルで初期化する．
2. l を順に 1 から m まで動かす．

100

3. 要素節点対応表により τl の 3 つの頂点の節点番号を求め，それを s1 , s2 , s3 とする．
(l)
(l)
(l)
(l)
(l)
(l)
(l)
(l)
(l)
4. as1 s1 , as1 s2 , as1 s3 , as2 s1 , as2 s2 , as2 s3 , as3 s1 , as3 s2 , as3 s3 , の 9 つを計算し，A に足し
(l)
(l)
(l)
(l)
(l)
(l)
込む（as1 s2 = as2 s1 , as2 s3 = as3 s2 , as3 s1 = as1 s3 なので，実際に計算するのは 6 つ）．
(l)
(l)
(l)
5. bs1 , bs2 , bs3 の 3 つを計算し，b に足し込む．
ただし，4. と 5. については s1 , s2 , s3 のうち境界上の節点については適用しない．
(l)

(l)

それでは具体的に as1 s1 , as1 s2 , · · · を求めてみよう．節点一覧表より座標を求め，ps1 =(x1 , y1 ), ps2 = (x2 , y2 ), ps3 = (x3 , y3 ) とする．φs1 は τl 上では一次関数なので

φ s 1 = a1 x + b 1 y + c 1と書ける．
ここで，φs1 は ps1 で 1，ps2 と ps3 では 0 となるので

a1 x 1 + b 1 y 1 + c 1 = 1
a1 x 2 + b 1 y 2 + c 1 = 0
a1 x 3 + b 1 y 3 + c 1 = 0となる．ここで
x1 y1 1P = x2 y2 1
x3 y3 1と置くと，クラメールの公式より

y2 − y3y3 − y1
y1 − y21 P −1 =
x1 − x3
x2 − x1  x3 − x2|P |x2 y3 − x3 y2 x3 y1 − x1 y3 x1 y2 − x2 y1
1

ϕs

0.8

1

0.6
0.4
0.2
0

ps

τl

2

ps

3

図 13: φs1 の様子

101

ps

1

となるので

 
 
a1
1
 
−1  
 b1  = P 0
c1
0

より a1 と b1 が求まる．具体的には

a1
b1

!
=

y2 − y3x3 − x2

1|P |!

となる．ここで，τl の面積を |τl | で表すと，ps1 , ps2 , ps3 が反時計回りに並んでいれば|P | = 2|τl |が成り立つので，結局，

a1 =y2 − y3,2|τl |b1 =x3 − x22|τl |

が成り立つ．同様に

y3 − y1,2|τl |
y1 − y2,a3 =2|τl |

x1 − x3
2|τl |
x2 − x1
b3 =2|τl |a2 =b2 =が成り立つ．これらを用いることにより，
Z(l)
as 1 s 1 =∇φs1 ∇φs1 dxdy = |τl |(a21 + b21 )Zτl∇φs1 ∇φs2 dxdy = |τl |(a1 a2 + b1 b2 )a(l)s1 s2 =
τl·········
等と計算することができる．
図 11，図 12，表 1，表 2 で定められる三角形分割について，剛性行列を求めると，以下のようになる．

4 −1 0
A = −1 4 −1
0 −1 4
プログラムを組む前に，まずは手計算で，これと同じ結果になるか確かめておくとよい．
Z(l)
次に，外力ベクトルを求めるためには bj =f φj dxdy を計算する必要がある．しかし，
τlf が単純な関数である場合を除くと正確に値を求めることは難しい．数値積分を用いるの102

も一つの方法であるが，それよりも f を区分一次関数で近似して外力ベクトルを求める方法がよく用いられる．具体的には，f の pk での値を fk とし，f を φk （境界上の基底も含む）の一次結合でXf≈fk φ k
k

と近似して

Z
(l)
bj =

f φj dxdy ≈

X

τl

k

Z
fk

φk φj dxdyτlと計算する．このためには，τl の 3 つの頂点の節点番号を s1 , s2 , s3 としたとき，
ZZ
φs1 φs1 dxdy,
φs1 φs2 dxdy, · · ·τlτl等がわかればよい．これは，以下を用いて計算できる（ちゃんとそうなることは確かめて
みよう）．


|τ |


(j = k)6
φsj φsk dxdy =τ |τ | (j =
̸ k)
12

Z

6.12

境界も含めて計算する方法

前節までに述べた方法は，領域の内点のみについて計算する方法であったが，この場合
は，新たに境界を設定したり，領域の形状を変更したりしたときに，内点だけに番号を振り直すという操作が必要になり，面倒である．そこで，内点も境界点も同時に計算する方法について説明する．
まず最初に，内点も境界点も区別せずに，既に説明した方法を用いて剛性行列と外力ベクトルを計算する．行列およびベクトルの次数は境界点も含めた分点の数となる．その上で，境界点に関しては値が零となるように剛性行列と外力ベクトルを変形する．具体的には，pk を境界上の点としたとき，以下のように書き換える．



第k列 ↓







A = 0 · · ·






0
..
.
0
0 1 0 ···
0
..
.








0 ,






0

103

 
 .. 
.
 
 
 
 
b = 0 ← 第 k 行
 
 
 
 .. 
.

これにより，境界上の点を有限要素基底から切り離し，零に固定することができる．

6.13

三角形分割の方法

本節では，与えられた領域を三角形分割する方法について概略を述べる．三角形分割することを，メッシュを切る，という．領域が直角だけで構成されるような単純な形状なら，図 14 のように直角三角形でメッシュを切るのが楽である．境界が曲がっている場合
でも，工夫すれば図 15 のようにメッシュを切ることができる．しかしながら，メッシュの良し悪しについては，三角形の大きさが出来るだけ均一で，しかも三角形の歪が小さい方が良いとされてるので，図 16 のような規則性のないメッシュも用いられる．このようなメッシュは，境界から中心に向けて三角形で埋めていったり，互いに反発する多くの点を領域内にばらまき，定常状態になったところで線で結んでメッシュを作る（ドロネー分割という手法を用いることが多い），等の方法によって作られているが，一長一短がある．FreeFem++ というフリーソフトを用いれば，境界のデータを与えるだけで手軽にメッシュを切ることができる．

6.14

実際の数値計算

有限要素解を求める問題は，最終的に連立一次方程式に帰着された．よって，実際に数値計算する上では，連立一次方程式の効率的な解き方についても考えなければならない．一番簡単なのはガウスの消去法を用いる方法である．しかし，基底関数の数が数万以上になってくると，SOR 法や CG 法などの反復法の方がガウスの消去法よりも効率的になってくる．また，有限要素法の性質から剛性行列は疎行列（非零成分が少ない行列）になる
ので，メモリを出来るだけ有効に使うテクニックも色々と考えられる．ちなみに，現実の図 14: 直角三角形によるメッシュ

104

図 15: 規則性のあるメッシュ

図 16: FreeFem++を用いて作成したメッシュ
問題の解析においては，基底関数の数が数百万個から数億個程度の計算は普通に行われて
いる．

6.15

有限要素法のプログラム

以下に有限要素法のプログラムを示す．汎用性を考えると，メッシュを切って節点一覧表と要素節点対応表を作るプログラムと，節点一覧表と要素節点対応表のデータを読み込
んで有限要素解を求める操作は分けた方が良いのだが，今回は簡単のため同時に行っている．領域は正方形 [0, 1]2 とし，直角二等辺三角形でメッシュで切っている．

105

#include <stdio.h>#include "Linear.h"// 右辺の関数
double f(double x, double y){return 1;}
int main(){int N = 20; // 分割数
int np = (N+1)*(N+1); // 分点数
int ne = 2*N*N;// 要素数// 節点一覧表
Vector px(np);
Vector py(np);
IntVector ps(np);// 要素接点対応表
IntVector ep1(ne);IntVector ep2(ne);IntVector ep3(ne);Vector ff(np); // 右辺の関数値を格納// 表の作成
int p = 1;
int e = 1;
for(int j=0; j<=N; j++){for(int i=0; i<=N; i++){double x = i/(double)N;
double y = j/(double)N;
次ページへ続く…

106ff[p] = f(x,y);px[p] = x;py[p] = y;if(i==0 || i==N || j==0 || j==N){
ps[p] = 1; // 境界点} else {ps[p] = 0; // 内点}
if(i<N && j<N){ep1[e] = p;//／|ep2[e] = p+1;//
／ |ep3[e] = p+N+2; // (x,y) ---e++;
ep1[e] = p;//---ep2[e] = p+N+2; //| ／ep3[e] = p+N+1; // (x,y) |／
e++;
}
p++;
}
}
Matrix A(np, np); // 剛性行列
Vector b(np); // 外力ベクトル
for(int i=1; i<=np; i++){ // 剛性表列と外力ベクトルを初期化
b[i] = 0;
for(int j=1; j<=np; j++){A[i][j] = 0;}}
次ページへ続く…

107

for(int e=1; e<=ne; e++){int s1 = ep1[e], s2 = ep2[e], s3 = ep3[e];double x1 = px[s1], x2 = px[s2], x3 = px[s3];double y1 = py[s1], y2 = py[s2], y3 = py[s3];double tau = ((x2-x1)*(y3-y1) - (x3-x1)*(y2-y1))/2; // 面積
double a1 = (y2-y3)/2/tau, b1 = (x3-x2)/2/tau;
double a2 = (y3-y1)/2/tau, b2 = (x1-x3)/2/tau;
double a3 = (y1-y2)/2/tau, b3 = (x2-x1)/2/tau;
A[s1][s1] += tau*(a1*a1 + b1*b1); // 剛性行列に足し込む
A[s1][s2] += tau*(a1*a2 + b1*b2);
A[s1][s3] += tau*(a1*a3 + b1*b3);
A[s2][s1] += tau*(a2*a1 + b2*b1);
A[s2][s2] += tau*(a2*a2 + b2*b2);
A[s2][s3] += tau*(a2*a3 + b2*b3);
A[s3][s1] += tau*(a3*a1 + b3*b1);
A[s3][s2] += tau*(a3*a2 + b3*b2);
A[s3][s3] += tau*(a3*a3 + b3*b3);// 外力ベクトルの計算
b[s1] += tau*(ff[s1]/6 + ff[s2]/12 + ff[s3]/12);b[s2] += tau*(ff[s1]/12 + ff[s2]/6 + ff[s3]/12);b[s3] += tau*(ff[s1]/12 + ff[s2]/12 + ff[s3]/6 );}// 境界値の分離と固定
for(int p=1; p<=np; p++){if(ps[p]>0){
for(int i=1; i<=np; i++){A[p][i] = 0;A[i][p] = 0;}A[p][p] = 1;b[p] = 0;}}
次ページへ続く…

108

Vector u = A.Gauss(b);// ガウスの消去法

FILE *fp;
fp = fopen("fem.dat", "w");for(int e=1; e<=ne; e++){int s1 = ep1[e], s2 = ep2[e], s3 = ep3[e];fprintf(fp, "%f %f %f\n", px[s1], py[s1], u[s1]);
fprintf(fp, "%f %f %f\n", px[s2], py[s2], u[s2]);
fprintf(fp, "%f %f %f\n", px[s3], py[s3], u[s3]);
fprintf(fp, "%f %f %f\n", px[s1], py[s1], u[s1]);
fprintf(fp, "\n\n");}fclose(fp);return 0;}
データは，空間内の三角形を順に出力するようになっているので，gnuplot の splot で表示すると，以下のように解が表示される．
メッシュ分割さえできれば，正方形領域に限らず，様々な形状の領域においても数値解を計算
することができる．例えば図 18 は，偏心円環領域 (x, y) x2 + y 2 < 1 (x − 1/4)2 + y 2 > 1/4
において f = 1 のときのポアソン方程式の解を有限要素法で求めた結果である．

0.1
0.08
0.06
0.04
0.02
0
1
0.8
0.6
0.4
0.2
0

0

0.2

0.4

0.6

0.8

1

図 17: 正方領域におけるポアソン方程式の数値解

109

0.1
0.08
0.06
0.04
0.02
0
-1
-0.5
0
0.5
-1

1

-0.5

0

0.5

1

図 18: 円環領域におけるポアソン方程式の数値解

6.16

補足

本資料では簡単のため，境界条件を零としたが，もちろん，境界条件が零以外の問題も解くことができる．興味のある方は調べてみるとよいだろう．
有限要素法の誤差解析については概略のみ触れることにする．有限要素法のガレルキン法的な構成において，有限要素解がZ|∇(u − uh )|2 dxdyΩ

を最小にする uh として求められたことを思い出そう．そこで，全ての分点で u と値が一致し，各三角形要素上では一次関数となるような関数 Πu を考える．Πu を，u の区分線
形補間，区分一次補間などと呼ぶ．Πu は有限要素基底の一次結合で書けるので，ZZ
2|∇(u − uh )| dxdy ≤|∇(u − Πu)|2 dxdyΩΩが成り立つ．これはつまり，真の解と有限要素解の間の誤差が，真の解と Πu の間の誤差
で押さえられることを示している．このように，有限要素法の誤差解析は補間誤差を介して解析することが多い．
有限要素法は関数解析をベースに構築されているので，厳密な構成と誤差解析には関数解析が必要不可欠である．本章では敢えて関数解析を用いていないので，曖昧でスッキリしない部分も多かった．ちゃんと理解したいという人は，頑張って勉強してもらいたい．

110