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

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

2.1 差分方程式

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

$$\left\{ \begin{aligned}% &\frac{\partial^2 u}{\partial x^2}= \f... ...& (0\leqq x \leqq 1)\\ % &u(0,t)=u(1,t)=0 % \end{aligned} \right.$$

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

波動方程式の他に，初期条件と境界条件がある．力学的状態は，ある時刻，ここでは $t=0$ の時の変位とその変位の速度が決まれば，それ以降を決めることができる．振動の 場合は，これに加えて更に，振動の境界条件を決める必要がある．これらが決まって初め て，波動方程式とともに，振動の状態，ある時刻と位置の変位の値が決まるわけである． 図4に初期条件と境界条件の様子を示す．

図 2: 時刻$t=0$ のときの弦の様子(スナップショット)．初期条件と境界 条件が表されており，y方向の速度が$\psi (x)$ になっている．

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

$$\begin{aligned}u(x+\Delta x,t)&=u(x,t) +\frac{\partial u}{\part... ...rac{\partial^4 u}{\partial x^4}(\Delta x)^4 -\cdots \end{aligned}$$

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

$$\left.\frac{\partial^2 u}{\partial x^2}\right\vert _{x,t} =\frac{1}{\Delta x^2}\left[ u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)\right]-O(\Delta x^2)$$ (6)

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

$$\left.\frac{\partial^2 u}{\partial t^2}\right\vert _{x,t} =\frac{1}{\Delta t^2}\left[ u(x,t+\Delta t)-2u(x,t)+u(x,t-\Delta t)\right]-O(\Delta t^2)$$ (7)

が得られる．

これらの式(6)と(7)を元の波動 方程式(4)に代入すれば，

$$\frac{1}{\Delta x^2}\left[u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)\right]= \frac{1}{\Delta t^2}\left[u(x,t+\Delta t)-2u(x,t)+u(x,t-\Delta t)\right]$$ (8)

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

$$u(x,t+\Delta t)= 2u(x,t)-u(x,t-\Delta t)+ \frac{\Delta t^2}{\Delta x^2}\left[ u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t) \right]$$ (9)

とすることができる．この式の右辺は，時刻$t$ と $t-\Delta t$ の値でである．そして， 左辺は時刻 $t+\Delta t$ の値である．このことから，式(9)を用いると，時 刻$t$ と $t-\Delta t$ の値から， $t+\Delta t$ の値が計算できることになる．

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

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

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

$$u_{i\,j+1}= 2u_{i\,j}-u_{i\,j-1}+ \alpha\left(u_{i+1\,j}-2u_{i\,j}+u_{i-1\,j}\right)$$ (11)

となり，数値計算し易い形になる．ただし，

$$\alpha = \frac{\Delta t^2}{\Delta x^2}$$ (12)

である．

この式を用いた計算の様子を図3に示す．

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

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

2.2 初期条件

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

$$\begin{aligned}% &u_{1\,0} && u_{2\,0} && u_{3\,0} && u_{4\,0}... ...N_t} && u_{5\,N_t} && \cdots && u_{N_x-1\,N_t}\\ % \end{aligned}$$

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

$$u_{0\,j}=u_{N_x\,j}=0 \qquad (j=0,\,1,\,2,\,3,\cdots,N_t)$$ (13)

を，忘れてはならない．

これを計算するためには，まず， $u_{i\,0}\quad(i=1,2,3,\cdots,N_x-1)$ の値を決める 必要がある．これ以前の状態が分からないので，式(11)は使えないが，式 (4)の初期条件が使える．すなわち，

$$u_{i\,0}=\phi(i\Delta x)$$ (14)

である．

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

$$\begin{aligned}% u(x,\Delta t)&=u(x,0) +\frac{\partial u}{\par... ...^2 u}{\partial x^2}(\Delta t)^2 +O(\Delta x^3)\\ % \end{aligned}$$

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

$$\begin{aligned}% u(x,\Delta t) % &\thickapprox u(x,0) +\psi(x... ...[ u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)\right] % \end{aligned}$$

となる．これは，めでたい式である．右辺は，$t=0$ のみの値で構成されている．これで， $u_{i\,1}\quad(i=1,2,3,\cdots,N_x-1)$ が計算可能になった．この式から，

$$u_{i\,1}=u_{i\,0} +\psi(x_i)\Delta t +\frac{\alpha}{2}\left[ u_{i+1\,0}-2u_{i\,0}+u_{i-1\,0}\right]$$ (17)

が得られる．

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

2.3 進行波の取り扱い

今までの議論で定在波の取り扱いは可能であろう．そこで，進行波の記述方法について， コメントしておく．進行波を数値計算すると面白いのでその方法を示す．進行波を記述す るためには，初期条件さえ記述すれば，後の差分方程式は同じである．その初期条件の記 述の仕方を示す．

元の波動方程式

$$\frac{\partial^2 u}{\partial x^2}= \frac{1}{c^2}\frac{\partial^2 u}{\partial t^2}$$ (18)

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

$$u(x,t)=f(x-ct)+g(x+ct)$$ (19)

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

初期条件

$$u(x,0)=\phi(x)$$ (20)

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

$$\frac{\partial u}{\partial t}(x,0)=\psi(x)$$ (21)

である．進行波になるように，$\psi (x)$ を決めればよい．$u(x,t)$ を進行波と仮定する と，式(20)から

$$u(x,\Delta t)=\phi(x-c\Delta t)$$ (22)

となる．この式を使って，$\psi (x)$ を求めることにする．$\psi (x)$ の定義より，

$$\begin{aligned}\psi(x) &=\frac{\partial u}{\partial t}(x,0)\\ &... ...-\phi(x-\Delta x)}{\Delta x}\\ &=-c\frac{d\phi}{dx} \end{aligned}$$

となる．進行波にするためには，$\psi (x)$ は$\phi(x)$ の導関数ににすればよいのである．

念のため言っておくが，後進波にするためには

$$\psi(x)=c\frac{d\phi}{dx}$$ (24)

とすればよい．
ホームページ: Yamamoto's laboratory
著者: 山本昌志