波動方程式

山本昌志^{∗}

2005 年 2 月 14 日

1 波動方程式とは

ラプラス方程式が済んだので、次に波動方程式に移ろう。その前に、 2 階の偏微分方程式の種類について説明しておく。 2 階の偏微分方程式は、ラプラス方程式のように楕円型、次に学習する波動方程式のような双曲型、学習はしないが拡散方程式のような放物型に分けられる。これが、 2 階の偏微分方程式の代表的な型である。これらの解法を知っておけば、自然現象の多くの問題を計算することができる。いうなれば、超基本の方程式である。

波動方程式は、名前が表しているように波の方程式である。自然科学では、波を扱うことが非常に多い。光、電磁波、量子力学等の問題は全て波を取り扱っている。いろいろな場面で出くわす波の方程式は簡単で、

\frac{∂^{2}f}{∂x^{2}} + \frac{∂^{2}f}{∂y^{2}} + \frac{∂^{2}f}{∂z^{2}} = \frac{1}{c^{2}} \frac{∂^{2}f}{∂t^{2}}  (1)

と書き表すことができる。 c は波の速度である。これは、 3 次元の場合で、時間を入れると 4 次元の方程式になり、ちょっと計算するには複雑である。そこで、ここでは空間 1 次元、時間 1 次元の 2 次元の方程式

\frac{∂^{2}f}{∂x^{2}} = \frac{1}{c^{2}} \frac{∂^{2}f}{∂t^{2}}  (2)

を数値計算で解くことを考える。

皆さんは、フーリエ級数を学習したときに、この方程式を解いたとはずである。ここでは、数値計算により近似解を得る方法を学習する。もちろん、フーリエ級数で解いた解は、解析解で完璧である。ただ、フーリエ級数が適用できるのは、空間が 1 次元の場合である。 2 次元以上になると境界条件が簡単な場合に限り、フーリエ級数を用いて計算できる。境界が複雑になると、数値計算で近似解を求めることが重要になる。数値計算は、空間が 2 次元以上の問題で威力を発揮することになるが、ここでは学習のため、空間が 1 次元の問題を解くことにする。

具体的な問題を例にして、学習を進める。比較的単純な問題として、図 1 のような弦の振動を考える。これは、ギターのように両端が固定された弦である。ある時刻 t の位置 x の変位を u(x, t) としている。この変位は波動方程式、

\frac{∂^{2}u}{∂x^{2}} = \frac{∂^{2}u}{∂t^{2}}  (3)

∗独立行政法人 秋田工業高等専門学校 電気工学科

を満たす。ただし、波の速度は c = 1 とした。こうしても、波動方程式を解くと言う意味はそうは変わらないし、計算が楽になるメリットはある。

図 1: 時刻 t の弦の様子。

2 差分法による 1 次元波動方程式の数値計算

このあたりの説明は、参考文献 [1] を大いに参考にした。これは分かりやすい教科書なので、読んでみると良いだろう。

2.1 差分方程式

1 次元波動方程式を数値計で解くことを考える。その前に、解くべき方程式と条件をきちんと書いておく。解くべき方程式と条件は、

\{ \begin{matrix}\frac{∂^{2}u}{∂x^{2}} = \frac{∂^{2}u}{∂t^{2}} & (0 ≦ x ≦ 1, 0 ≦ t) \\ u(x, 0) = ϕ(x),  \frac{∂u}{∂t}(x, 0) = ψ(x) & (0 ≦ x ≦ 1) \\ u(0, t) = u(1, t) = 0 \end{matrix}  (4)

となる。弦を伝わる波の速度は 1、弦の長さも 1 としている。この最初の式は波動方程式であるが、 2 番目を初期条件、 3 番目を境界条件と言う。 2 番目の初期条件は、 t = 0 の時の弦の状態を示しており、 ϕ(x) はそのときの弦の形 (変位)、 ϕ(x) は弦の変位の速度である。

波動方程式の他に、初期条件と境界条件がある。力学的状態は、ある時刻、ここでは t = 0 の時の変位とその変位の速度が決まれば、それ以降を決めることができる。振動の場合は、これに加えて更に、振動の境界条件を決める必要がある。これらが決まって初めて、波動方程式とともに、振動の状態、ある時刻と位置の変位の値が決まるわけである。図 4 に初期条件と境界条件の様子を示す。図 2: 時刻 t = 0 のときの弦の様子 (スナップショット)。初期条件と境界条件が表されており、 y 方向の速度が ψ(x) になっている。

まずは、波動方程式を差分方程式に書き直すことからはじめる。これも、いつものように、解 u(x, t) をテイラー展開する。 x 方向の微小変位を δx、時間軸方向の微小変位を δt とする。すると、

u(x + Δx, t) = u(x, t) + \frac{∂u}{∂x}Δx + \frac{1}{2!} \frac{∂^{2}u}{∂x^{2}}(Δx)^{2} + \frac{1}{3!} \frac{∂^{3}u}{∂x^{3}}(Δx)^{3} + \frac{1}{4!} \frac{∂^{4}u}{∂x^{4}}(Δx)^{4} + ⋯
u(x − Δx, t) = u(x, t) − \frac{∂u}{∂x}Δx + \frac{1}{2!} \frac{∂^{2}u}{∂x^{2}}(Δx)^{2} − \frac{1}{3!} \frac{∂^{3}u}{∂x^{3}}(Δx)^{3} + \frac{1}{4!} \frac{∂^{4}u}{∂x^{4}}(Δx)^{4} − ⋯

(5)

となる。これらの式の辺々を足し合わせえると、

\frac{∂^{2}u}{∂x^{2}} |_{x,y} = \frac{1}{Δx^{2}} [u(x + Δx, t) − 2u(x, t) + u(x − Δx, t)] − O(Δx^{2})  (6)

が得られる。このことから、 2 階の偏導関数の値は微小変位 Δx の場所の関数の値を用いて、 (Δx)^{2} の精度で近似計算ができることが分かる。すなわち、式 ( 6) の右辺の第 1 項を計算すればよいのである。ラプラス方程式と同じである。同様なことを時間軸方向についても行うと

\frac{∂^{2}u}{∂t^{2}} |_{x,t} = \frac{1}{Δt^{2}} [u(x, t + Δt) − 2u(x, t) + u(x, t − Δt)] − O(Δt^{2})  (7)

が得られる。これらの式 (6) と (7) を元の波動方程式 (4) に代入すれば、

\frac{1}{Δx^{2}} [u(x + Δx, t) − 2u(x, t) + u(x − Δx, t)] = \frac{1}{Δt^{2}} [u(x, t + Δt) − 2u(x, t) + u(x, t − Δt)] (8)

となる。これが、 1 次元波動方程式の差分の式である。この式を計算し易いように、もう少し変形すると、

u(x, t + Δt) = 2u(x, t) − u(x, t − Δt) + \frac{Δt^{2}}{Δx^{2}} [u(x + Δx, t) − 2u(x, t) + u(x − Δx, t)]  (9)

とすることができる。この式の右辺は、時刻 t と t − Δt の値でである。そして、左辺は時刻 t + Δt の値である。このことから、式 (9) を用いると、時刻 t と t − Δt の値から、 t + Δt の値が計算できることになる。

実際に式 (9) を数値計算する場合、 x 方向には Δx、時間軸方向には Δt 毎に分割する。ラプラス方程式を格子点で分割したのと同じである。格子点に分割し数値計算する場合、 u(x, t) や u(x + Δx, y) と表現するよりは、 u_{i j} と表現したほうが便利である。そこで、

\begin{aligned}u(x, t) & = u(iΔx, jΔt) \\ & = u_{i j} \end{aligned}  (10)

と表現を改める。このようにすると、式 (9) は

u_{i j+1} = 2u_{i j} − u_{i j−1} + α (u_{i+1 j} − 2u_{i j} + u_{i−1 j})  (11)

となり、数値計算し易い形になる。ただし、

α = \frac{Δt^{2}}{Δx^{2}}  (12)

である。この式を用いた計算の様子を図 3 に示す。

図 3: 差分方程式の計算の様子

波動方程式というけったいな偏微分方程式が、ただ単に数値を順番に代入していく式に変換されたわけである。この計算は非常に簡単である。ただ、時間領域を 1000 分割 (N_{t} = 1000)、 x 軸領域も 1000 分割 (N_{x} = 1000) すると、 100 万回の計算が必要であるが、コンピューターにとって、その程度の計算は大したことはない。

2.2 初期条件

式 (11) を計算すると、 t = 0 の状態から、時間の経過によって弦の様子がどうなるか分かる。以下のように、芋づる式に、弦の変位が計算できるわけである。

u_{1 0} u_{2 0} u_{3 0} u_{4 0} u_{5 0} ⋯ u_{N_{x}−1 0}
⇓
u_{1 1} u_{2 1} u_{3 1} u_{4 1} u_{5 1} ⋯ u_{N_{x}−1 1}
⇓
u_{1 2} u_{2 2} u_{3 2} u_{4 2} u_{5 2} ⋯ u_{N_{x}−1 2}
⇓
u_{1 3} u_{2 3} u_{3 3} u_{4 3} u_{5 3} ⋯ u_{N_{x}−1 3}
⋮
u_{1 N_{t}} u_{2 N_{t}} u_{3 N_{t}} u_{4 N_{t}} u_{5 N_{t}} ⋯ u_{N_{x}−1 N_{t}}

このように、計算を盲目的に進めれば、弦の振動の式 (4) の数値計算の結果である近似解が得られる。当然、境界条件

u_{0 j} = u_{N_{x} j} = 0  (j = 0, 1, 2, 3, ⋯, N_{t})  (13)

を、忘れてはならない。

これを計算するためには、まず、 u_{i 0} (i = 1, 2, 3, ⋯, N_{x} − 1) の値を決める必要がある。これ以前の状態が分からないので、式 (11) は使えないが、式 (4) の初期条件が使える。すなわち、

u_{i 0} = ϕ(iΔx)  (14)

である。

次に、 u_{i 1} (i = 1, 2, 3, ⋯, N_{x} − 1) を計算するわけであるが、まだ、式 (11) は使えない。なぜならば、この式は 2 つ前の状態まで必要なので、これまでのところ、一つ前の状態しか分かっていないからである。そこで、 2 番目の初期条件 (変位の速度) を使うことになる。計算したい量は u(x, Δx) なので、とりあえずテーラー展開してみる。これを、 t = 0 の周りでテーラー展開すると、

u(x, Δt) = u(x, 0) + \frac{∂u}{∂t}Δt + \frac{1}{2!} \frac{∂^{2}u}{∂t^{2}}(Δt)^{2} + O(Δx^{3})

初期条件と波動方程式より

= u(x, 0) + ψ(x)Δt + \frac{Δt^{2}}{2} \frac{∂^{2}u}{∂x^{2}}(Δt)^{2} + O(Δx^{3})

(15)

となる。この右辺の第 1 と 2 項は簡単に計算できる。問題は第 3 項であるが、これは見覚えのある式である。式 6 と同じである。これを代入すると、

u(x, Δt) ≈ u(x, 0) + ψ(x)Δt + \frac{Δt^{2}}{2Δx^{2}} [u(x + Δx, t) − 2u(x, t) + u(x − Δx, t)]

≈ u(x, 0) + ψ(x)Δt + \frac{α}{2} [u(x + Δx, t) − 2u(x, t) + u(x − Δx, t)]

(16)

となる。これは、めでたい式である。右辺は、 t = 0 のみの値で構成されている。これで、 u_{i 1} (i = 1, 2, 3, ⋯, N_{x} − 1) が計算可能になった。この式から、

u_{i 1} = u_{i 0} + ψ(x_{i})Δt + \frac{α}{2} [u_{i+1 0} − 2u_{i 0} + u_{i−1 0}]  (17)

が得られる。以上より、 u_{i 0} と u_{i 1} が得られたわけである。 u_{i 2} 以降は、式 (11) に従い、計算すればよい。

3 実際の計算方法

ここで、実際のプログラムを作成する上で、一つだけ注意を与えておく。プログラムが安定に動作するためには、位置と時刻の微小変化の比の 2 乗は、

α ≤ 1  (18)

である必要がある。この値は、繰り返し計算をすると、その回数だけ乗算される。即ち、 n 回の繰り返し計算があれば、 α^{n} が現れる。もし、 α が 1 以上であれば、早急に発散するであろう。従って、これは 1 以下になるように、微小変化を考えなくてはならない。

4 練習問題 (定在波)

ラプラス方程式のプログラムを参考にして、練習問題のプログラムを作成せよ。

図 4 のようにギターの弦を留め金でとめている。 t = 0 の瞬間に留め金をはずした場合、その振動はどうなるか?。 t = 1 まで、振動の様子を数値計算で求め、それをアニメーションで表示せよ。予め頭で想像したものと結果は大きく食い違うはずである。非常に、興味深い結果が得られるはずである。この問題は、フーリエ級数の時間に学習した??。

ヒントを与えておく。 t = 0 のとき、止められていた弦が動き始める。従って、このときの速度はゼロであるので、初期条件

ψ(x) = 0  (19)

となる。図 4: 問題の弦

7

5 進行波の取り扱い

先の練習問題、弦を三角形に張った後の様子は、定在波である。ここでは、進行波の記述方法について、コメントしておく。進行波を数値計算すると面白いのでその方法を示す。進行波を記述するためには、初期条件さえ記述すれば、後の差分方程式は同じである。その初期条件の記述の仕方を示す。

元の波動方程式

\frac{∂^{2}u}{∂x^{2}} = \frac{1}{c^{2}} \frac{∂^{2}u}{∂t^{2}}  (20)

には、明らかに、ダランベールの解

u(x, t) = f(x − ct) + g(x + ct)  (21)

というものがある。これは元の波動方程式に代入すれば、それを満足していることは直ちに理解できる。ここで、 f(x − ct) は x 軸を正の方向に進む進行波 (forward wave) で、 g(x + ct) は負の方向に進む後進波 (backward wave) である。

初期条件

u(x, 0) = ϕ(x)  (22)

の波が x 軸を正の方向に進む進行波として取り扱うには、どうしたらよいだろうか?。のこる条件は、

\frac{∂u}{∂t}(x, 0) = ψ(x)  (23)

である。進行波になるように、 ψ(x) を決めればよい。 u(x, t) を進行波と仮定すると、式 (22) から

u(x, Δt) = ϕ(x − cΔt)  (24)

となる。この式を使って、 ψ(x) を求めることにする。 ψ(x) の定義より、

ψ(x) = \frac{∂u}{∂t}(x, 0)

= lim_{Δt→0} \frac{u(x, Δt) − u(x, 0)}{Δt}

= lim_{Δt→0} \frac{ϕ(x − cΔt) − ϕ(x)}{Δt}

cΔt = Δx とおくと

= −c lim_{Δx→0} \frac{ϕ(x) − ϕ(x − Δx)}{Δx}

= −c\frac{dϕ}{dx}

(25)

となる。進行波にするためには、 ψ(x) は ϕ(x) の導関数ににすればよいのである。念のため言っておくが、後進波にするためには

ψ(x) = c\frac{dϕ}{dx}  (26)

とすればよい。

6 練習問題 (進行波)

前回の問題の弦の上を進む進行波と後進波について、計算してみよう。弦の長さ L = 1、波の速度は c = 1、両端は固定されているとする。その条件のもとで、以下 2 つの波が衝突する様子を計算せよ。

後進波は、

\begin{matrix}ϕ_{1}(x) & = a exp [−\frac{(x − x_{0})^{2}}{σ^{2}} ] \\ ただし、 \\ a = 0.5 \\ σ_{1} = 0.02 \\ x_{1} = 0.8 \end{matrix}  (27)

とする。そして、進行波は、

\begin{matrix}ϕ_{2}(x) & = a(x − x_{0}) exp (−\frac{(x − x_{0})^{2}}{σ^{2}} ) \\ ただし、 \\ a = 20 \\ σ_{1} = 0.02 \\ x_{0} = 0.4 \end{matrix}  (28)

とする。 t = 0 のときの様子を、図 5 に示す。波の衝突で、その形はどうなるか?。壁に衝突するとどうなるか?計算せよ。

図 5: t = 0 の時の進行波と後進波の様子。

参考文献^{[1] 高橋大輔. 数値計算. 岩波書店, 1996.}

10