ここでは,3元連立常微分方程式系として与えられる,カオス的な振る舞いをする2つの力学系を紹介 しよう.その振る舞いが数値解析的に捉えられるかどうかを見たいのである.
レスラー系とは,次で与えられる3元連立の常微分方程式系である.
dx
dt =−y−z dy
dt =x+αy dz
dt =b−cz+xz
この常微分方程式系を,すべての微分項を前進差分によって差分化することによって,数値的に解く ことを考えよう.差分方程式は,以下で与えられる.
xn+1 =xn−(∆t)(yn+zn) yn+1 =yn+ (∆t)(xn+ayn) zn+1 =zn+ (∆t)(b−czn+xnzn)
例えば,a= 0.2, b = 0.2, c= 5.7というパラメータの値を設定すると,面白い形が浮かび上がってきま す.注目すべきなのは,レスラー方程式においては,非線形項が最後の方程式のxzという項しかないこ とである.いかにカオス的という現象が微妙な現象であるかを物語っていると言えるだろう.ちなみに,
次のような方程式系だと,形は似たようになりますが,やがて軌道は周期的になってしまい,カオス的 な振る舞いをしません.これも前進差分で差分化することによって,数値的に解くことで捉えることが できます.最後の方程式でbがbxになっているだけです.
dx
dt =−y−z dy
dt =x+αy dz
dt =bx−cz+xz やはり前進差分で差分化することによって得られる方程式は,
xn+1 =xn−(∆t)(yn+zn) yn+1 =yn+ (∆t)(xn+ayn)
zn+1 =zn+ (∆t)(bxn−czn+xnzn) で与えられる.
この例9もカオス的という現象の微妙さを表していると言えるだろう.
9実は,[山口2]にあった例だが,山口昌哉氏のミスプリントであった.弘法も筆の誤りである.
ローレンツ系とは,次で与えられる3元連立の常微分方程式系であり,気象学者ローレンツが,カオ ス的な振る舞いをする方程式として導いたものである10.
dx
dt =σ(x−y) dy
dt =γx−y−xz dz
dt =−bz+xy やはりすべての微分項を前進差分によって差分化すると,
xn+1 =xn+ (∆t)σ(xn−yn)
yn+1 =yn+ (∆t)(γxn−yn−xnzn) zn+1 =zn−(∆t)(bzn−xnyn)
となる.パラメータをσ = 10, γ = 26, b = 2.667とすると,二つの目を持つような軌道が現れる.この方 程式の解の図は,意外と多くの本に載っている.
このようなカオス的な振る舞いをする微分方程式系についても,差分法は比較的有用な手法を提供し ていると言えるであろう.
3.8 1 次元楕円型境界値問題
常微分方程式の差分法による数値解析を見る上で,最後に,1次元楕円型境界値問題を扱うことにし よう.
1次元楕円型境界値問題とは,以下で与えられる常微分方程式を,しかるべき境界条件のもとで解く という問題である.
−d2u
dx2(x) +q(x)u(x) =f(x) (0< x <1) ここで設定する境界条件とは,以下の3つの場合を考えることにする.
u(0) = 0, u(1) =a (Derichlet境界条件)
du
dx(0) =b, du
dx(1) =c (N eumann境界条件) du
dx(0) =b, u(1) =a (混合境界条件)
1次元楕円型境界値問題は,今まで扱ってきた常微分方程式とはかなり趣の違うものである.それは,
今まで考えてきた常微分方程式が,ある量の時間発展を記述するという意味で,非定常問題であったの に対し,1次元楕円型境界値問題は,1次元区間上の定常分布を決定するという定常問題だからである.
従って,差分法によって数値的に解を求めるとは言っても,今までのもののように,∆tを含んだ漸化式
10導出はここでは省略する.例えば[丹羽]を参照されたい.
を用いて逐次値を決めていくなどということをするわけではない.むしろ,考えている区間を細分し,そ の有限個の点での値を決定すれば良い.具体的に差分法を試みればおのずとわかることであるが,この 種の問題の場合には,連立一次方程式を1度解けば,各分点上での関数の値が一斉に求まる.しかし,分 点の個数だけ方程式があるので,例えば100元連立一次方程式を解いたりしなければならず,ここでは,
大規模連立一次方程式を数値計算する必要性に迫られる.本節では,連立一次方程式の解法を説明する ことはせず,付録にまわすことにして,解くべき連立一次方程式を得ることを主眼として差分法を試み よう.
さて,例によって,常微分方程式を差分化しよう.ここでは,空間差分すなわちxに関して差分近似 を行うのであるが,差分近似式としては最も代表的な
d2u
dt2 ∼ u(x+ ∆x)−2u(x) +u(x−∆x) (∆x)2
を採用しよう.また,以下では,特に簡単な場合として,q(x) = λ, f(x) = −µxという場合を考えるこ とにしよう.
d2u
dx2(x) =λu(x) +µx (0< x <1) を考えるわけである.これを差分化すると,
u(x+ ∆x)−2u(x) +u(x−∆x)
(∆x)2 =λu(x) +µx
となる.これを漸化式の形にすると,
um+1−2um+um−1
(∆x)2 =λum+µm∆x となる.変形すると,
um+1−(2 +λ(∆x)2)um+um−1 =µm(∆x)3
となる.次に境界条件を差分化する.ここでNeumann境界条件は前進差分で差分化する.
u0 = 0, uN =a (Derichlet境界条件)
u1−u0
∆x =b, uN −uN−1
∆x =c (N eumann境界条件) u1 −u0
∆x =b, u(1) =a (混合境界条件)
となる.
さて,実際にどのようなメカニズムで差分解が決定されるのかを少し丁寧に各場合ごとに見てみるこ とにしよう.
のDerichlet境界条件の場合.このとき,分点{j∆x} (j = 0,1,· · · , N)での値は,次のような連立 一次方程式を解くことによって得られる.
−(2 +λ(∆x)2) 1
1 −(2 +λ(∆x)2) 1 . ..
. .. 1
1 −(2 +λ(∆x)2)
·
u1 u2 ... ... uN−1
=
µ(∆x)3 2µ(∆x)3
... ...
(N −1)µ(∆x)3−a
これを図式的に表すと以下のようになる.
0 · · · 1 連立一次方程式によって一斉に値が決定される.
Derichlet境界条件によって決定される.
つまり,両端のu0, uNの2点はDerichlet境界条件によって値が決定されており,u1,· · ·uN−1では,3点 の間の関係式が与えられており,各点での値は,連立一次方程式を解くことによって一斉に決定される というわけである.(このメカニズムは,今まで見てきた常微分方程式の差分解析によるメカニズムとは 異なっていることに注意.)
次に のNeumann境界条件の場合.この場合も,以下で与えられる連立一次方程式によって各分点で
の値が決定される.
−1 1
1 −(2 +λ(∆x)2) 1
1 −(2 +λ(∆x)2) 1 . ..
. .. 1
1 −(2 +λ(∆x)2) 1
−1 1
·
u1
u2
... ... ... uN−1uN
=
b∆xµ(∆x)3 2µ(∆x)3
... ...
(N−1)µ(∆x)3 c∆x
やはりこれを模式図で表してみよう.
0 · · · 1
連立一次方程式によって一斉に値が決定される.
Neumann境界条件による関係式が与えられている.
ここでは,両端の2点に対してはNeumann境界条件から決まる関係式が与えられており,u1,· · · , uN−1 には3点の間の関係式が与えられている.これらをまとめて連立一次方程式を解くことによって一斉に 決定することになる.
最後に の混合境界条件の場合.この場合は, と の変形版である.連立一次方程式は,以下で与
えられる.
−1 1
1 −(2 +λ(∆x)2) 1
1 −(2 +λ(∆x)2) 1 . ..
. .. 1
1 −(2 +λ(∆x)2)
·
u1
u2 ... ... ... uN−1
=
b∆x µ(∆x)3 2µ(∆x)3
... ...
(N −1)µ(∆x)3−a
これをやはり図式的に表してみる.
0 · · · 1
連立一次方程式によって一斉に値が決定される.
Derichlet境界条件によって決定される.
Neumann境界条件による関係式が与えられている.
これは と の混合版であるから再度説明する必要はないであろう.
さて,1次元楕円型境界値問題を考える場合には,次のような著しい事実が知られている.
事実:1次元楕円型の差分スキームは,無条件安定であり,∆xを小さくすれば,差分解は厳密解に 収束する.
この意味で1次元楕円型境界値問題はそれほど心配することなく差分法を試みることが出来る.問題 点は,すでに述べたように,分点の数が多ければ多いほど,解くべき連立一次方程式が大きくなるので,
大規模連立一次方程式を精度良くしかも効率的に解く必要があるという点である.
また,2次元以上の問題になると,これはまた様相が変わってくる.その理由は,考えている領域の形 が問題になるからである.これについては,後の偏微分方程式の数値解析のところでも多少触れること にするが,有限要素法の内容にむしろ近い内容であるので,差分法の節ではごく簡単に触れるにとどめ る.
4 偏微分方程式の差分法による数値解析
概 要
本章では,(ほぼすべて1次元の問題ではあるが)いくつかの偏微分方程式を取り上げて差分法に よって数値的に解くことを試みる.扱う偏微分方程式は,線形・非線形から様々なものを選ぶが,常 微分方程式以上に微妙な問題が多いことを実感してもらうと同時に,差分法による数値解析が威力を 発揮する面もあることも実感して頂きたい.前者の微妙な問題としては,常微分方程式のところです でに見たように,安定性と収束性の問題であると言える.理論のところで見たようにLaxの同等定理 などが収束性と安定性の同値性を保証してくれる局面もあるが,そうでない場合もある.ここでは非 線形拡散方程式のひとつを例に,数値的安定性を持っていても収束性を持たない例も扱う.また,波 動方程式などでは,厳密解では本来波束が拡散しないはずなのに,拡散効果が現れてしまう例も扱う.
一方,後者の例としては,爆発解,衝撃波解,2波相互作用,ソリトン解などが挙げられる.微分方 程式の数値解析は,こうした解をシミュレーションすることで,微分方程式の理論(厳密解の挙動の 研究)にも重要な役割を果たしうるのである.
4.1 偏微分方程式についての一般的注意
偏微分方程式を差分解法によって数値的に解くメカニズムは,すでに熱伝導方程式を例に取ってみて いた.その意味で基本的なところは押さえてあるわけで,すぐに様々な方程式の数値解析に突入しても 良いのだが,ここでは,偏微分方程式について,少し一般的な注意を極めて簡単にだが述べておくこと にしよう.
偏微分方程式は,未知関数の偏導関数たちの間の関係式のことであった.やはり線形と非線形という 分類は可能である.中でも2階線形偏微分方程式には典型的な3つの型,すなわち,放物型・双曲型・楕 円型が知られている.重要なことは,偏微分方程式の解の性質が,型ごとに特徴的であるということで ある.まず,非定常現象を扱う方程式という意味で放物型・双曲型,定常現象を扱うという意味で楕円 型とに大別できる.放物型と双曲型の違いとしては,例えば,解の不可逆性の点にある.熱伝導現象と いうのは不可逆現象であるが,波動は,例えば摩擦による減衰などを考えずに理想化して考えると,可 逆現象である.また放物型では一般に最大値原理とよばれる定理が成り立つが,双曲型では一般に成り 立たない.また楕円型では最大値原理が成り立ち,調和関数との関係が深い.偏微分方程式の解の性質 や特徴的な解を差分法で捉えることが出来るかどうかは,差分法による数値解析入門講義としては,極 めて重要な視点であるから,そうしたことについては,ある程度は触れざるを得ない.そこで,こうし た偏微分方程式の解析的理論については,本文および付録で若干初等的な部分を解説するが,もちろん 本格的な議論は,他書に譲らざるを得ない11.
ここで次に述べておくべきことは,偏微分方程式の代表的な解法として知られている求積法と変数分離法 である.もちろん偏微分方程式の解法はこれだけではなく,実に多種多様なものがあるわけだが,ここ ではひとまずこの2つの解法が重要なので,特に取り上げる.
求積法とは,(偏)微分方程式を積分して,一般解を直接求める方法のことをさす.しかし,すでに何 度か繰り返して来たように,偏微分方程式が求積法によって解ける場合は極めて稀であって,その可能 性はほぼ絶望的だと思ったほうが良い.求積法で解ける代表的な偏微分方程式としては,2階波動方程式 およびその「縮小版」とでも言うべき1階波動方程式が挙げられる.この詳細については,波動方程式 の箇所と付録で特性曲線法との関わりを中心に説明する.
変数分離法とは,一般解を特殊解の重ね合わせとして求める方法で,例えば,未知関数をu(t, x) =
T(t)X(x)などとおいて常微分方程式を解くことに帰着し,特殊解を求めた上で,それらを重ね合わせる
(線形重畳)という方法が代表的である.この方法は,比較的多くの偏微分方程式を解く際に用いられて おり,汎用性は少なくとも求積法よりは良いと言えるが,最終的に線形重畳することによって,解を無
11例えば,[金子],[井川],[MKI]あたりを参照されたし.