Introduction to Numerical Analysis of Differential Equations Naoya Enomoto (Kyoto.univ.Dept.Science(math))

微分方程式の数値解析入門

─

Introduction to Numerical Analysis of Differential Equations

─

Naoya Enomoto

∗

(Kyoto.univ.Dept.Science(math))

2001

はじめに

本稿は，微分方程式の数値解析手法を概説するものである．主に，差分法，有限要素法，境界要素法 の基礎的な事項を扱う．微分方程式の数値解析的な手法は，主に非線形微分方程式という一般解の振る 舞いが理論的に解析できていない対象について，その振る舞いを近似的に考察しようとするものである． 物理系の学生の多くが，理論系・実験系を問わず，どこかで数値シミュレーションとして上にあげたよう な解法を用いたことがあるだろう．しかし，それらはあくまでも近似解法でしかない．与えられた微分 方程式の厳密解を近似した解を得ているにすぎないのである．しからば，本来，数値的に得た解が厳密 解を本当に近似しているのかどうか，詳細に検討する必要がある．物理系と言えどもここはおろそかに できない．本稿では数値解析例を多く紹介するが，その中にも，微分方程式の厳密解ではないような解 を数値解として求めてしまう場合が数多く存在する．まずは，数値解析的な手法を用いる際，上のよう なことに留意せねばならないことを，理論・例を通して知って欲しい．一方，数値解析を円滑に行うた めには，安定したスキームで，かつ効率的に行えるかどうかもひとつのポイントとなる．これを数度の 実験を繰り返して得るという「実験的数値解析」ではなく，ある程度統一的な理論を重視してアプロー チする「理論的数値解析」の立場から，本稿では安定性と効率性について考察を加える． なお，数値解析例として紹介するものの中には，解析的数学理論が存在するものが多い．いわば「応 用解析的数理物理学」1が，様々な微分方程式について解析を行った結果があるのである．それらについ ても触れなければ，円滑に数値解析を行ったり，数値解析的手法の意味を理解することは難しい．そこ で，そうした理論も軽視せず，本文中に盛り込んでいく．多少難があるものは付録の形を取ることもあ る．また，筆者の紹介する数値解析例を鵜呑みにすることなく，読者は各自，数値解析を試みて欲しい． そうしなければ効果は半減するであろう．数値解析手法を手法として理解することと，実際にプログラ ムを書くこととは，それぞれに異なった困難があるのであり，手法として理解してもプログラムにおこ せるとは限らないのである． 本稿の概要について簡単に述べよう． まず，微分方程式というものがどういうものであるのか，具体的な方程式を「展示」しつつ，概説し， 数値解析手法についても簡単な説明を行う．その上で，続いて，差分法を重点的に概説する．ここでは， 差分近似式の作り方といった基本的なことを一通り述べた上で，熱伝導方程式を例にとり，差分法によ る数値解析を実際に行ってみる．それを見たあとで，さらに踏み込んで差分スキームの安定性や収束性 ∗_{e-mail:[email protected]} 1 数理物理学のもうひとつの側面として「代数解析的数理物理学」があると筆者は思う．それについては，拙稿『数理物理 学入門講義』を参照されたい．

といった理論的な事柄に言及する．常微分方程式・偏微分方程式を差分法によって数値解析する．扱う 方程式は，指数関数の方程式，単振動の運動方程式，ロジスティック方程式，Van der Pol 方程式，レス ラー方程式，ローレンツ方程式，1 階楕円型境界値問題といった常微分方程式，熱伝導方程式，1 階波動 方程式，2 階波動方程式，2 階楕円型境界値問題，非粘性および粘性 Bergers 方程式，数理生態学に現れ る非線形拡散方程式，非線形拡散方程式の爆発解，2 波相互作用方程式，KdV 方程式といった偏微分方 程式である．続いて，有限要素法・境界要素法についても，その基礎について触れ，いくつかの方程式 については，数値解析を行いたい．

目 次

1 微分方程式およびその数値解法について 5 1.1 微分方程式とは？ . . . . 5 1.2 なぜ数値解析か？ ─理論と応用の両面から─ . . . . 12 1.3 微分方程式の数値解法概観 ─基礎的 3 手法について─ . . . . 13 2 微分方程式の差分解法 15 2.1 差分解法とは何か？ . . . . 15 2.2 差分近似式のつくり方 . . . . 16 2.3 とりあえずやってみよう ─ 1 次元熱伝導方程式を例に─ . . . . 19 2.4 差分解の安定性と収束性 . . . . 22 2.4.1 数値的安定性の定義 . . . . 23 2.4.2 Fourier分解による安定性の判定 . . . . 24 2.4.3 Laxの同等定理 . . . . 26 2.4.4 Kreissの行列定理 . . . . 27 3 常微分方程式の差分法による数値解析 29 3.1 常微分方程式の解法についての一般的注意 . . . . 29 3.2 指数関数の常微分方程式 . . . . 29 3.3 単振動の運動方程式 . . . . 32 3.4 ロジスティック方程式 . . . . 34

3.5 van der Pol方程式 . . . . 38

3.6 ロトカ・ヴォルテラ方程式 . . . . 39 3.7 レスラー系・ローレンツ系 . . . . 40 3.8 1次元楕円型境界値問題 . . . . 41 4 偏微分方程式の差分法による数値解析 45 4.1 偏微分方程式についての一般的注意 . . . . 45 4.2 熱伝導方程式 . . . . 46 4.2.1 増幅率による安定性条件の導出 . . . . 47 4.2.2 安定/不安定の生じる理由は？ . . . . 47 4.2.3 陰差分スキーム . . . . 51 4.2.4 その他のスキーム . . . . 51 4.3 1階波動方程式 . . . . 53 4.3.1 1階波動方程式とは？ . . . . 53 4.3.2 主要な 3 つの差分スキーム . . . . 53 4.3.3 安定/不安定はなぜ生じるか？ . . . . 56 4.3.4 波束の拡散 . . . . 58 4.4 2階波動方程式 . . . . 59 4.5 楕円型境界値問題 . . . . 64 4.5.1 楕円型境界値問題の難しさ . . . . 64 4.5.2 2次元長方形領域に対する Poisson 方程式の境界値問題 . . . . 65 4.5.3 ラプラス作用素の差分化と楕円型境界値問題の安定性 . . . . 67 4.6 数理生態学に現れる拡散方程式 . . . . 70 4.6.1 ∆tに関して分岐現象を生じる例 . . . . 70

4.6.2 ∆xに関して分岐現象を生じる例 . . . . 73 4.7 粘性・非粘性 Burgers 方程式 . . . . 74 4.7.1 粘性・非粘性 Burgers 方程式とは？ . . . . 74 4.7.2 非線形性による波の崩壊 . . . . 75 4.7.3 衝撃波解 . . . . 78 4.8 爆発解を持つ非線形拡散方程式 . . . . 80 4.9 2波相互作用方程式 . . . . 82 4.10 KdV方程式 . . . . 83 4.10.1 ソリトン解の数値的発見 . . . . 83 4.10.2 ソリトン解を捉えよう． . . . . 83 5 微分方程式の有限要素法 86 6 微分方程式の有限要素法による数値解析 86 7 微分方程式の境界要素法 86 8 微分方程式の境界要素法による数値解析 86 9 付録 87 9.1 代数方程式および連立方程式の解法 . . . . 87 9.2 特性曲線法 ─ 1 階偏微分方程式の求積法の一例として─ . . . . 87 9.3 波動方程式の解の基本的性質 . . . . 87 9.4 熱伝導方程式の解の基本的性質 . . . . 87 9.5 Laplace・Poisson 方程式の解の基本的性質 . . . 87 9.6 数値的安定性の理論 ─ Lax の同等定理・Kreiss の行列定理─ . . . . 87 9.7 非粘性・粘性 Bergers 方程式 ─弱解としての衝撃波解を中心に─ . . . . 87 9.8 ある非線形拡散方程式における差分スキームの安定性と収束性 . . . . 87 9.9 ある非線形拡散方程式における爆発解について . . . . 87 9.10 KdV方程式のソリトン解について . . . . 87

1

微分方程式およびその数値解法について

概 要 本章では，微分方程式の数値解法全般について簡単に概説することを目標とする．そのためには， まず微分方程式がどういうものであるのかを理解しなければならない．そこで，まず第1節では，微 分方程式とは何かを定義し，その上で，微分方程式を「解く」ための諸条件について説明する．その 上で，大まかな分類を意識しながら，いくつかの微分方程式を選んで，「展示」する．続いて，第2節 では，なぜ微分方程式を数値的に解かなければならないのか，微分方程式の数値解析の必要性や意味 および特徴的な性格などを述べる．また，数値解法を行う際，どのような点に注意を怠るべきではな いのかについても述べる．さらに第3節では，微分方程式の具体的な数値解法として3つの手法，す なわち差分解法・有限要素法・境界要素法を取り上げて，その特徴を簡単に概観する．

1.1

微分方程式とは？

微分方程式とは何か．これに定義を与えることで答えるのは難しくない． 定義 1.1. 一変数関数 y(t) にある階数までの導関数 d i_y dti(i = 0, 1,· · · , p) の間に与えられた

関数関係を y に関する常微分方程式(ordinary differential equation) という．また，多変数関数

u(x1, x2,· · · , xn)のある階数までの偏導関数

∂ui

∂xj1∂xj2· · · ∂xji

(i = 0, 1,· · · , p) の間に与えられた関数 関係を u に関する関数関係を偏微分方程式(partial differential equation) という．

あまり定義のようなことにがちがち拘らなくても，微分方程式と言えば，例を思い浮かべることで，こ んなものだというイメージを持っている方も多いだろう．その意味でも，また数値解析という視点から 言えば，問題なのは，「微分方程式とは何か？」ということよりも，「微分方程式を 解く とはどういうこと か？」ということなのだとも言える．このことを考えるために，ひとまず微分方程式を大まかに分類し てみることを考えよう． 常微分方程式も偏微分方程式も，線形であるか非線形であるかによって，方程式の性質は異なってく るわけで，もちろんそのような分類は有効である．ただ，ここでは常微分方程式について，次のような 分類を与えてみよう． 常微分方程式の分類 ¶ ³ ：特殊関数を定義する常微分方程式． ：力学系を与える常微分方程式． µ ´ 常微分方程式を「解く」ことの意味は，教育のカリキュラム上でも，また数学理論の上でも変遷を遂げ てきた．例えば，古典力学の講義などで，常微分方程式を考えるときには，d 2_r dt2 =−ω 2_{xというような運} 動方程式を解くことが重要であった．つまり，一般解を具体的な表式として表すことが重要であったわ けである．しかし，これをひとたび数学として眺めてみると，そもそも解を具体的に書き下すことがで きる常微分方程式は，そもそもそれほどない（ほとんどないと言ってもいい）ことがわかってきた2．具 体的に解を書き下すことを求積 というが，求積できない常微分方程式の方がごく一般的であるというこ とがわかってきたわけである． そこで問題になるのが，そもそも「解は存在するのか」ということである．微分という概念は一体ど 2 そもそも具体的に書き下すというのが，なかなか定義に困る表現である．いわゆる初等関数の範囲でと考えておけばひと まず十分だろう．

のようなものであったかということを考えると，それは大域的な概念ではなく，局所的な（一次）近似 であった．その意味で，そもそも「局所的な解が存在するか」ということをまず考えなければならない であろう．その上で，そのような局所解の定義域がどこまで広げられるかという問題が出てくる．さら には，常微分方程式の場合，初期条件を与えて解の一意性が成り立つかというようなことが問題になっ てくる．そのような問題をひっくるめて「解の存在と一意性」の問題へと関心が移ってくる． 解の存在と一意性の問題 ¶ ³ (1)：局所的な範囲で解は存在するか？ (2)：局所的な解は，どこまで定義域を拡張できるか？ (3)：初期値問題で解の一意性は成立するか？ µ ´ この問題に対する結論をここで定理の形で述べることは避けよう3_{．ひとこと注意しておくとすれば，解} の存在しない常微分方程式，局所解には極大延長解なるものが存在するが，その中には，例えば任意の t ∈ R で解が存在するとは限らないようなものもあるし，初期値問題の解が一意ではない常微分方程式 さえあるということである． さて，求積可能な常微分方程式がそう多くは無いという事実を前にして，解の存在と一意性に関する 問題に関心が向いたわけであるが，もうひとつの方向性もある．すなわち，常微分方程式が，例えば初 等関数の範囲で具体的に書き下せないのなら，むしろ常微分方程式の解であるということを定義として， あらたな関数を定義してしまうという方向性である．これが，先の常微分方程式の分類における に相 当する．すなわち，常微分方程式の解関数として新たな特殊関数を定義しようというわけである．もち ろん，これには，ある常微分方程式が，本当に「新しい」特殊関数を定義しているのかという問題がつき まとうので，それに関しては厳密な証明が必要である．しかし，現在までに多くの常微分方程式をもと に新たな特殊関数が定義されることがわかっている．また，そのような特殊関数が物理においては，日 常的に用いられている．例えば，以下のようなものがあげられるだろう． 特殊関数を定義する常微分方程式 ¶ ³ 指数関数 : dx dt = x 三角関数 : d 2_x dt2 =−x リッカチ関数 : dx dt = x 2 + p(t)x + q(t) ベッセル関数 : d 2_x dt2 + 1 t dx dt + ( 1− n 2 t2 ) x = 0 エルミート関数 : d 2_x dt2 + 2t dx dt + 2nx = 0 ルジャンドル関数 : (1− t2)d 2_x dt2 + 2t dx dt + n(1 + n)x = 0 超幾何関数 : x(1− x)d 2_y dx2 + (γ− (α + β + 1)x) dy dx − αβy = 0 µ ´ これらの常微分方程式は，古典的なもので，性質もよく調べられている．いわば現代の特殊関数の立場 にたちつつあるのが，Painlev´e 関数である．それは次のような常微分方程式で記述される． 3 詳細は，[高橋] あるいは [伊藤] などを参照されたい．

Painlev´e方程式 ¶ ³ P : d 2_y dt2 = 6y 2_{+ t} P : d 2_y dt2 = 2y 3_{+ ty + α} P : d 2_y dt2 = 1 y ( dy dt )2 −1 t dy dt + 1 t(αy 2_{+ β) + γy}3₊ δ y P : d 2_y dt2 = 1 2y ( dy dt )2 + 3 2y 3_{+ 4ty}2_{+ 2(t}2_{− α)y +} β y P : d 2_y dt2 = ( 1 2y + 1 y− 1 ) ( dy dt )2 − 1 t dy dt + (y− 1)2 t2 ( αy + β y ) P : d 2_y dt2 = 1 2 ( 1 y + 1 y− 1 + 1 y− t ) ( dy dt )2 − ( 1 t + 1 t− 1 + 1 y− t ) dy dt +y(y− 1)(y − t) t2_(t− 1)2 ( α + β t y2 + γ t− 1 (y− 1)2 + δ t(t− 1) (y− t)2 ) µ ´ この方程式は，いわば「代数解析的数理物理学」と密接に関わる方程式であると言えよう． しかし，こうした特殊関数を定義する常微分方程式は，数値解析の範疇には入りにくい．まさにその 関数の（荒っぽく言うのであれば，代数的・幾何的）性質が重要であり，数値解析として近似的に数値 的に解を求める意味が薄いのである．このような関数が物理において重要な役割を果たしていることは 確かであるが，例えば表現論などの数学の分野が，こうした重要性を基礎付けてくれるであろう．この 数値解析入門では，こうした方程式は扱わない． では，数値解析的な興味の対象にあたる常微分方程式はどのようなものであろうか．先に，解の存在 と一意性が常微分方程式の関心事になったと書いたが，やがて時代が進むと，また個々の解関数の振る 舞いに対する関心が生じてきた．様々な非線形常微分方程式の解の定性的振る舞いを論じることで，物 理現象との接点となる事柄が出てきたわけである．その代表格がいわばカオス理論のような力学系の話 題であると言えるだろう．例えば，解の存在と一意性がわかっていても，解を何らかの形で具体的に書 き下すことができないということは重要な制約である．常微分方程式の形から，解関数の定性的性質を 調べなければならないのである．ここに，数値解析の果たすべき役割が見えてくる． 少しこの観点から常微分方程式を挙げてみよう． まず，いわゆる Newton の運動方程式のようなものが挙げられる． Newtonの運動方程式 ¶ ³ N ewtonの運動方程式 : md 2_r dt2 = F (r) 単振動の運動方程式 : d 2_x dt2 =−ω 2 x µ ´ 次に，電気的な振動を記述するような方程式もある．

リエナール方程式 ¶ ³ リエナール方程式 : d 2_x dt2 + f (x) dx dt + xg(x) = 0

van der Pol方程式 : dx

2 dt2 − a(1 − x 2 )dx dt + x = 0 (a > 0) ダフィング方程式 : d 2_x dt2 + a dx dt + x− x 3 _{= 0} µ ´ 力学におけるパラメータ共振の方程式としては，Hill の方程式とよばれるものがある． ヒルの方程式 ¶ ³ ヒルの方程式 : d 2_u dt2 + q(t)u = 0 マシューの方程式 : d 2_u dt2 + (ω 2_{− 2² cos t)u = 0} µ ´ 数理生態学でも興味深い常微分方程式がある． 数理生態学の方程式 ¶ ³ ロトカ・ヴォルテラ方程式 :      dx dt = αx− γxy dy dt =−βy + δxy µ ´ カオスと関連した方程式を挙げてみよう． カオスに関係する方程式 ¶ ³ ロジスティック方程式 : dN dt = αN (1− λN) レスラー方程式 :            dx dt =−y − z dy dt = x + αy dz dt = b− cz + xz ローレンツ方程式 :            dx dt = σ(x− y) dy dt = γx− y − xz dz dt =−bz + xy µ ´ このような方程式について数値解析的な方法を用いることで，数値解析の「できること」を実感するの がよいと思う．それを以下の章では試みることにしたいのである． 一方，偏微分方程式はどうであろうか．そもそも偏微分方程式には，包括的な形での「解の存在と一

意性を保証する定理」はあまりない．線形偏微分方程式の 3 タイプ，放物・双曲・楕円については，それ なりの解析的理論が存在するが，非線形偏微分方程式については，まだ包括的な解析的理論はないよう である．しかし，物理現象から近似という操作を経て導出される微分方程式の多くは偏微分方程式であ る．このあたりに数値解析の果たすべき姿が見えてくる． 2階線形偏微分方程式の 3 つの型，放物・双曲・楕円を中心に，偏微分方程式をいくつかみてみること にしよう． 放物型方程式の代表的な例が熱伝導方程式であるといえる． 放物型偏微分方程式 ¶ ³ 熱伝導方程式（or 拡散方程式） : ∂u ∂t = ν4u [熱源ありの場合] : ∂u ∂t = µ4u + f(x, t) 移流拡散方程式a : ∂u ∂t = ν ∂2u ∂x2 − β ∂u ∂x 反応拡散方程式 : ∂u ∂t = ν ∂2_u ∂x2 + ku 2 生態学に現れる反応拡散方程式b _: ∂u ∂t = ν ∂2_u ∂x2 + au− ku 2 a_{ドリフトを伴う Brown 運動の方程式としても導出される．} b_{1930年に Fisher が優勢種のはびこりを記述する方程式として提出したもので，Kolmogorov-Petrovskii-Piskov によ} る進行波解の研究もなされ，燃焼面の進行を記述する方程式としても知られている． µ ´ 反応拡散方程式の一般形は，u をベクトルとして，連立偏微分方程式形 ∂u ∂t =4u + f(u) で与えられる．これは，ある場合においてフラクタル的なパターン形成が見られることすら知られてい る． 次に，波動方程式に代表される双曲型方程式について見てみることにしよう．

双曲型方程式 ¶ ³ 波動方程式 : ∂ 2_u ∂t2 =4u [外力項のある場合] : ∂ 2_u ∂t2 =4u + f(x, t) 電信方程式a [速度に比例する抵抗のある場合の波動方程式] : ∂ 2_u ∂t2 + k ∂u ∂t − c 2∂2u ∂x2 1階波動方程式 : ∂u ∂t + c ∂u ∂x = 0 非粘性バーガーズ方程式 [1階波動方程式の非線形化] : ∂u ∂t + cu ∂u ∂x = 0 バーガーズ方程式b _: ∂u ∂t + cu ∂u ∂x − ν ∂2_u ∂x2 = 0 KdV 方程式 : ∂u ∂t + 6u ∂u ∂x + ∂3_u ∂x3 = 0 a_{Heavisideが直線状ケーブルの電磁波の減衰の様子を記述する方程式として用いたことによる．} b_{Navie-Stokes方程式の圧力項を圧縮したものとみることもできる．} µ ´ このように波動方程式の仲間たちは，数多く存在するわけであるが，多くは物理現象との関わりで興味 深い性質を持っている．しかも数値解析的にもなかなか難しい問題をはらんでいるものが多い． Laplace・Poisson 方程式というような楕円型方程式についてみてみよう． 楕円型方程式 ¶ ³ Laplace方程式a _{: 4u = 0} P oisson方程式 [Laplace方程式の非斉次化]   : 4u = f Helmholtz方程式b : 4u + k2u = 0 変形 Helmholtz 方程式 : 4u − k2u = 0 a_Laplace 方程式が支配する現象は，定常現象に数多く見られる．例えば，重力ポテンシャル・静電ポテンシャル・定 常熱伝導・定常物質拡散・非圧縮性渦なし流れなど． b_{時間によらない Schrodinger 方程式} ~ 2 2m4φ + Eφ = 0 はこの形になっている． µ ´ 放物・双曲型が非定常現象を記述する偏微分方程式であったのに対し，楕円型方程式は定常現象を記述 する偏微分方程式であることに注意しておこう． さて，高階偏微分方程式を見てみよう．

高階偏微分方程式

¶ ³

均質等方的弾性体の運動方程式 : ρ∂

2_u

∂t2 = µ4u + (λ + µ) grad div u + ρf

均質等方的弾性体の平衡方程式 : −µ4u − (λ + µ) grad div u = ρf 薄板の曲げに対する運動方程式 : ρ∂ 2_ν ∂t2 =−(λ + µ)4 2_{ν + ρg} 棒の微小横振動を表す方程式 : ρ∂ 2_u ∂t2 =− ∂2 ∂x2 ( EI∂ 2_u ∂x2 ) + f µ ´ 最後に，物理学における重要な偏微分方程式を見て置くことにしよう．このあたりの数値解析は，もち ろんこの入門では扱えないし，そもそも扱えるのかどうかすら筆者はしらないが，いずれは扱えるよう になるのであろうか・・・． 数理物理学における偏微分方程式 ¶ ³ Schrodinger方程式 : i~∂ ∂tφ = ( − ~2 2m4 + V (r) ) φ M axwell方程式 :                div E = ρ ε div B = 0 rot E + ∂B ∂t = 0 rot B− εµ0 ∂E ∂t = µ0i Eular方程式 : ∂v ∂t + (v· ∇)v + 1 ρgrad p = f N avie− Stokes 方程式 : ∂v ∂t + (v· ∇)v + 1 ρgrad p− ν4v = f Eular− Lagrange 方程式 : d dt ∂L ∂ ˙qi − ∂L ∂qi = 0 Hamilton方程式 : { ˙ qi = ∂H_∂p i ˙ pi =−∂H_∂qi Hamilton− Jacobi 方程式 : H (q, ∇qS, t) + ∂S ∂t = 0 Klien− Gordon 方程式 : ( 4 − 1 c2 ∂2 ∂t2 − m2_c2 ~2 ) φ = 0 Dirac方程式 : 4 ∑ µ=1 γµ ∂φ ∂xµ + iκφ = 0 Einsteinの重力方程式 : Ri_j(x)− 1 2δ i jR(x) = κT i j (1≤ i, j ≤ 4) [宇宙項を入れたもの] : Ri_j(x)−1 2δ i jR(x) + Λδ i j = κT i j (1≤ i, j ≤ 4) µ ´ さて，こうした微分方程式を解くためには，いくつかの条件が必要である．それについて若干注意 しておこう．なお，これについては，それぞれの方程式によって多様であり，また数値解析的な視点か ら言ってもいろいろ複雑な問題があるので，ごくごく一般的な注意をするだけにとどめ，各論のところ

で必要なことは補うようにしたい． 一般に，常微分方程式を解く際には，しかるべき初期条件が必要である．例えば，2 階の常微分方程式 であるような単振動の方程式 d2_x dt2 =−ω 2_x を考えると，初期条件として，時刻 t0における位置 x(t0)と，そのときの速度 dx dt(t0)が与えられていな ければならない．一般に n 階常微分方程式を解くためには，n 個の初期条件が必要である．これは，そ れほど複雑な話ではない．しかし，偏微分方程式になるとかなり話はややこしくなる． 偏微分方程式を解くための諸条件 ¶ ³ (1):初期条件（ある時刻における未知関数あるいはその導関数を与える条件） (2):境界条件（考えている領域の境界での未知関数あるいはその導関数を与える条件） (3):適合条件（初期条件と境界条件が矛盾しないための条件） µ ´ 初期条件を与えて偏微分方程式を解くことを初期値問題，初期条件＋境界条件によって解くことを 初期・境界値問題という． 境界条件については，いくつかの形がある． 偏微分方程式を解くための諸条件 ¶ ³ (1):Derichlet境界条件（境界に置ける未知関数そのものを指定する境界条件） (2):Neumann境界条件（境界での未知関数の導関数を指定する境界条件） (3):周期境界条件（境界の両端での未知関数が一致することを指定する境界条件） µ ´ もちろん境界条件には，ある領域の「境界」というものの数学定式化の難しさによる困難があり，その 数学的な扱いは，考えるほど容易にはいかないものがある．しかし，ここではそのことには触れない． これまでのことを踏まえて，数値解析の果たすべき役割について，節を改めて考えてみよう．

1.2

なぜ数値解析か？ ─理論と応用の両面から─

前節で述べてきたことを敷衍して，なぜ数値解析なのかという問題を考察しておくことにしよう． 前節でみたことを復習しておこう．微分方程式というものを考えたとき，常微分方程式・偏微分方程 式を問わず，まず解を具体的に表示する，すなわち求積することが考えられた．しかし，実際に多くの 微分方程式を扱ってみると，求積可能な微分方程式はむしろ少数派であって，解を具体的に表示できな い微分方程式が多数派であることがわかってきた．ここで，例えば特殊関数を微分方程式によって定義 するという方向性が生まれ，現代では Painlev´e 方程式などのようなものが興味深い対象として，代数解 析的数理物理学を中心に扱われている．一方で，解の存在と一意性の問題が，解析学としての微分方程 式論では重要になった．それについては，常微分方程式・偏微分方程式を問わず様々な理論がつくられ てきたが，その後は，解を具体的に表示できないとしても，何らかの解析的方法によって，その定性的 性質を調べるという方向が重要になってきたというわけである． 上のような背景のもとで微分方程式を考えてみると，線形の微分方程式については，常微分方程式・ 偏微分方程式を問わず比較的綺麗な理論があるが，非線形微分方程式については，なかなか包括的な理 論がないばかりか，個々の方程式についても有効な方法がなかなかないのが現状であろう．そこで数値

解析という手法が重要性を帯びて来るのである． 数値解析学を文章として定式化するとすれば， 「厳密解を近似する数値解を構成し，その挙動を調べることで，厳密解の挙動を予測すること」 だともいえるだろう．非線形微分方程式について，解析的な方法にまだ限界がある現状では，このよう な方法を試みることは重要であり，計算機の高度な発達も，このような方法を可能にしてくれているの である．だだし，ここで注意しておくべきことの一つは，ただ単に微分方程式をなんらかの形で数値的 に「解く」ことだけが，数値解析では ない ということである．すなわち，数値解析は，むしろ数値的に シミュレーションすることで，もとの微分方程式の解析的理論をつくることも含んでいるのであり，た だ単なる計算手段ではないということである． また，注意すべきことの 2 点目は，数値解析には数値解析なりの弱点もあり，それに対する注意を怠れ ば，とんでもない結果を得ることになりかねないということである．これがどのようなことなのか，一 般的な観点から述べておこう．具体的なことについては，あとの章でじっくりと見ることになる． 数値解析とは言っても，またいかに高性能の計算機を用いるとはいっても，計算を行える桁数はあく までも有限桁であって，そこには切り捨て・切り上げという操作が行われる．これによって，数値解を 計算する際には，「丸め誤差」と呼ばれる数値的な誤差が不可避である．したがって，何度も計算を繰 り返すうちに，この誤差が累積するようなことになれば，数値計算は意味をなさないことになる．これ を数値的安定性 の問題という．また，同時に，いかに高性能の計算機を用いるとはいえ，その数値計算 の効率性を度外視して進めるわけにはいかない．一度の計算にあまりに時間がかかるような数値計算は 避けなければならない．したがって，数値計算を行う際には，効率の良い数値計算を行うと同時に精度 が良く，しかも数値的安定性を持つようなプログラムを考えなければならないのである． さらに注意すべきこととして，数値解の収束性 の問題がある．あくまでも数値解は近似解であって， 厳密解ではない．したがって，数値解が厳密解に収束しているかが重要なのである．ある数値解を数値 的にシミュレーションしてみて得られた結果が本当にもとの微分方程式の厳密解の挙動を表しているの かどうかについて，数学的な証明が必要なのである．これを軽視していると，かなりとんでもない解を 厳密解の性質を反映した解であると判定しかねない．そのような例についてもあとで見る．このことは， 数値解析がただ単なる数値シミュレーションを行うことではないということと密接に関係している．数 値解が厳密解の性質を性格に反映していることを証明するには数学の解析的理論が必要なのである． 以上に述べたことを，本稿の目的に即してまとめておこう． 数値解析学に関して本稿で強調する点 ¶ ³ (1):数値解析とは言っても，数学における解析的理論が重要であること． (2):数値解析とは言っても，微分方程式に関する一般論が重要であること． (3):差分法を中心にして，数値解析的手法のメカニズムを説明すること． (4):いくつかの微分方程式について差分法による数値シミュレーションを試みること． µ ´ さて，次節では，数値解析的手法の特徴などを見てみることにしよう．

1.3

微分方程式の数値解法概観 ─基礎的

3

手法について─

本節では，差分法，有限要素法，境界要素法について，その特徴や用いられる問題などについて概 観する．まず，3 つの方法を大まかに説明してみよう．

ある領域上で微分方程式を解こうとするとき，その解である関数は，領域上のすべての点で値を決め なければならないが，[0, 1] 区間には，（非可算）無限個の点があるわけで，その点での値を，人の手で計 算して決めることはできない．数値計算行うということは， 「領域内の ´有 ´限 ´個 ´の ´点での値を決めることによって，解である関数を近似する操作」 であると言える．差分法も有限要素法も境界要素法も，その手法として上に述べたような操作である． 差分法 とは， 考えている領域を ´正 ´方 ´形（ ´な ´い ´し ´は ´長 ´方 ´形） ´の ´メ´ッ´シ´ュで覆い， 各格子点上での解の値を数値的に計算するメカニズム である．それに対して，有限要素法 は， 考えている領域を ´三 ´角 ´形（ ´な ´い ´し ´は ´多 ´角 ´形） ´の ´メ´ッ´シ´ュで覆い， 各格子点上での解の値を数値的に計算するメカニズム であると言える．境界要素法 は，領域の境界で有限要素法を行うものである． 方法の詳細は各章で説明するので，ここでは，3 つの方法を表にして比較してみることにしよう4_． 差分法  有限要素法 境界要素法 問題の定式化 差分近似 弱形式表現 （汎関数の極値問題） 積分方程式表現 離散化 格子点 有限要素 境界要素 解法 領域型 領域型 境界型 連立一次方程式の 係数行列 （対称）疎・帯性 （対称）疎・帯性 （対称）疎・帯性を示さず 特徴 ・離散化が簡便 ・計算時間が早い ・離散化の精度やスキ ームの安定性が明確 ・任意の形状に対す る取り扱いが容易 ・境界条件（自然境 界条件）の処理が 容易 ・計算の自動化が達 成し易い ・問題の次元が 1 次 元低下 ・外部問題と内部問 題を同時に扱える ・境界の微分量も未 知量として定式化 できる 適用性 ・非線形問題も含め て極めて広範 ・特に流体解析 ・非線形問題も含め て極めて広範 ・主に構造解析（最 近は流体解析にも） ・基本解の存在する 問題に限定 ・線形 3 次元問題、 特異応力問題、 移動境界問題 続く第 2 章では，特に差分解法についてじっくりと見ることにする． 4 以下の表は，[大西・登坂] による．

2

微分方程式の差分解法

概 要 本章では，微分方程式の差分解法について，理論的な側面を中心に述べる．まず第1節では，差分 解法とはどのようなものであるかを説明する．次に第2節で，差分近似式のつくり方を中心に述べた あと，第3節で，熱伝導方程式を題材に取り，ひとまず差分解法を行ってみることで，差分解法を扱 う上でいかなる問題点があるのかについて，注意を促す．その上で第4節において，差分解の安定性 と収束性について，Laxの同等定理とFourier分解の方法と行列の固有値評価法を中心に述べる．

2.1

差分解法とは何か？

本節では，差分解法がどういうものであるかを，より詳しく定式化しておく．すでに，差分法 とは， 考えている領域を ´正 ´方 ´形（ ´な ´い ´し ´は ´長 ´方 ´形） ´の ´メ´ッ´シ´ュで覆い， 各格子点上での解の値を数値的に計算するメカニズム である，と前節で定式化していたわけだが，ここでは，それをより詳細な言葉で説明し，どういう Step を踏んで近似解を得るのかを述べるわけである． ある微分方程式を差分解法によって解くことを考える時，まずしなければならないことは，差分近似 式によって，もとの微分方程式を離散化し，差分方程式を得ることが重要である．差分近似式とは， ∂f ∂x ∼ f (x + ∆x)− f(x) ∆x のように，微分を，差分の形で置き換えたものを言う．もちろん，近似式のつくり方は，多種多様である し，また 2 階以上の微分についてもいろいろなやり方が考えられる．これについては次節で詳細に論じ るし，差分近似式のつくり方の難しさは，あとで見るような種々の例でわかるであろう．ここでは，何 らかの方法によって，微分方程式を差分方程式に離散化するのである． さて，次に行うことは，離散化した際の ∆x, ∆t の間隔を決定し，「時間柱」を格子上のメッシュで覆う． ここで言う時間柱とは，考えている領域 D とそれに時間方向を付け加えたものである（下の図を参照)． D t Dが 2 次元平面内の領域の場合 Dが 1 次元の区間の場合 t ただし，あとで見るように，楕円型などの定常問題を扱う場合の差分法は，時間的に定常なものを扱っ ているので，時間軸を考える必要はなく，領域 D だけを考えれば良い． メッシュで覆ったならば，次に行うことは，初期条件や境界条件など，与えられた諸条件によって値が 決定している格子点をチェックする．このような条件をすべてまとめたものを差分スキーム という．そ の上で差分スキームを解けば良い．この解く段階では，関数の値が逐次決まっていく場合もあるが，連 立一次方程式や代数方程式を解かなければならない場合もある．前者を陽差分法，後者を陰差分法 と呼 ぶ．「逐次決まる」とか「連立一次方程式や代数方程式を解く」ということの意味やなぜそのような違い

が生じるのかは，あとで，具体例を見てもらうことにする． さて，これで数値解が差分近似式を解くことによって構成されたわけである．もちろん，ただ単なる 数値シミュレーションだけが目的ならば，これで終わりであるが，数値解析と言う以上は，差分法によっ て得られた結果の吟味が必要不可欠である．ここで安定性や収束性の問題が出てくる．これについては， 概括的なことを述べるだけでは難しいので，あとで熱伝導方程式を例にとって，少し具体的な状況がわ かるようにしてから述べることにする． 今までのことをまとめておこう． 微分方程式の差分解法 とは，与えられた微分方程式を差分近似式によって離散化し，与えられた 諸条件を加えて差分スキームを構成し，それに基づいて時間柱を覆う格子点での値を決定すること によって，与えられた微分方程式の近似解を構成する方法を言う． 差分法は，以下のようなメカニズムに従って行われる． (1):与えられた微分方程式の微分項を差分近似式によって離散化する． (2):得られた差分近似式とあらかじめ与えられた初期条件や境界条件とを合わせて，差分スキーム を構成する． (3):得られた差分スキームをもとに数値計算を行って，格子点での値を計算する． (4):場合によっては、差分法の安定性や収束性を吟味し，再度数値計算を試みるなど，検証を行う．

2.2

差分近似式のつくり方

本節では，差分近似式のつくり方を詳細に見ることにしよう．差分近似式がつくれなければ，そもそ も差分法を行うことはできないし，また差分近似式としてどのようなものを選択するかは，それほど簡 単な問題ではなく，やや慎重に行う必要もある． さて，差分近似式のつくり方を説明するために，まず，関数の微分（導関数）の定義を復習しておこ う．それは，次のような極限 df dx(x)≡ limh→0 f (x + h)− f(x) h で与えられるのであった．そこで， ∂f ∂x ∼ f (x + ∆x)− f(x) ∆x と近似して，もとの微分方程式を書き換えようというのが，差分近似式によって，微分方程式を差分方 程式に転換しようというのが，差分法の第一段階である． これだけかくと，例えば， ∂f ∂x ∼ f (x)− f(x − ∆x) ∆x も差分近似なのではないかと思われるであろう．もちろんその通りである．1 階の（常）微分ですら，す でに 2 つの差分近似式がある．実は，代表的な差分近似式には，もうひとつ ∂f ∂x ∼ f (x + ∆x)− f(x − ∆x) 2∆x

というものもある．さらにもっといろいろあるのだが，ひとまず，以上のことを定義としてまとめておこう． 定義 2.1. 常微分項 df /dx あるいは偏微分項 ∂f /∂x を， f (x + ∆x)− f(x) ∆x , f (x)− f(x − ∆x) ∆x , f (x + ∆x)− f(x − ∆x) 2∆x で近似することを，それぞれ，前進差分近似，後退差分近似, 中心差分近似 と呼ぶ． これらは一体どのような近似なのかを図式的に見てみることにしよう． x x− ∆x x + ∆x f (x) 中心差分近似 前進差分近似 後退差分近似 0 xにおける f の接線 前進・後退・中心ともに点 x における関数 f (x) の接線を線分で近似していることになるわけだが，明 らかに前進・後退差分よりも中心差分の方が，よい近似であるように見えるであろう．差分近似式とい うくらいだから，われわれは近似式がどの程度，良い近似であるかを判定せねばならない．そのために 使える道具が，Taylor 展開 f (x + ∆x) = f (x) + (∆x)df dx(x) + (∆x)2 2! d2_f dx2(x) + (∆x)3 3! d3_f dx3 + O((∆x) 4 ) である．この式を用いて差分近似式の誤差評価をしてみることにしよう． 前進差分の式にこの Taylor 展開を代入してみると， f (x + ∆x)− f(x) ∆x = df dx(x) + (∆x) 2! d2_f dx2(x) + O((∆x) 2₎ となる．ここで ∆x の項が存在するので，これを前進差分近似は ∆x に関して近似度1 である，という． 後退差分と中心差分についてもやってみると， f (x)− f(x − ∆x) ∆x = df dx(x)− ∆x 2 d2_f dt2(x) + O((∆x) 2₎ f (x + ∆x)− f(x − ∆x) 2∆x = df dt(x) + (∆x)2 6 d2f dx2(x) + O((∆x) 3 ) であることから，後退差分は ∆x に関して近似度 1 であり，中心差分は ∆x に関して近似度 2 であること がわかる．このことは，前進差分や後退差分よりも中心差分の方が近似の度合いがよいことを示してい る．しかし，近似度の良い差分近似式をつくることは，いろいろな方法によってできる．ここでは，さ らに二つの例を紹介しよう． 3f (x)− 4f(x − ∆x) + f(x − 2∆x) 2∆x

なる差分近似式を考えると， 3f (x)− 4f(x − ∆x) + f(x − 2∆x) 2∆x = df dx(x)− (∆x)2 3 d3_f dx3(x) + O((∆x) 3 ) であるから，∆x に関して近似度 2 の差分近似式であることがわかる． また， 2f (x + ∆x) + 3f (x)− 6f(x − ∆x) + f(x − 2∆x) 6∆x なる差分近似式を考えると， 2f (x + ∆x) + 3f (x)− 6f(x − ∆x) + f(x − 2∆x) 6∆x = df dx(x) + (∆x)3 12 d4_f dx4(x) + O((∆x) 4 ) であるから，∆x に関して近似度 3 の差分近似式であることがわかる． このようにして f (x± s∆x) (s = 0, 1, 2, · · · ) というように，使う項を多くすれば，近似度の高い差分 近似式をつくることは，多少の計算を行うことによってできる．しかし，ここで，次のことには注意し なければならない． Remark． 微分方程式の差分解法においては，ただ単に近似度の高い差分近似式を用いれば，良 い近似解が得られるという保証はないし，同じ近似度を持つ差分近似式を用いても，近似解の振る 舞いが大きく異なることがある． 上の Remark の具体的な実例は，具体例を取り扱う各章でみることにするので，ここではひとまず注 意をめぐらしておいて欲しい． ここまで，1 階導関数の差分近似式を作ってきた．2 階以上の導関数の差分近似式のつくり方も述べて 置く必要があろう． 通常，2 階導関数の差分近似式としてよく用いられるのは， d2f dx2 ∼ f (x + ∆x)− 2f(x) + f(x − ∆x) (∆x)2 というものであろう．例によって Taylor 展開の式を代入すると， f (x + ∆x)− 2f(x) + f(x − ∆x) (∆x)2 = d2f dx2(x) + (∆x)2 12 d4f dx4 + O((∆t) 4₎ となるので，∆x に関して近似度 2 であることがわかる． 先に述べたように，近似度の高い差分近似式をつくることはできる．例えば， −f(x + 2∆x) + 16f(x + ∆x) − 30f(x) + 16f(x − ∆x) − f(x − 2δx) (∆x)2 は，∆x に関して，近似度 4 の差分近似式になっている．（検証は読者に委ねよう．） さて，3 階以上の導関数を含む微分方程式は，もちろん存在するのだが，3 階以上の導関数の差分近似 式はもうかなり形が複雑になることもあり，一般にどのようなものが良いかを論じることは難しい．こ こではそのつくり方を述べることはしないが，KdV 方程式のところでひとつの例をみることになるであ ろう． ひとまず，差分近似式をつくれるようになったところで，1 次元熱伝導方程式を例にとって，具体的に 差分法のステップに従ってやってみることにしよう．

2.3

とりあえずやってみよう ─

1

次元熱伝導方程式を例に─

本節では，1 次元熱伝導方程式を例に取り上げて，差分解法を実際に試みてみることにしよう．なぜ 1次元熱伝導方程式かを少し説明すると，2 次元以上の場合には，考えている領域の形などを問題にしな ければならず，扱いが多少面倒であるためであり，熱伝導方程式・波動方程式・Laplace 方程式という 3 つを比べてみると，やはり熱伝導方程式が一番，初学者にとっての例には適切であろうとの判断である． 波動方程式は，意外と差分法を使うのにも，テクニックが必要であり，Laplace 方程式は定常問題である から，多少初学者にとっての例としては適切でないであろうと思われる． 1次元熱伝導方程式とは，以下のように定式化される．ここでは，初期・境界値問題を考える． ∂u ∂t(x, t) = ∂2_u ∂x2(x, t) x∈ [0, 1], t ∈ [0, T ] （熱伝導方程式） u(0, t) = u(1, t) = 0 t∈ [0, T ] （境界条件） u(x, 0) = sin x x∈ [0, 1] （初期条件） ここでは，区間 [0, 1] に置かれた針金を考え，その上で初期の温度分布 sin x が与えられているとき，そ の後の時刻での温度分布を考える問題に相当している．その際，境界条件として，針金の両端での温度 が 0 とされているわけである5． さて，上で定式化された熱伝導方程式の初期・境界値問題を，先に述べていた差分解法の Step に従っ て解いてみることにしよう． まず，差分近似式を用いて離散化し，差分方程式を得る のであった．すでにいくつか見ていた差分近 似式を用いて，熱伝導方程式を離散化してみよう．左辺∂u ∂t(x, t)は，1 階偏微分の差分化として， u(x, t + ∆t)− u(x, t) ∆t なる前進差分を採用することにしよう．右辺の 2 階偏微分については，

u(x + ∆x, t)− 2u(x, t) + u(x − ∆x, t)

(∆x)2 なる差分近似式を採用することにしよう．すなわち，熱伝導方程式∂u ∂t(x, t) = ∂2_u ∂x2(x, t)を u(x, t + ∆t)− u(x, t) ∆t =

u(x + ∆x, t)− 2u(x, t) + u(x − ∆x, t)

(∆x)2 なる差分方程式として離散化したことになる． 次に，与えられた諸条件を用いて差分スキームを構成する のであった．ここでは，1 次元区間に対応す る時間柱を考えねばならない．この時間柱をメッシュで覆うためには，次のようにして考える．すなわち， 区間 [0, 1] を，N = 1/∆x 個の区間に分割する．このとき，格子点 (j, k) を時間柱内の点 (j∆x, k∆t) (0≤ j ≤ N, k = 0, 1, 2, · · · ) として定義することにすれば，時間柱が，辺の長さが ∆x, ∆t の長方形メッシュに よって覆われたことになる． 5 ここで 0 というのは，何を基準に取るかで変わってしまうけれど，基本的には，針金の両端を，例えば氷で常に冷却する などして，一定温度に保っておくことを意味している．

格子点 (i, j) における解 u の値 u(j∆x, k∆t) を uk j とかくことにする．これを元に，すでに得ていた差 分近似式を書き変えよう．すると， uk+1_j − uk j ∆t = uk j+1− 2ukj + ukj−1 (∆x)2 となる．これを整理すると，

uk+1_j = λuk_j+1+ (1− 2λ)uk_j + λuk_j−1 となる．ここで， λ = ∆t (∆x)2 である．この λ に関する条件が後に安定性などについて重要な役割を果たすことになるので注意を止め ておいて欲しい． 上で得られた差分方程式は，次図に示されるようなメカニズムに相当している． uk j−1 ukj ukj+1 uk+1_j 時刻 (k + 1)∆t 時刻 k∆t 前の時刻の 3 点の値の重み付き平均 それでは，数値解全体はどのようにして決定されるのであろうか．これを時間柱全体を見ることによっ て調べてみよう．

・・・ ・・・ ・・・ ∆t 0 2∆t 3∆t 4∆t 5∆t 6∆t 7∆t 0 ∆x 2∆x 3∆x _{1 = N ∆x} 初期条件から決定される格子点たち 境界条件から決定される格子点たち ・・・ ・・・ ・・・ ・・・ まず，初期条件を uk j を用いて表すと， u0_j = sin (j∆x) (j = 0, 1, 2,· · · , N − 1, N) となる．これによって時間柱の最下辺の格子点たちの値が決定される． 次に境界条件を uk j を用いて表すと， uk₀ = uk_N = 0 (k = 0, 1, 2,· · · ) となる．これによって時間柱の左右の辺上の格子点たちの値が決定されるわけである． ここまでで，与えられた諸条件を組み込んだ差分スキームが決定されたことになる．すなわち，すで に見た 3 点から次の時刻の 1 点での値が決まるメカニズムによって，図では○で描かれた各格子点の値 が逐次決まっていくことになる．すなわちこの差分スキームは陽的差分スキーム であることがわかる． 例えば，ここで，境界条件がなければ，差分スキームは成立しないことに注意しよう．なぜなら，最 下辺の格子点が初期条件によって決定されると，1 段目の格子点のうちで，j = 0, j = N の各点を除いた 各点での値がだけである．両端での値は決まらない．さらに 2 段目には，両端の 2 点を除いた格子点で の値が決まるだけである．このようにして，時間柱の全格子点での値を決定するためには，境界条件が 必要不可欠なのである． さて，実際に，差分解法によって数値計算を行うと，どのようなことが起こるのであろうか．残念なが ら，筆者のスキル不足により，C++で実行した数値計算例を TeX に取り込むことが出来ないため，Visual に結果をお見せすることができない． ここでは，次図のようなイメージ図でひとまず我慢してもらおう．

λ ≤ 1/2 の場合 λ > 1/2の場合 このイメージの意味するところは，λ = (∆t)/(δx)2_{の値が，1/2 を境にして，数値解の振る舞いが大} きく変化してしまうことである．λ≤ 1/2 では，初期関数の sin 曲線が，徐々に 0 に向かって落ちていく という数値解が得られ，λ < 1/2 では，数回ステップを繰り返すと，数値解は徐々に振動を始め，最終 的には，画面全体が縦線でびっしりと覆われてしまうほど激しく振動してしまう．熱伝導という物理現 象として考えれば，両端が一定温度に保たれている場合には，やがて，針金自身もその温度に落ち着く であろうことが予測される．このようなことを考えれば，λ > 1/2 の場合は，何か変なことが起こって いると判断できるであろう．こうのような場合のことを，数値的安定性がない解という．これは，各ス テップを繰り返すごとに累積していく誤差がやがて無視出来なくなって，数値解の値が発散してしまう ことから来る典型的現象である． 数値解析で注意すべきことの第一のポイントが，この数値的安定性 である．λ の選び方が重要だとい うことは，すなわち ∆t と δx の値をどの程度として選ぶかが重要だということに他ならない． 一方，あとで見るように，数値的に安定な解であるにも関わらず，実際には厳密解の近似と言えない ものが出現する場合がある．これを数値解の収束性 といい，第 2 の注意点である．これについては，状 況がやや複雑であるので，偏微分方程式の数値解析のところで，具体例をもとに触れることにする． ひとまず，数値解析には数値解析特有の難点があるだろうということはおわかり頂けたと思う．あと で具体的な常微分方程式や偏微分方程式を差分法によって解く際に，さらに詳しくいろいろな例を見る ことになるが，それらは置いておいて，次節では，数値的安定性と収束性について，やや理論的な見地 から見て置くことにしよう．

2.4

差分解の安定性と収束性

本節では，数値的安定性と収束性の問題に関して，（必ずしも数学的とは言いがたいが，）簡便な判定法 として Fourier 分解の方法を与え，さらに Lax の同等定理と Kreis の行列定理などについて述べるが，詳 しい理論は付録に譲る．実際，ここで論じるようなことを厳密にやろうとすると，例えば Lax の同等定 理の証明には，関数解析学の基本定理である「一様有界性の原理」が必要となる．なぜ関数解析的な知 識が必要となるのかについてだけ，冒頭で若干注意しておこう． 一般に微分方程式の解を考える時には，例えば 2 階の偏微分方程式であれば，未知関数は少なくとも 2階微分可能でなければならない，と考えるであろう．そうでなければ，2 階微分項が計算出来ない．こ のような立場から見た解のことを古典解 と呼ぶ．しかし，実際には，このような解のクラスだけを考え ていたのでは，物理的にも満足な結果を得ることができない．例えば，いくつかの点で微分可能性が崩 れていたり，場合によっては連続性すら崩れていても，物理的に意味のある現象を与えていることはあ る．例えば熱伝導方程式の場合，初期の温度分布が台形（方形）をしていてもなんら問題はない．この ような場合，台形の端の点では微分可能性が崩れる．また波動方程式においては，初期関数として方形

波を考えることがある．これは物理的にはなんら不自然ではない．さらに，あとで Burgers 方程式を扱 うときに見る衝撃波解は，1 点で不連続な解である．こうした通常の古典解の範囲では捉えきれない解 を扱うために，微分などの概念を（多くは積分を用いて）多少拡張して，しかるべき関数空間を設定し， その中で解の存在と一意性や解の挙動を考えることが必要になるのである．関数空間は一般には無限次 元であるから，無限次元空間の解析学が必要であり，それが関数解析学と結びつくわけである．考えて いる関数空間には，何らかの意味でノルムが入り，Banach 空間と呼ばれるクラスに属していることが多 い．「一様有界性」は，次節で述べるように，数値的安定性の Key Word であるが，Banach 空間であるよ うな関数空間に属する関数たちの族で，しかるべき条件を満たすものは，関数空間のノルムに関して一 様有界であることを保証するのが，「一様有界性の原理」なのである．ここで重要なのは，考察する関数 空間を指定しなければ，微分方程式の解に関する考察ができないことである．どのような関数空間で考 えるのかによって，解の存在や一意性，性質などが変わってくる．関数空間をどのように設定するのか は，関数解析的微分方程式論の重要なポイントである． とは言え，ここで関数空間についてごまかさない記述をしても，冗長になるので，そのあたりのこと は付録にまわそうというわけである． 2.4.1 数値的安定性の定義 今まで繰り返し数値的安定性と言ってきたが，まだその厳密な定義はしていなかった．ここでは，ま ずそれについて，考えてみることにしよう．数値的安定性と言っても，強安定や弱安定といった概念で 細かく定義することもできる6_{．しかし，ここではごく素朴に次のように定義しておこう．} 定義 2.2. ある関数空間の中の関数族として差分解{uk j} を考える時，maxj(ukj) = uk∆とかくこと にすると，ある定数 C が存在して，任意の k に対し， kuk ∆k ≤ C が成り立つとき，差分解は安定 であるという．（ここでk · k は，考えている関数空間のノルムであ り，この定義は，uk∆がこのノルムに関して一様に有界 であるとき，差分解は安定であると定義す ると主張している．） _こ こで，一様有界の意味に注意しよう．各 uk∆は有限，すなわち有界でも，一様に有界ではない．熱伝導方 程式の場合，λ > 1/2 とすると，各ステップで差分解は有界でも，だんだん振れ幅が大きくなり，発散し てしまう．どんな k に対しても，すなわちどんなに時間がたっても，差分解がある定数 C で押さえられ ていることが，一様有界の意味である7_． また念のため付け加えておくが，前節で熱伝導方程式の場合には，∆t/(∆x)2_{と 1/2 の大小で安定/不} 安定が決まったが，いつでもそのような条件であるという保証はない．それぞれの差分スキームごとに 吟味しなければならない． 6_{付録および [矢嶋・野木],[山口・野木] などを参照．} 7 ここで定数でおさえられるといっても，関数の大きさを定義するためには，ノルムを使わなければならず，ノルムを考え るためには，考えている関数空間をある程度指定しなければ意味がない．このようなことから関数解析的な議論が必要なの である．

2.4.2 Fourier分解による安定性の判定 本小節では，Fourier分解 による安定性の判定法を紹介する．この判定法の原理を簡単に言ってしま うと，任意の関数が波数の異なる三角関数の重ねあわせでかけることに注目し，どんな波数の三角関数 も時間発展で増幅されないことを，増幅率を計算することで示してやる方法である．そのために，差分 スキームが uk+1 = Suk なる形で書かれているとしよう．ここで S は差分演算子 であって， S = ∑ m cmTm という形をしている．T は平行移動作用素 とでも言うべきもので，任意の関数 f (x) に対し， T (f (x)) = f (x + ∆x) なる作用をする演算子である．ここで m∈ Z であり，cmは ∆t, ∆x を含む定数である． さて，ここで， uk_j = gnexp(iξj∆x) i =√−1 とおく．これは波数 ξ の正弦波である．g は一般に ∆x, ∆t を含んで良い． まず，平行移動演算子 T の exp(iξj∆x) への作用を調べる． T exp(iξj∆x) = T exp(iξx) (j∆x = xとおいた．) = exp(iξ(x + ∆x)) (T の定義) = exp(iξ∆x) exp(iξj∆x) であるから，

Tmexp(iξj∆x) = exp(iξm∆x) exp(iξj∆x) となる．したがって，uk+1_{= Su}k_{に S の式を代入すると，} gk+1exp(iξj∆x) =∑ m cmgmexp(iξm∆x) exp(iξj∆x) より， g =∑ m cmexp(iξm∆x) となる．この g は演算子 S によって時間を 1 ステップ ∆t だけ進めたとき，正弦波がどれだけ拡大あるい は縮小されるかを示しており，S の増幅率 と呼ばれる． いま，差分間隔 ∆x と ∆t について， ∆x = h(∆t) なる関係を保ちながら k を大きくしたとき，ある波数 ξ の増幅率が 1 よりも大きければ，その波数成分 は増幅されて，やがて発散してしまうであろう．したがって，次のようなことが言える．

¶ ³ 0≤ k∆t ≤ T をみたすどのような n, ∆t に対しても |g|n _{≤ C} を満たすような定数 C が存在すれば，差分スキームは安定であろう． µ ´ これを少し変形する． C = exp(KT )≥ (1 + KT/n)n≤ (1 + K∆t)n と見ることによって，von Neumannの安定性条件 が得られる． von Neumannの安定性条件 ¶ ³ 差分スキーム uk+1 _{= Su}k_{の差分作用素 S の増幅率 g に対して，} |g| ≤ 1 + K(∆t) をみたすような定数 K が存在すれば，差分スキームは安定である． µ ´ gは元来 ∆x, ∆t の関数だから，|g| ≤ 1 + K∆t で押さえられるために必要な関係式 ∆x = h(∆t) が決ま る． これらのことには数学的証明が必要であるが，ここでは，少なくとも気分は理解してもらえることと 思う．詳細は付録に委ねる． 2つ注意をしておこう． 陰的差分スキームは，S1uk+1 = S0ukという形をしている．ここで， S0 = ∑ p C_p(0)Tp, S1 = ∑ q C_q(1)Tq の形をしている．このとき，同様に，uk j = gkexp(iξj∆x)を代入して，差分作用素 S = S1−1S0の増幅率 を求めると， g = ∑ pC (0) p exp(iξp∆x) ∑ qC (1) q exp(iξq∆x) が得られる． 第 2 に，もし連立偏微分方程式系を差分化して，連立差分方程式系を得た場合には，差分演算子は行 列となり，したがって，cmや増幅率 g も行列となるが，これを増幅行列 という．これが k 乗しても発散 しないためには，増幅行列 g のスペクトル半径(固有値の絶対値の最大値)γ が一様に有界であればよい． このことから，多次元版の von Neumann の条件として， γ = 1 + K∆t なる定数 K が存在することが必要である． 実は，1 次元の場合は，von Neumann の安定性条件は必要かつ十分であり，2 次元以上では一般に十 分ではないことが知られている． Fourier分解の具体的な適用は，偏微分方程式の数値解析のところで，繰り返し見るであろうから，こ