4.2 熱伝導方程式
4.3.4 波束の拡散
さて，安定な差分スキームとして， のタイプのものが得られた．これでめでたしめでたしである．あ とは数値実験を試みよう．そう上手くはいかないのが数値解析の難しさである．ここまでの段階で考察 したことを元に，数値実験を行うと，もちろん ・ の差分スキームでは不安定となるが， では，安 定条件を満たす∆t,∆xを取っても，ほ´と´ん´ど´の´場´合，波束の拡散´ と呼ばれる現象が起きる．
波束の台が広がり振幅が低くなっていく．
この現象は上図のように，波束の台（正値の領域）が徐々に広がっていき，波束の振幅も徐々に低くなっ ていく現象である．いま，1階波動方程式の一般解は任意関数fを用いて，
u(x, t) = f(x−ct)
とかけていた．これは，初期関数fが形を変えずに，右向きに進行して行く解である．従って，このよ うな波束の崩れは起きないはずである．これは，物理現象と同時に，数学の解析的理論と矛盾した数値 解なのである．
このようなことはどうして起こるのであろうか．これも前節の時間柱を利用した考察によって理解で きる．確かに のスキームでは，依存領域の初期関数のデータは取り込んでいた．しかし，もし∆t/∆x が，1/cよりも真に小さいと，特性曲線よりも左側，すなわち，本来依存領域に含まれないような領域の 初期関数のデータまで取り込んでしまうことになる．このことが，波束の拡散という現象を引き起こし たのだと考えられる．
本来依存領域に属さない初期データ 特性曲線 傾きλ = ∆t/∆xの直線
このことを考えれば，ちょうどλ= 1/2のときには，波束の拡散のない解が得られるはずである．実 際数値計算をしてみれば，そのようになることがわかる．しかし，このような比較的簡単なスキームで，
比較的単純な偏微分方程式を解く場合にはともかくも，複雑な偏微分方程式系を解く場合には，差分ス キームの安定性条件の限界ぎりぎりに設定することは，数値的誤差などの点であまり良くない．
このような波束の拡散を抑制する差分スキームを考えるには，少し複雑な表式を考えねばならない．以 下のような3つのものが知られている．
uk+1j = 1
2(ukj+1+ukj−1)− c 2
∆t
∆x(ukj+1−ukj−1) uk+1j = ukj−1−c∆t
∆x(ukj+1−ukj−1) uk+1j = ukj − 1
2 (∆t
∆x )
(ukj+1−ukj−1) + c2 2
(∆t
∆x )2
(ukj+1−2ukj +ukj−1)
ここで挙げたものは，どれもλ≤1/cで安定である．2番目のものは時間に関して中心差分を取ってい るので，第1ステップだけは別の方法で計算せねばならないことに注意しておく．
4.4 2 階波動方程式
2階波動方程式の初期・境界値問題を考えよう．2階波動方程式は，以下で与えられる．
∂2u
∂t2 =4u
この2階波動方程式を，簡単な方法で差分化することを考えよう．すなわち，時間・空間について，2
階微分項を一番典型的な差分近似するのである．
u(t+ ∆t, x)−2u(t, x) +u(t−∆t, x)
(∆t)2 =cu(t, x+ ∆x)−2u(t, x) +u(t, x−∆x) (∆x)2
として差分方程式に直すと，
uk+1j =λ2cukj+1+ 2(
1−λ2c)
ukj +λ2cukj−1−ukj−1
となる．2階波動方程式は，時間に関して2階の微分項を持つために，差分スキームは，例えば熱伝導方 程式などと比べて複雑なものになる．メカニズムを見ると，次図のようになっている．
時刻(k+ 1)∆t
時刻k∆t 時刻(k−1)∆t
(j+ 1)∆x (j−1)∆xj∆x
つまり，時刻について，ひとつ前だけではなく二つ前の値が必要になるわけである．従って，格子点 での値を逐次決めていくためには，もう少し吟味が必要である．それには，波動方程式の初期・境界値 問題について少し詳しく見ておく必要がある．
時間についての2階微分項が存在するために，初期関数だけを与えても偏微分方程式は解けない．一 般に，2階波動方程式の初期・境界値問題は，以下のように定式化される．
∂2u
∂t2 =4u t∈[0,∞), x∈[0,1]
u(x,0) = φ(x), ∂u
∂t(x,0) =ψ(x) u(0, t) = 0, u(1, t) = 0
ここでは両端を固定して考える固定端の問題に限定する．（Neumann境界条件などを課せば，自由端の 問題も論じることはできる．）
初期関数についての条件，および境界条件を差分化して，
uk0 =ukN = 0, u0j =φ(j∆x)
となるのは良い．問題は，初期時刻における未知関数の微分に関する条件をどう処理するかが問題であ る．いずれにしても，差分化が必要である．これによってu1が計算できるようにしたいのである．
u(t+ ∆t)−u(t)
∆t =ψx
と差分近似すると，
u1j =u0j + (∆t)ψ(j∆x) なる差分方程式が得られる．
これらを利用すると，差分スキームが構成される．
境界条件により定まる．
初期条件により定まる．
未知関数の微分に関する差分方程式から求まる．
注目すべきことは，未知関数の微分によって，1段目から2段目の格子点での値が決定されなければ，波 動方程式を差分化した差分方程式が適用できないということである．この初期関数の微分に関する条件 が本質的なのである．
このように考えると，初期関数の微分に関する条件を差分化した時の近似度を考えたくなる．いま，2 階の微分項を近似した差分近似式の近似度は2であったのに対し，初期関数の微分に関する条件は，近 似度1でしかない．これは，数値計算をする上ではかなり危ういことをしていると思った方が良い．こ れを解決するためには，例えば，以下のような方法を用いる．
まず，u(x, t)をtに関して，Taylor展開すると，
u(x, t) =u(x,0) + ∆t∂u
∂t(x,0) + (∆t)2 2
∂2u
∂t2(x,0) +O((∆t)3) となる．ここでt= 0でも波動方程式が成り立っているとみなして，
∂2u
∂t2(x,0) =c∂2u
∂x2(x,0) =cu0j+1−2u0j +u0j−1 (∆x)2 と近似することによって，
u1j =u0j + ∆tψ(j∆x) + c 2
(∆t
∆x )2
(u0j+1−2u0j +u0j−1)
という差分方程式を得る．これを用いて1段目から2段目を計算することにすれば，精度はあがると考 えられるわけである12．
さて，以上の考察で得られた差分スキームの安定性を吟味しておかねばならない．記号の簡単のため c= 1として計算しよう．（得られた結果を，安定性条件として1/cが効いてくるように変更すればよいだ けである．）例によって，
uk+1j =λ2cukj+1+ 2(
1−λ2c)
ukj +λ2cukj−1−ukj−1
の増幅率を計算しよう．これは少し計算が面倒ではある．ukj =gkexp(iξj∆x)を代入することによって，
gに関する2次方程式
g2−2g+ 1 = 4g (∆t
∆x )2
sin2
(ξ∆x 2
)
となる．これを解きたい．
a= ∆t
∆xsin
(ξ∆x 2
)
とおくと，
g2−2(1−a2)g+ 1 = 0 となる．|a| ≤1であれば，
g = 1−2a2±i2a√ 1−a2 となって，
|g|= 1
となる．一方，|a|>1であれば，gは実数でg1 6=g2かつg1g2 = 1であることから|g1|,|g2|のいずれかは 1より大きくなる．
|a|= (∆t
∆x )
·¯¯
¯¯sin
(ξ∆x 2
)¯¯¯¯
であるから，以上のことをまとめると，次の結論を得る．
(1) ∆t
∆x ≤1なら|a| ≤1であって，|g1|=|g2|= 1となって安定．
(2) ∆t
∆x >1なら|a|>1となるようなξが存在するので，max (|g1|,|g2|)>1となって不安定．
このようにして2階波動方程式の数値解の安定性が調べられた．実際に数値実験してみても，このよ うな安定/不安定性が見られる．その理由は，1階波動方程式の場合と同様であるが，ここで注意すべき ことは，1階の場合と違い，2階の場合には，2本の特性曲線が存在するということである．影響領域・
依存領域はいずれも2本の直線ではさまれた領域となる．
12ただ，筆者の数値実験では，どちらの手法でも，それほど違いが見られなかった．[金子]には各段に違うと書いてあるの だが・・・
依存領域 影響領域
bababababababababababababababababababababababab
以上の節が差分法を用いた偏微分方程式の数値解析の初等的で一般的な議論である．以下の節 は各論である．ここまでで行ったような詳細な検討はしないが，それぞれの各論を取り上げたの には意味がある．各節冒頭でそれについて説明してあるので注意されたい．