「応用数学」第 7 章『発展系の数値解析』
桂田 祐史

1995 年〜, 2024 年 8 月 17 日

目次
この文書は、放送大学の講義科目「応用数学」 (藤田 [?]) のテキスト第 7 章 (ここだけ桂田が執筆した) に少し手を入れたものである。一部第 6 章の内容に言及している。
桂田のゼミ生で 6 章の内容を見たい人はメールで依頼すること (他の資料で代用できるので、7 章を理解するために、必ずしも 6 章を見る必要はないが、言及していない部分に、有益なものが多くある)。

コンピューター・シミュレーションの意義

1

さまざまな現象が、数学的には微分方程式の形で表現されることは、すでにこの科目でも、何度も見て来た通りである。ところが、微分方程式の解を具体的に求めることは、常に可能であるとは限らない。たとえば第 4 章1 「惑星の運動の数理」において、太陽の周りを万有引力の法則に従って運動する一つの惑星の軌道を求めることはできたが、二つ以上の惑星を同時に考慮した場合には、同様の手段で解を求められないことが証明されている (三体問題の不可能性)。この「解を求めることができない」というのは、解が存在しないということでも、まして現在の人類の持っている数学の知識が不十分であるために解けない、という意味でもなく、良く知られている関数を使って解を表現する公式がない、ということである。
この冷厳な事実に気付いた後に数学者が行ったことは、一つには微分方程式の解として、逆に新しい関数を定義することであった。この努力の結果、一般に特殊関数とよばれる多くの関数が導入されることとなり、応用数学の中で一定の位置を占めている。
もう一つの方向は、解が存在するかどうか、解が存在したとしてそれが一意 (ただ一つに限る) かどうか、さらには解がどのような定性的な性質を持つか、を論ずることである。この努力は、現在でも盛んに続けられていて豊富な理論が築きあげられている。
しかしながら、微分方程式の解の定性的な性質がわかっただけで、定量的な情報が一切得られないのでは、満足できないことが多いのも、また確かである。さらに微分方程式を見ただけ
で、解の性質を知るのもそう容易なことではない。そのため、様々な近似解法が工夫されることになった。その中には数値解法もあったが、人力での数値計算でできることにはおのずと限界があった。
コンピューターの出現により、大規模な数値計算が可能になったことで、このあたりの事情は大きく変わることになった。コンピューター・シミュレーションによって、実際の応用に大1

藤田宏、「応用数学」、放送大学出版会の第 4 章のこと。

1

いに役立つ、定量的な情報を入手できるようになった。それだけにとどまらず、豊富なコン
ピューター・シミュレーションによって、それなしではわかり難かった、問題の定性的な性質が明らかになることも多い。現在の応用数学にとって、コンピューター・シミュレーションは必要不可欠なものとなっている。
この章では、微分方程式、特に偏微分方程式の代表的な数値解法である差分法について学ぶことにする。差分法とは

1. 問題を考えている領域を格子で区切って、各格子点上での値を求めることを目標にする。
2. 微分方程式に現れる微分係数を差分商で置き換えて (微分方程式を近似する) 差分方程式を作る。
という方針で離散的な問題 (差分方程式と呼ばれる — 後述) を導き、それを数値計算で解いたもの (差分解と呼ばれる) を、元の問題の近似解として採用する、というのが基本的な考え方であり、極めて汎用性が高い。差分法が数値解法の基礎として、応用数学に役立つものであることは当然であるが、差分法自身にひそむ数理の探求も、重要な応用数学の主題である。特に、微分方程式の解の存在、一意性、解の定性的性質の研究を背景にしている部分が少なくな
いことも注意しておこう。

差分近似

2

ここでは、後の議論に必要になる、いくつかの差分近似式を導いておこう。
関数 f の点 x における微分係数 f ′ (x) の定義

f ′ (x) = lim

(1)

h→0

f (x + h) − f (x)hを見れば、|h| が十分小さいとき、f ′ (x) の近似値として前進差分商

f (x + h) − f (x)h(2)

を採用することを真っ先に思いつく。これは f ′ (x) に対する前進差分近似と呼ばれる。ここでは f の Taylor 展開を用いて、さらに色々な近似式を導いてみよう。
f が点 x の近傍で十分滑らかとすると、(3)f (x + h) = f (x) + f ′ (x)h +f ′′ (x) 2 f ′′′ (x) 3 f (4) (x) 4h +h +h + ···23!
4!

が成り立つ。h を −h で置き換えると、(4)f (x − h) = f (x) − f ′ (x)h +f ′′ (x) 2 f ′′′ (x) 3 f (4) (x) 4h −h +h − ··· .
23!
4!

(3) から(5)f ′ (x) =f (x + h) − f (x)+ O(h)
h2(h → 0)が導かれる。すなわち前進差分近似の誤差は O(h) であることがわかる。ここで現れた O(h)はいわゆる Landau の記号である (補足 (1) (p. 15) で説明してある)。一方 (4) からは

f ′ (x) =(6)f (x) − f (x − h)+ O(h)h(h → 0)が得られる。この右辺の第 1 項を後退差分商と呼ぶ。f ′ (x) に対する後退差分近似の誤差はやはり O(h) である。
次に (3) から (4) を引いて整理すると

f ′ (x) =(7)f (x + h) − f (x − h)+ O(h2 )2h(h → 0)が得られる。この右辺第 1 項を 1 階中心差分商と呼ぶ。f ′ (x) に対する 1 階中心差分近似の精度は O(h2 ) で、前進差分近似、後退差分近似よりも近似の精度が高いことがわかる。
最後に (3) と (4) を辺々加えて整理すると、次式を得る:

f ′′ (x) =(8)f (x + h) − 2f (x) + f (x − h)+ O(h2 )2h(h → 0).

この右辺第 1 項を 2 階中心差分商とよぶ。
問

上の各式を厳密に証明せよ。Landau の記号の定義に現れる定数 C はどう取ればよいか？

問

f が x の近傍で C 2 級ならば、lim

明せよ。

f (x + h) − 2f (x) + f (x − h)
= f ′′ (x) であることを証
h→0h2

余談 英語で o(·) は “littel-o notation” といい、それと対比させる場合 O(·) の方は “Big-O
notation” というらしい。英語で、大きい、小さいを表す形容詞はたくさんあるけれど、どれを使うのかは日本人には分かりにくい。

熱方程式の陽的差分法による解法

3

ここでは、第 6 章で取り扱った熱伝導方程式 (熱方程式) の初期値境界値問題

(10)

∂u
∂ 2u
(x, t) =(x, t) (t > 0, 0 < x < 1)∂t∂x2
u(0, t) = u(1, t) = 0 (t > 0)(11)u(x, 0) = f (x) (0 ≦ x ≦ 1)

(9)

を考える。
差分法では、解 u = u(x, t) の定義域 [0, 1] × [0, +∞) を図 1 のような格子に切り、各格子
点 (xi , tn ) における u の値

(12)

uni := u(xi , tn )の近似値 Uin を求めることを目標にする。

3

t

O

1 x

図 1: 閉領域 [0, 1] × [0, ∞) を格子で分割する
空間変数 x の区間 [0, 1] を N 等分して、等分点を xi (i = 0, · · · , N ) としよう。すなわち、刻み幅 (stepsize) を h = 1/N として、

xi = ih (i = 0, 1, · · · , N ).

(13)

とおく。次に時間変数 t に関する刻み幅 τ (> 0) を 1 つ定めて、

tn = nτ(14)(n = 0, 1, 2, · · · )とおく。
格子点 (xi , tn ) において、∂u∂ 2u
を前進差分近似、 2 を 2 階中心差分近似すると、
∂t
∂x

∂u
un+1
− uni
(xi , tn ) = i+ O(τ ) (τ → +0),∂tτ
uni+1 − 2uni + uni−1
∂ 2u
(x
,
t
)
=
+ O(h2 ) (h → +0)
i n
∂x2
h2

(15)
(16)

となる。u は熱伝導方程式 (9) を満たすので、h と τ が十分小さいとき、

un − 2uni + uni−1
un+1
− uni
i
≒ i+1
τ
h2

(17)
がなりたつ。
そこで(18)n
U n − 2Uin + Ui−1
Uin+1 − Uin
= i+1
τ
h2

という方程式を考え、この方程式の解 {Uin } を、{uni } の近似値 (もとの熱方程式の初期値境界値問題の近似解) に採用することにする。
このような、未知数列に関する方程式で、差分を含む方程式を差分方程式という2 。また差分方程式の解 Uin を差分解という。
差分解 = 差分方程式の解 = {uni } の近似値 = もとの方程式の近似解2「差分」とは何？という素朴な疑問に対して、きちんと回答するのは結構難しい。広辞苑によると「変数のとる一連の数値に対応して関数のとる値の差」とのことだが、例えば Fn = Fn−1 + Fn−2 のような Fibonacci 数
列を定義する「漸化式」も差分方程式と呼ばれる。

4

n
n
が現
, Ui+1
ここで添字 n は 0 以上の任意の整数に対して意味があるが、i については Ui−1
れるため i = 0, N に対して (18) を考えることはできず、 1 ≦ i ≦ N − 1 となることに注意しよう。分母を払って λ = τ /h2 とおくとnn).
− 2Uin + Ui−1
Uin+1 − Uin = λ(Ui+1

(19)
移項して

(20)

n
n
+ Ui+1
Uin+1 = (1 − 2λ)Uin + λ Ui−1(1 ≦ i ≦ N − 1; n = 0, 1, · · · )となる。
一方境界条件 (10) に対応して

U0n = UNn = 0 (n = 1, 2, · · · ),(21)
初期条件 (11) に対応して

Ui0 = f (xi ) (0 ≦ i ≦ N )(22)を考えるのは自然であろう。
以上をまとめると、熱伝導方程式の初期値境界値問題 (9), (10), (11) に対応して、未知数列{Uin ; 0 ≦ i ≦ N, n ≥ 0} に関する方程式系

(23)

n
n
+ Ui+1
Uin+1 = (1 − 2λ)Uin + λ Ui−1

(24)

U0n = UNn = 0(25)

Ui0 = f (xi )(1 ≦ i ≦ N − 1; n = 0, 1, · · · )(n = 1, 2, · · · )(0 ≦ i ≦ N )が得られる。
この方程式の解 {Uin ; 0 ≦ i ≦ N, n ≧ 0} は次のようにして計算することができる。

1. 右肩の添字 n が 0 の部分 {Uin ; 0 ≦ i ≦ N } については、(25) により値を求める。
2. n = 0, 1, · · · の順に次の手順で {Uin+1 ; 0 ≦ i ≦ N } を求める。
(a) 1 ≦ i ≦ N − 1 について Uin+1 を (23) で求める。
(b) i = 0, N については (24) により U0n+1 = UNn+1 = 0.
この方法は前進 Euler 法と呼ばれる。あるいは、差分方程式が Uin+1 について閉じた式 (それを求めるために方程式を解く必要がない) になっていることを前面に出して、陽解法または陽的差分法と呼ぶことも多い。

4

数値実験例
前節で導いた計算式に基づいて、実際にコンピューターで計算してみた結果を紹介しよう。

5

u(1/2, 0.1) の近似値
0.36654433
0.37118820
0.37232923
0.37261327
0.37268420
0.37270193
0.37270636

N
10
20
40
80
160
320
640

相対誤差 (%)
1.653709
0.407728
0.101583
0.025374
0.006342
0.001585
0.000396

表 1: (x, t) = (1/2, 0.1) で値の変化 (λ = 1/2).
10

1

Error

0.1

0.01

0.001

0.0001
10

100
N

1000

図 2: (x, t) = (1/2, 0.1) での相対誤差の変化 (λ = 1/2).

(1) 解析解との比較
偏微分方程式の問題に対して、差分法で求めた近似解、すなわち差分方程式の解のことを差分解とよび、それに対して元の偏微分方程式の解のことを厳密解とよぶ。厳密解は、それを解析的に表示した (式で表した) 場合解析解ともよばれる。
ここでは初期値境界値問題 (9)-(11) の解析解 u が簡単に求まる場合について、その解析解と差分解を比べてみよう。初期条件の関数 f を f (x) = f0 (x) ≡ sin πx とすると、解 u はu(x, t) = e−π t sin πx2(26)となる。このことは、第 6 章で扱った公式(27)

u(x, t) =

∞
X

bk e

−k2 π 2 t

Z 1
sin kπx,f (x) sin kπx dx (k = 1, 2, · · · )bk = 20k=1中の係数 bk を計算して b1 = 1, bk = 0 (k = 2, 3, · · · ) となることを示してもよいし、 (26)で定義される u が (9), (10), (11) を満足することを直接代入して確かめてもよい。表 1 は(x, t) = (1/2, 0.1) における値 u(1/2, 0.1) = exp(−0.1π 2 ) sin(π/2) = 0.372707838 · · · に対する近似値を並べてみたものである。これを見ると、誤差が N の増加とともに減少して行っているのがよくわかる。図 2 はその様子を両側対数グラフに表したものである (対数グラフの読み方については、補足の (2) にまとめておいた)。グラフの横軸が N で、縦軸が相対誤差を表す。このグラフから、誤差は N −2 に比例していると推察できる。

6

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5
t = 0.01

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5
t = 0.02

図 4: t = 0.01

図 5: t = 0.02

図 3: 初期条件 f = f1

(2) 解の漸近挙動
時刻 t を t → +∞ とした場合の u(·, t) の振る舞いのことを、解の漸近挙動とよぶ。
＜ディリクレ (Dirichlet) 境界条件の場合＞
初期値境界値問題 (9)-(11) の解の公式(28)

u(x, t) =

∞
X

bk e

−k2 π 2 t

Z 1
sin kπx,f (x) sin kπx dx (k = 1, 2, · · · )bk = 20

k=1

から、t → ∞ とすると、u は定数関数 0 に収束することがわかる:(29)

lim u(x, t) = 0 (0 ≦ x ≦ 1).

t→∞(式 (29) の厳密な証明は、補足 (3) にまとめておいた。)
このことを、初期条件の関数 f が



1



0≦x≦
 x
2
(30)
f1 (x) =



1

 1−x
≦x≦1

2
である場合 (グラフは図 3) に、差分法で (9)-(11) を解いて、確かめてみよう。N = 50, λ = 1/2
として解いた結果を描いたのが図 4 〜7 である。それぞれの図は、一つの固定された n の値に対する {Uin ; 0 ≦ i ≦ N } の値をもとにして描かれている。すなわち、t = nτ = 0.01, 0.02,
0.03, 0.04 としたときの x の関数 u(x, t) のグラフに対応するものである。グラフの横軸が位置、縦軸が温度を表していることを考えると、時間が経つにつれて、温度の高いところから低いところへ熱が流れていく様子がわかる。図 8 は t = 0 から t = 0.3 まで 0.01 刻みの時刻におけるグラフをまとめて描いてみたものである。時間が経つにつれて、温度が 0 に近付いていくのがわかる。

7

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5
t = 0.03

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5
t = 0.04

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.5
Tmax = 0.3, interval = 0.01

図 6: t = 0.03

図 7: t = 0.04

図 8: ディリクレ条件, t =0 ∼ 0.3, ∆t = 0.01

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5
t = 0.01

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5
t = 0.02

図 9: t = 0

図 10: t = 0.01

図 11: t = 0.02

＜ノイマン (Neumann) 境界条件の場合＞
ここでは境界条件を同次ディリクレ境界条件 (10) から、同次ノイマン境界条件(31)

ux (0, t) = ux (1, t) = 0 (t > 0)に変えた問題を考えよう。第 6 章と同様の計算により、
Z 1
∞
X
1
−k2 π 2 tu(x, t) = a0 +
(32)
ak e
cos kπx, ak = 2
f (x) cos kπx dx20(k = 0, 1, · · · )k=1

を導くことができる。この公式から、t → ∞ とすると、u はある定数関数に収束することが
わかる:

(33)

lim u(x, t) =

t→∞

a0
2

(0 ≦ x ≦ 1).

図 9〜13 では、初期条件は f1 のままで、境界条件を同次ノイマン境界条件にした問題を、
N = 50, λ = 1/2 の差分法で解いた場合の計算結果を示したものである。確かに時間が経つにつれて u(·, t) は定数関数に収束して行く様子がわかる。さらにその値についても調べてみよう。 f = f1 に対して、
Z 1
a0
1
=
f1 (x) dx = 図 3 のグラフと x 軸の囲む三角形の面積 =
(34)
2
4
0

8

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5
t = 0.03

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5
t = 0.04

heat equation, Neumann B.C.
theta=0, N=50, lambda=0.5
Tmax = 0.2, interval = 0.01

図 12: t = 0.04

図 13: t = 0.05

図 14: t = 0 ∼ 0.2

x
t
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40

0

0.1

0.2

0.3

0.4

0.5

0.000000
0.221849
0.246121
0.249466
0.249926
0.249990
0.249999
0.250000
0.250000

0.100000
0.227225
0.246862
0.249568
0.249940
0.249992
0.249999
0.250000
0.250000

0.200000
0.241301
0.248801
0.249835
0.249977
0.249997
0.250000
0.250000
0.250000

0.300000
0.258699
0.251199
0.2501650.250023
0.250003
0.250000
0.250000
0.250000

0.400000
0.272775
0.253138
0.250432
0.250060
0.250008
0.250001
0.250000
0.250000

0.500000
0.278151
0.253879
0.250534
0.250074
0.250010
0.250001
0.250000
0.250000

表 2: 定数定常解への収束 (ノイマン境界条件, θ = 0, N = 40, λ = 0.5)であるが、数値解の値は表 2 のようになっていて、解析解と同様に 1/4 = 0.25 へ収束して行く様に見える。 (問題が対称なので、区間の左半分の値のみを示してある。)(3) 熱伝導方程式の平滑化作用、t → ∞ の際の漸近形
既に何回か、「熱は温度の高いところから低いところへ流れる」ので、温度の高いところは低くなり、温度の低いところは高くなる、ことを見て来ているが、より詳しく分析することにしよう。
微分可能性とフーリエ係数の減衰
第 6 章で関数 f のフーリエ級数∞
a0 X+(ak cos kπx + bk sin kπx),2k=1の収束には、f の滑らかさ (微分可能性) が関係することに言及した。実は関数 f が滑らかなほど、フーリエ係数 ak , bk の絶対値 |ak |, |bk | は k が増加するにつれてより速く減衰することがわかっている。たとえば、具体的には次のような事実がある:(35)

f が C ℓ 級 (ℓ 回連続的微分可能)=⇒
9

ak , bk = o(k −ℓ )(k → +∞).

また、これとは逆に、f のフーリエ係数がより速く減衰するならば、f はより高い微分可能性をもつ (=より多い回数微分できる) こともわかっている。実際
∞
X

(36)

k ℓ (|ak | + |bk |) < ∞=⇒f は C ℓ 級.

k=1
かど

それでは、f が図 15 のような、たくさんの角のあるグラフを持つ関数 f2 の場合にはどうなるであろうか？この場合は、角となっている点で f は微分できないので、k → +∞ のときの |ak |, |bk | の減衰はゆっくりである3 。この f = f2 を初期値としたときの (9)-(11) の解を考えよう。解の公式 (27) を見ると、第 k 項には、時間に関係する因子 exp (−k 2 π 2 t) がかかっていて、これは t > 0 のとき、k の増加につれ急激に小さくなる。従って、初期条件 f = f2が低い微分可能性しか持っていなくても、少しでも時間が経過すると、u は高い微分可能性をもつ (実は無限回微分可能になる) ことが導かれる。
実際に f = f2 を初期データとした問題を差分法で解いた結果が図 15 〜 20 である。これによると、すでに t = 0.001 という早い段階で角が消滅していることがわかる。
凹凸の消滅
さらに時間が経過するにつれ、凹凸の激しいところほど急速に変形して消滅し、グラフが平滑な形に近づくことが見てとれる。凹凸が激しい (変化が急) なところがあるのは、大きな番
号 k に対する sin kπx の項が大きな割合を占めているからであると推測できる。しかし k が大きいほど、時間に関係する因子 exp (−k 2 π 2 t) は、t の増加につれて速く小さくなる。そのため、凹凸の変化が激しいところほど早く消滅する、と考えられる。
漸近形
公式 (27) において、t が非常に大きくなると、 k = 2, 3, · · · の項は、k = 1 に対応する項
b1 exp(−π 2 t) sin πx に比べて無視できるほど小さくなることがわかる。すなわち t が大きいところでは、u(x, t) ≒ b1 e−π t sin πx.
2(37)

これはグラフでいうと、形がサイン・カーブの半周期分に近づくことを意味する。実際、図 20
ではそうなっている。

陰解法

5(1) 安定性
前節で示した計算結果は、素朴に検討したに過ぎないが、なかなか満足の出来るものであった。ところが、うまく行かない場合もあることを見てみよう。
前節では、いずれの計算例も λ = τ /h2 を 1/2 あるいは λ = 1/4 として計算したのであるが、 λ = 0.51 として、f = f1 の場合の初期値境界値問題 (9)-(11) を解いたのが図 21〜24 である。解こうとしている熱伝導方程式の初期値境界値問題自体は図 4〜7 と同じものなのだが、計算結果は似ても似つかない悲惨なものになっている。このような現象を不安定現象とよぶ。
3

具体的には、ak , bk = o(1) (k → ∞) であるが、∞X

k (|ak | + |bk |) = ∞ である。

k=1

10

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25
t = 0.001

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25
t = 0.002

図 15: t = 0

図 16: t = 0.001

図 17: t = 0.002

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25
t = 0.003

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25
t = 0.004

heat equation, Dirichlet B.C.
theta=0, N=100, lambda=0.25
Tmax = 0.2, interval = 0.01

図 18: t = 0.003

図 19: t = 0.004

図 20: t = 0 ∼ 0.2, ∆t =0.01

これに対して、差分解の大きさが初期値、境界値で定まる定数 (N にも依存しない) で押えられるとき、差分スキームは安定である、という。
実は第 3 節で導いた差分方程式の解が、安定に計算できるための必要十分条件は、λ が0 < λ ≦ 1/2 を満たすことである、ということがわかっている。実際、0 < λ ≦ 1/2 のとき(38)|Uin | ≦ max |f (x)| (0 ≦ i ≦ N ; n = 0, 1, · · · )x∈[0,1]

が成り立つ (補足 (5) 定理 ?? 参照)。
さらに、この条件 0 < λ ≦ 1/2 が満たされる場合は、任意に固定した正数 T に対して

(39)

lim

N →∞

max |u(xi , tn ) − Uin | = 00≦i≦N0≦n≦T /τが成り立つこと、すなわち差分解の厳密解への収束も証明できる (これについても、補足 (5)定理 ?? を参照)。条件 0 ≤ i ≤ N , 0 ≤ n ≤ T /τ は、(xi , tn ) が [0, 1] × [0, T ] に含まれることを意味する。

(2) 後退 Euler 法
前項では、第 3 節で導いた差分方程式の安定な計算のためには 0 < λ ≦ 1/2 があれば十分であると述べたが、この条件を保証しようとすると、計算量が大きく (=計算に要する時間が長く) なりがちである。これを緩和するための、重要な手法を紹介しておこう。

11heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.51
t = 0.009996

図 21: t = 0.009996

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.51
t = 0.019992

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.51
t = 0.029988

図 22: t = 0.019992

図 23: t = 0.029988

heat equation, Dirichlet B.C.
theta=0, N=50, lambda=0.51
t = 0.039984

図 24: t = 0.039984
第 3 節では、時間に関する偏導関数
き換えて得られる(40)∂uを前進差分近似したが、これを後退差分近似に置∂tn
U n − 2Uin + Ui−1
Uin − Uin−1
= i+1
τ
h2

(1 ≤ i ≤ N − 1, n ≥ 1)を考えよう。番号 n を一つずらすと (n − 1 を n と書き直すと)、これは

(41)

n+1
n+1
Ui+1
− 2Uin+1 + Ui−1
Uin+1 − Uin=τh2(1 ≤ i ≤ N − 1, n ≥ 0)となる。移項して、二番目の添字が n であるものと、n + 1 であるものとを分離すると

n+1
n+1
(42)
(1 + 2λ)Uin+1 − λ Ui−1
+ Ui+1= Uin (1 ≦ i ≦ N − 1, n ≥ 0)のような差分方程式が得られる。この差分方程式にもとづいて差分解を定める方法を後退 Euler 法とよぶ。
(42) は第 3 節の方程式 (23) とは異なり、上付き添字が n + 1 の項が複数あることに注意しよう。(23) は、直接右辺の値を計算して、それを左辺に代入するという手順で、差分解 Uin+1(i = 1, 2, · · · , N − 1) が順次求まったが、 (42) を解くには、すべての i (1 ≦ i ≦ N − 1) についての方程式をまとめて考えた、連立 1 次方程式を解く必要がある。このことを区別するため、 (23) を陽 (的) 公式、(42) を 陰 (的) 公式とよぶ。
計算法の詳細は後述するが、後退差分近似公式で前項と同じ問題を、やはり λ = 0.51 に対して解いたのが図 25〜28 である。不安定現象が起こっていないのがわかる。実はこの公式 (42)では、λ (> 0) の値を何にしても、安定であることがわかっている。実際、図 29 〜32 では同じ問題を λ = 5 で解いてあるが、まったく不安定現象は起こっていない。

12heat equation, Dirichlet B.C.
theta=1, N=50, lambda=0.51
t = 0.009996

図 25: t = 0.009996

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=0.51
t = 0.019992

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=0.51
t = 0.029988

図 26: t = 0.019992

図 27: t = 0.029988

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=0.51
t = 0.039984

図 28: t = 0.039984
heat equation, Dirichlet B.C.
theta=1, N=50, lambda=5
t = 0.01

図 29: t = 0.01

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=5
t = 0.02

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=5
t = 0.03

図 30: t = 0.02

図 31: t = 0.03

heat equation, Dirichlet B.C.
theta=1, N=50, lambda=5
t = 0.04

図 32: t = 0.04

13

(3) θ 法, クランク-ニコルソン法
ここでは後退差分近似に基づく差分方程式の解の求め方を説明するが、せっかくだから少し一般化して、(41) の代わりに、(18) と (41) をいわば 1 − θ : θ の重みで重ね合わせた公式

(43)

n+1
n+1
n
n
Ui+1
− 2Uin+1 + Ui−1
− 2Uin + Ui−1
Ui+1
Uin+1 − Uin= (1 − θ)
+θ
τ
h2
h2

を考えよう。ここで θ はパラメーターで 0 ≦ θ ≦ 1 なる範囲から選ぶものとする。この式
は、θ = 0 とすると式 (18) に、θ = 1 とすると式 (41) に、それぞれ一致する。これまでと同様に、移項によって、左辺には添字が n + 1 の項が、右辺にはそれ以外の項があるようにすると


n+1
n+1
n
n
(44) (1 + 2θλ)Uin+1 − θλ Ui−1
+ Ui+1+ Ui+1= {1 − 2(1 − θ)λ} Uin + (1 − θ)λ Ui−1

(1 ≤ i ≤ N − 1, n ≥ 0)となる。これは ({Uin ; 0 ≦ i ≦ N } が既にわかっているとした場合) N +1 個の未知数 {Uin+1 ; 0 ≦i ≦ N } に関する N − 1 個の方程式である。境界条件より

U0n+1 = UNn+1 = 0(45)であるから、2 個の未知数 U0n+1 , UNn+1 は消去することができて、次のような N − 1 個の未知数 {Uin+1 ; 1 ≦ i ≦ N − 1} に関する N − 1 個の連立 1 次方程式を得る。



U1n+1
1 + 2θλ−θλ  n+1   U2
 −θλ
1
+
2θλ
−θλ






.
.
.
.

..
..
..
  ..   n+1 


 U

(46) 
−θλ 1 + 2θλ −θλ

 i



.
.
.
.

..
..
..
  .. 




n+1 


−θλ 1 + 2θλ−θλ   UN −2 

UNn+1
−1
−θλ
1 + 2θλ{1 − 2(1 − θ)λ}U1n + (1 − θ)λ(U0n + U2n ) {1 − 2(1 − θ)λ}U2n + (1 − θ)λ(U1n + U3n )




···


n
n
n

=
 {1 − 2(1 − θ)λ}Ui + (1 − θ)λ(Ui−1 + Ui+1 )  .


···




n
n
n
 {1 − 2(1 − θ)λ}UN −2 + (1 − θ)λ(UN −3 + UN −1 ) {1 − 2(1 − θ)λ}UNn −1 + (1 − θ)λ(UNn −2 + UNn )00

この方程式は後述するガウスの消去法などの方法によって解くことができる。
この公式を用いる方法を θ 法とよぶ。特に θ = 1/2 の場合は Crank-Nicolson (クランクニコルソン) の公式とよばれ、よく使われる (ルーツは記念碑的論文 Crank-Nicolson [?])。

6

熱伝導方程式以外の問題

ここまで熱伝導方程式のみ扱ってきたが、差分法の原理そのものは、より広い範囲の問題に適用できることは明らかであろう。熱伝導方程式を解くために第 6 章で学んだフーリエの変14数分離の方法は強力であるが、重ね合わせの原理が成立するような問題、すなわち線形同次の問題にしか適用できない。このうち同次性自体は、成立しなくても取り扱う方法がないでもな
いが、線形性の要請は本質的である。たとえば∂u∂ 2u(x, t) =(x, t) + ku(x, t)(1 − u(x, t)) (t > 0, 0 < x < 1)∂t
∂x2
ux (t, 0) = −1, u(t, 1) = 0 (t > 0)(47)u(0, x) = f (x) (0 ≦ x ≦ 1)のような非線形の問題をフーリエの方法で解くことは不可能であるが、差分法を用いて、とりあえず数値的に解くことはさして困難ではない。たとえば、n
U n − 2Uin + Ui−1
Uin+1 − Uin= (1 − θ) i+1
τ
h2
n+1
n+1
n+1
U
− 2Ui + Ui−1
+ kUin (1 − Uin )+ θ i+1h2(48)のような差分方程式を作ればよい。非線形偏微分方程式の研究にはコンピューター・シミュ
レーションの果たす役割が大きくなる所以である。

補足

7(1) ランダウ (Landau) の記号 O(·), o(·) について
この章の中で何回か、 O(h) (h → 0) とか o(n−k ) (n → ∞) のような記号を用いた。これはランダウの記号と呼ばれ、広く使われているが、念のために解説しておこう。
大きな O(·)

十分大きい正の定数 C を取ると、α に十分近い任意の x に対して

|f (x)| ≦ C|g(x)|(49)
が成り立つとき4 、このことを記号

f (x) = O(g(x)) (x → α)(50)

で表し5 、α において f が g で押さえられる、という。
小さな o(·)

条件

(51)

f (x)
=0
x→α g(x)limが成り立つことを、記号(52)f (x) = o(g(x)) (x → α)で表す。
D を f の定義域とするとき、論理式で書くと、(∃C > 0) (∃ε > 0) (∀x ∈ D : |x − α| < ε) |f (x)| ≤ C|g(x)|.
蛇足: 由緒正しい記号であるが、本当は f (x) ∈ O(g(x)) のように書く方が筋が通っている、という人もいる。O(g(x)) は 1 つのイデアルだ、とも。そのうちこの辺の記号は「現代化」されるかも知れない。
4

5

15

例

f が a において微分可能で、微分係数が A であるとは、
f (a + h) − f (a)
h→0
h

(53)

A = limが成り立つことであるが、この条件はf (a + h) − f (a) − Ah = o(h) (h → 0)(54)のように書き直すことが出来る。また関数 f が C k 級であるとき、テイラーの定理から(55)

f (a + h) = f (a) + f ′ (a)h +f ′′ (a) 2
f (k−1) (a) k−1 f (k) (a + θh) k
h + ··· +h+h2!
(k − 1)!
k!

となる θ (0 < θ < 1) が存在することが示されるが、この式から導かれる(56)

f (a + h) = f (a) + f ′ (a)h +f ′′ (a) 2
f (k−1) (a) k−1
h + ··· +h+ O(hk ) (h → 0)2!
(k − 1)!

はしばしば使われる。
注意

これらの記号は α が実数でなく無限大 ±∞ の場合にも使われる。例えば6x2 + x + 1= O(x) (x → ∞).
x(2) 対数グラフの読み方
第 4 節 (1) で、理工学の世界では常套手段である、変数の対数を目盛りに取ったグラフ —対数グラフ — を用いた。これについて簡単に説明しておこう。
二つの変数 x, y の間に y = ax + b (a, b は定数) のような 1 次の関係がある場合は、そのことを確認するのにはグラフを書いて直線になることをみればよく、グラフの傾きが a となる、というのは中学校以来の常識である。
それでは y = bxa (a, b は定数, b > 0) のような関係のある場合は、どうすればよいだろう
か？こういう場合に役立つのが対数グラフである。y = bxa の両辺の対数を取ると、(57)

log y = log bxa = a log x + log b.

そこで Y = log y, B = log b, X = log x とおくと、

Y = aX + B.
これは X, Y の間に 1 次の関係が成り立つことを示している。よって x, y の対数 X, Y をプロットすると、傾き a の直線になるはずである。
熱伝導方程式を差分法で解いた場合、分割数を N としたときの誤差 E は、近似的に E = C/N 2(C は定数) の形の関係を満たすので、この場合、横軸を log N , 縦軸を log E として、データをプロットすると、傾き −2 の直線上に乗ることになる。
6

主張を論理式で書くと、(∃C > 0)(∃R > 0)(∀x > R)x2 +x+1x≤ C ということになる。C = 3, R = 1 とす

れば、x > R を満たす任意の x に対して確かに不等式が成り立つ。

16

(3) 熱伝導方程式の解の指数関数的減衰
ディリクレ境界条件の場合の解 u について、(58)

lim u(x, t) = 0,t→∞

ノイマン境界条件の場合の解 u について、

a0
lim u(x, t) =
≡
t→∞
2

(59)

Z 1
f (x) dx0であると述べたが、証明をしておこう。どちらでも同様であるから、(??) だけ考える。解の減衰 (0 へ収束すること) は、公式

u(x, t) =

(60)

∞
X

bn e−n π t sin nπx
2 2

n=1

において、 lim e−n π t = 0 であることから、明らかのようであるが、
2 2

t→∞

(61)

lim

t→∞

∞
X

bn e

−n2 π 2 t

sin nπx =n=1のような無限和を取る操作

∞
X

∞
X
n=1

lim bn e−n π t sin nπx
2 2

t→∞

と、極限を取る操作 lim の順序の可換性は一般には成り立た
t→∞

n=1

ないことであるから、証明を要することである。ただ 0 に収束することを示すだけならば、級数が一様収束していることを証明すればよいが、せっかくだから、より詳しい結果を導こう。
関数 f (t) が、ある正数 C, τ に対し、

|f (t)| ≦ Ce−t/τ(63)(十分大きい任意の t)のように減少する指数関数 Ce−t/τ で上から評価されるとき、「f (t) は指数関数的に減衰する」という。
補題を認めれば t → ∞ の時、u(x, t) が 0 に収束するだけでなく、それが指数関数的に減衰していることが分かる。
補題の証明|f | は連続であるから、M := max |f (x)| が定まる。bn の定義式
x∈[0,1]Z 1(64)

bn = 2f (x) sin nπx dx
0

から

Z 1
(65)

|bn | ≦ 2

Z 1
|f (x) sin nπx|dx ≦ 2

0Z 1
M · 1 dx = 2M0dx = 2M.
0

したがって、|bn | は n によらない定数 2M で上から押さえられることが分かった。次に(66)

u(x, t) =∞X

bn e−n π t sin nπx
2 2

n=1

17

であるから|u(x, t)| ≦

(67)

∞
X

bn e

−n2 π 2 t

sin nπx ≤
!

n=1

=

2M

∞
X

2M · e−n2 π 2 t

· 1 = 2Mn=1∞
X

e−n π t2 2n=1e−π t .

e−(n −1)π t
2

∞ 
X

2

2

n=1

ゆえに t ≧ t0 とすると

|u(x, t)| ≦

(68)

2M

∞
X

!
e

−(n2 −1)π 2 t

0e−π t .
2

n=1

そこで
∞
X

(69)

e−(n −1)π t0
2

2

n=1

が収束することを示せばよいが、これはe−(n −1)π t0 ≦ e−(n−1)π t0 = rn−1 ,2(70)より、収束する等比級数

∞
X

2

2

r = e−π t0 < 12

rn−1 , 0 < r < 1 を優級数に持つことから明らかである。

n=1(余談: t0 = 0 とできないか？と疑問に感じる人に。例えば初期値 f が C 1 級の場合、
∞
X
2
C :=
|bn | とおくと、C < +∞ かつ |u(x, t)| ≤ Ce−π t (x ∈ [0, 1], t ≥ 0) が成り立つこと
n=1が分かる。)(4) 連立 1 次方程式 (46) が解けること
5 節 (3) の陰解法を実行するためには、連立 1 次方程式 (46) を解いて {Uin+1 } を求めなければならないが、これが可能であること (可解性, より強く一意可解性) を確かめよう。これはN − 1 個の未知数に関する N − 1 個の 1 次方程式であるから、係数行列が正則であることを示せばよい。
一般に、n 次の正方行列 A = (aij ) に対して
X
|aii | >|aij | (i = 1, 2, · · · , n)
(71)
1≦j≦n,j̸=iが成り立つとき、すなわち、行列の各行において、対角成分の絶対値が他の成分の絶対値の和よりも大きいとき、A は狭義優対角であるという。
5.3 節の連立 1 次方程式の係数行列は狭義優対角である。実際、最初の行と最後の行については明らかである。また、それ以外の各行については、対角成分は 1 + 2θλ で、他の成分で 0でないものは −θλ, −θλ の二つだけであり、(72)|1 + 2θλ| − | − θλ| − | − θλ| = 1 + 2θλ − 2(θλ) = 1 > 0となるから。それゆえ、次の定理を用いると、係数行列の正則性、ひいては連立 1 次方程式の可解性が分かる。

18

証明 背理法による。狭義優対角行列 A = (aij ) が正則でないと仮定すると、Ax = 0 を満たすベクトル x (6= 0) がある7 。

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

(73)

xn
とし、絶対値最大の成分を xk とする。方程式 Ax = 0 の第 k 成分を書くと
n
X

(74)

akj xj = 0.

j=1

この等式の左辺の第 k 項のみ残して、残りは右辺に移項すると
X(75)
akk xk = −akj xj .
j̸=k

絶対値を取って、(76)|akk ||xk | = −Xakj xj ≦X|akj ||xj | ≦

j̸=k

j̸=k

X

!
|akj | |xk |.

j̸=k

両辺を |xk | (> 0) で割算して(77)|akk | ≦X|ajk |.

j̸=k

これは行列 A が狭義優対角であるという仮定に矛盾する。

(5) 差分解の厳密解への収束の証明
ここでは熱伝導方程式の初期値境界値問題 (9)-(11) に差分法を適用することの正当性を保証する定理を与える。具体的には差分解の厳密解への収束と差分解の安定性を示す。
条件 (??) のことを、差分解の厳密解への収束とよぶ。
なお、差分解の安定性は、次の定理8 から分かる。

「任意の n 次の実正方行列 A と、A により定まる線型写像 f : Rn 3 x 7→ Ax ∈ Rn に対して、以下の条件は互いに同値である。(i) A は正則である (A−1 が存在する), (ii) f は全単射である, (iii) f は単射である,(iv) f は全射である, (v) ker f = {0}, (vi) rank f = n, (vii) det A 6= 0.」— この定理の証明の鍵は、次元定理「n − dim ker f = rank f 」である (これは準同型定理の系ですね)。
8
熱方程式の初期値境界値問題の解について、「最大値原理」と呼ばれる定理が成立するが、この定理はその離散化版と言うべきものになっている。
719

定理??の証明

第 2 節に述べたことから

(82)

un+1
− uni
i
= ut (xi , tn ) + O(τ ),τ(83)

uni−1 − 2uni + uni+1= uxx (xi , tn ) + O(h2 ).
h2

これから

un − 2uni + uni+1
un+1
− uni
i
− i−1
τ
h2
= ut (xi , tj ) − uxx (xi , tj ) + O(τ + h2 )(84)

= O(τ + h2 ).
ここで u は熱伝導方程式 ut = uxx を満たすことを用いた。同様にして

un+1 − 2un+1+ un+1
un+1
− uni
i
i+1
i
− i−1
= O(τ + h2 ).
τ
h2

(85)

以上から、{uni } は差分方程式を近似的に満たすことが分かる:

(86)

n+1
(1 + 2θλ)un+1− θλ(un+1i
i−1 + ui+1 )= [1 − 2(1 − θ)λ]uni + (1 − θ)λ(uni−1 + uni+1 ) + O(τ 2 + τ h2 ).
(この式は 1 ≦ i ≦ N − 1, 0 ≦ n ≦ J − 1 なる (i, n) に対して成り立つ。以下現れる方程式、不等式でも同様であるが、一々断らないことにする。) もちろん差分解 Uin は (その定義によって) 差分方程式
(87)

n+1
n+1
n
n
(1 + 2θλ)Uin+1 − θλ(Ui−1)+ Ui+1+ Ui+1) = [1 − 2(1 − θ)λ]Uin + (1 − θ)λ(Ui−1を厳密に満たすから、eni := Uin − uni とおくと、eni は

(88)

n+1
(1 + 2θλ)en+1− θλ(en+1i
i−1 + ei+1 )= [1 − 2(1 − θ)λ]eni + (1 − θ)λ(eni−1 + eni+1 ) + O(τ 2 + τ h2 )を満たす。パラメーター θ, λ に関する仮定から、係数 1 + 2θλ, θλ, 1 − 2(1 − θ)λ, (1 − θ)λ は
いずれも 0 以上になることに注意しよう。すると、

(89)

n+1
(1 + 2θλ)|en+1| − θλ(|en+1ii−1 | + |ei+1 |)≦ [1 − 2(1 − θ)λ]|eni | + (1 − θ)λ(|eni−1 | + |eni+1 |) + C(τ 2 + τ h2 ).
ここで C は十分大きな定数である。
いま E n := max |eni | とおくと、
0≦i≦N

n+1
(1 + 2θλ)|en+1| − 2θλE n+1 ≦ (1 + 2θλ)|en+1| − θλ(|en+1iii−1 | + |ei+1 |)(90)≦ [1 − 2(1 − θ)λ]|eni | + (1 − θ)λ(|eni−1 | + |eni+1 |) + C(τ 2 + τ h2 )≦ [1 − 2(1 − θ)λ]E n + 2(1 − θ)λE n + C(τ 2 + τ h2 )= E n + C(τ 2 + τ h2 ).
20

en+1
= en+1= 0 に注意して、i を 1 ≦ i ≦ N − 1 の範囲で動かして、左辺の最大値をとると
0N(1 + 2θλ)E n+1 − 2θλE n+1 ≦ E n + C(τ 2 + τ h2 ).

(91)
すなわち

E n+1 ≦ E n + C(τ 2 + τ h2 ).

(92)
この式を繰り返し用いると、

E n ≤ E n−1 + C(τ 2 + τ h2 )(93)

≤ E n−2 + 2C(τ 2 + τ h2 )···≤ E 0 + nC(τ 2 + τ h2 ) = nC(τ 2 + τ h2 ).
ここで、任意の i に対して e0i = Ui0 − u0i = f (xi ) − f (xi ) = 0 であるから E 0 = 0 となることを用いた。
ゆえに nτ ≦ Jτ ≦ T に注意すれば

E n ≦ CT (τ + h2 ).

(94)
これから

max |Uin − uni | ≦ CT (τ + h2 ).

(95)

0≦i≦N
0≦n≦J

定理 ?? の証明 (??) において、左辺はすべての Uin の最大値で、右辺は一部分の最大値であるから、左辺 ≥ 右辺は明らかである。以下、左辺 ≤ 右辺を示そう。このために、まず 0 ≦ n ≦ J −1
なるすべての n に対して


n+1
n+1
n+1
n
max Ui(96).
≤ max max Uk , U0 , UN
0≦i≦N

0≦k≦N

が成り立つことを示す。そのために

max Uin+1 ≤ max Ukn(97)

1≦i≦N −10≦k≦Nが成り立つことを示す。左辺の最大値が i = i0 で達成される、すなわちmax Uin+1 = Uin+101≦i≦N −1であると仮定しよう。このとき (??) を i = i0 に対して考える:(98)(1 + 2θλ)Uin+1− θλ(Uin+1+ Uin+1) = (1 − 2(1 − θ)λ)Uin0 + (1 − θ)λ(Uin0 −1 + Uin0 +1 ).
0
0 −1
0 +1

パラメーター θ, λ に関する仮定から、係数 1 + 2θλ, θλ, 1 − 2(1 − θ)λ, (1 − θ)λ はいずれも 0以上になることに注意しよう。すると、(99)(??) の右辺 ≤ {1 − 2(1 − θλ)} max Ukn + (1 − θ)λ( max Ukn + max Ukn )0≦k≦N{[1 − 2(1 − θλ)] + 2(1 − θ)λ} max Ukn
0≦k≦N
max Ukn .
0≦k≦N
0≦k≦N

=
=

21

0≦k≤N

一方

(??) の左辺 ≥ (1 + 2θλ)Uin+1− θλ{Uin+1+ Uin+1} = Uin+1
.
0
0
0
0

(100)

(??), (??) によって Uin+1≦ max Ukn . すなわち (??) が得られる:
00≦k≦N

max Uin+1 ≦ max Ukn .

(101)

1≦i≦N −10≦k≦N

これから (??) が成り立つことが分かる。(??) を n = 0 として適用すると、


1
1
0
1
max Ui ≤ max max Uk , U0 , UN .
(102)
0≦k≦N

0≦i≦N

さらに (??) を n = 1 として適用すると、


max Ui2 ≤ max
0≦i≦N

(103)

max Uk1 , U02 , UN2
0≦k≦N

.

(??), (??) から



(104)

max Ui2 ≤ max
0≦i≦Nmax Uk0 , U01 , U02 , UN1 , UN2
0≦k≦N




= max

max Uk0 , max U0ℓ , max UNℓ
1≦ℓ≦2
1≦ℓ≦2
0≦k≦N

以下この操作を繰り返すと、n = 1, 2, · · · , J に対して


ℓ
ℓ
0
n
(105)
max Ui ≤ max max Uk , max U0 , max UN
0≦i≦N

0≦k≦N

1≦ℓ≦n

1≦ℓ≦n

が示される。このことから (??) の左辺 ≤ 右辺が分かる。

(6) ガウス (Gauss) の消去法のアルゴリズム
連立 1 次方程式の解法として、線形代数の教科書には クラーメル (Cramer) の公式や掃ヨ ル ダ ン

き出し法 (Jordan の消去法ともいう) が説明されていることが多いが、ガウスの消去法は、掃き出し法を改良したものである9 。
例として次の方程式を取りあげて説明しよう。

 2x1 + 3x2 − x3 = 5(106)
4x1 + 4x2 − 3x3 = 3 −2x + 3x − x = 1.
12
3
掃き出し法では係数行列と右辺のベクトルを並べた行列を作り、それに1. ある行に 0 でない定数をかける。
2. 2 つの行を入れ換える。
3. ある行に別の行の定数倍を加える。
9

見方によっては、ガウスの消去法は中学校で習う加減法 (初めて習う解法！) そのものであり、大学の線形代数で習う解法は、実用性では退化していると言えなくもない (？)。

22

のような操作 — 行に関する基本変形と呼ぶ — をほどこして、連立方程式の係数行列に相当する部分を単位行列にするのであった。






2 3 −1 5
1 32 − 21 52
1 32 − 21 52






 4 4 −3 3  →  4 4 −3 3  →  0 −2 −1 −7  →−2 3 −1 1−2 3 −1 1
0 6 −2 6






5
11
1
0
−
−
1 32 − 12 52
1 0 − 54 − 11
4
4
4






7
1
7
→
→  0 1 21 72  →  0 1 12
0
1


→
2
2
2
5
15
0 6 −2 6
0 0 −2 − 2
0 0 1
3



  
1 0 0 1
x1
1



  
→  0 1 0 2  , ゆえに  x2  =  2  .
0 0 1 3x3
3
ガウスの消去法も、前半の段階はこの方法に似ていて、同様の変形を用いて掃き出しを行なうのだが、以下のように対角線の下側だけを 0 にする。






2 3 −1 5
2 3 −1 5
2 3 −1 5






 4 4 −3 3  →  0 −2 −1 −7  →  0 −2 −1 −7  .
−2 3 −1 1
0 6 −2 6
0 0 −5 −15
最後の行列は

2x1 + 3x2 − x3 = 5,−2x2 − x3 = −7,−5x3 = −15ということを表しているので、後の方から順に

x3 =

−15
= 3,
−5

x2 =

−7 + x3
= 2,
−2

x1 =

5 − 3x2 + x3
5−3×2+3==1
22と解くことが出来る。前半の対角線の下側を 0 にする掃き出しの操作を前進消去 (forwardelimination)、後半の代入により解の値を求める操作を後退代入 (backward substitution) と呼ぶ。
以下簡単に 3 つの方法の比較をしよう。
クラーメルの公式を適用するには n + 1 個の行列式を求める必要があるため、計算の手間がかかる (大きな計算量が必要になる、という)。実際、行列式を一つ計算するための手間は、連立方程式を一つ解くための手間と本質的に同等であることが分かっているので、クラーメルの公式を使うことに固執すると、本来必要な計算量の n 倍程度の計算をする羽目になり、大変な損をすることになる (差分法に応用する場合には、n が非常に大きな数になることに注意しよう)。そのため、実際の数値計算では、ごく特殊な例を除いて、クラーメルの公式が利用されることはない。クラーメルの公式は、理論的な問題を扱う場合に、真価が発揮されるものである。
このクラーメルの方法に比べれば、掃き出し法は、かなりの実用性を持っているが、第 5 節の方程式のように、係数行列が三重対角行列の場合には、ガウスの消去法の方が断然有利である。それは、ガウスの消去法を採用すると、掃き出しの途中に現れる行列が三重対角のままであることから、計算量が少なくてすむためである。
問 第 5 節の方程式で、n が 5, 6 程度の小さい場合に、簡単のため θλ = 1 として、掃き出し法とガウスの消去法を手計算で実行して計算の手間について比較せよ。

23

(7) (35) の証明(工事中)
任意の連続な周期 2π の周期関数 f に対して、
Z
Z
1 π
1 π
f (x) cos kx dx, bk (f ) =f (x) sin kx dx
ak (f ) =
π −π
π −π
とおく。
f が C ℓ 級であれば、部分積分によって、(
1
ak (f (ℓ) ) (ℓ が偶数)
ak (f ) = ± ℓkbk (f (ℓ) ) (ℓ が奇数),

1
bk (f ) = ± ℓk(bk (f (ℓ) ) (ℓ が偶数)
ak (f (ℓ) ) (ℓ が奇数)が導かれる。
f (ℓ) は連続であるから、Riemann-Lebesgue の定理により、

lim ak (f (ℓ) ) = lim bk (f (ℓ) ) = 0

k→∞

k→∞

が成り立つ。ゆえにlim k ℓ ak (f ) = lim k ℓ bk (f ) = 0.

k→∞

8

k→∞

文献案内
1. 非同次境界条件、Neumann 境界条件の場合の差分方程式→ 『熱方程式に対する差分法 I』 [?]2. 差分スキームの、行列を使った安定性の解析, von Neumann の安定性解析→ 『発展系の数値解析の続き』 [?], 『熱方程式に対する差分法 I』[?]
3. LU 分解の理論、Gauss の消去法の実際→ 『発展系の数値解析の続き』 [?], 『連立 1 次方程式 1』 [?]4. 波動方程式に対する差分法→ 『波動方程式に対する差分法』 [?]5. この文書を書くにあたって主に参考にしたもの
菊地・山本, 『微分方程式と計算機演習』 [?]6. 半線型放物型方程式に対する差分法について、特に勧めたいもの齊藤 宣一, 『線形・非線形放物型方程式の差分解法と可視化』 [?]7. 古典的 (?) 文献
Richtmyer-Morton [?], フォーサイス・ワソー [?], Smith [?], 山口・野木 [?]

24

参考文献
[1] 藤田宏：三訂版 応用数学, 放送大学出版協会 (2000), 第 7 章『発展系の数値解析』はhttp://nalab.mind.meiji.ac.jp/~mk/labo/text/heat-fdm-0.pdf から入手できる。
[2] Crank, J. and Nicolson, P.: A Practical Method for Numerical Evaluation of Solutions of Partial Diﬀerential Equations of the Heat-Conduction Type, Proc. Camb.
Phil. Soc., Vol. 43, pp. 50–67 (1947), http://www.springerlink.com/content/
g43649440307j807/ から reprint が入手可能.
[3] 桂田祐史：熱方程式に対する差分法 I — 区間における熱方程式 —, http://nalab.mind.
meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf (1998 年〜).
[4] 桂田祐史：発展系の数値解析の続き, http://nalab.mind.meiji.ac.jp/~mk/labo/text/
heat-fdm-0-add.pdf (1997 年〜).
[5] 桂 田 祐 史：連 立 1 次 方 程 式 I, http://nalab.mind.meiji.ac.jp/~mk/labo/text/
linear-eq-1.pdf (2002 年〜).
[6] 桂田祐史：波動方程式に対する差分法, http://nalab.mind.meiji.ac.jp/~mk/labo/
text/wave.pdf (1999 年?〜).
[7] 菊地文雄, 山本昌宏：微分方程式と計算機演習, 山海堂 (1991).
[8] 齊藤宣一：線形・非線形放物型方程式の差分解法と可視化, 明治大学大学院理工学研究科
で行われた集中講義の講義ノートで、 http://infsup.jp/ から入手可能 (2006).
[9] Richtmyer, R. D. and Morton, K. W.: Diﬀerence Methods for Initial-Value Problems,Krieger Pub Co, 2nd edition (1994).
[10] フォーサイス, ワソー：偏微分方程式の差分法による近似解法 (上), (下), 吉岡書店 (1968),
Finite-Diﬀerence Methods for Partial Diﬀerential Equations (1960) の翻訳 (藤野 精一訳).
[11] Smith, G. D.: Numerical solution of partial diﬀerential equations third edition, Clarendon
Press Oxford (1986), 第一版の邦訳が G. D. スミス著, 藤川洋一郎訳, コンピュータによる偏微分方程式の解法 新訂版, サイエンス社 (1996) である.
まさや

[12] 山口昌哉, 野木達夫：数値解析の基礎, 共立出版 (1969), 2009 年「復刊 数値解析の基礎」というタイトルで出版されました。

25