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

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

2.1 差分方程式

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

\begin{equation*}\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.\end{equation*}

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

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

図 2: 時刻$ t=0$のときの弦の様子(スナップショット)。初期条件と境界 条件が表されており、y方向の速度が$ \psi (x)$になっている。
\includegraphics[keepaspectratio, scale=0.85]{figure/wave_init.eps}

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

\begin{equation*}\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}\end{equation*}

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

$\displaystyle \left.\frac{\partial^2 u}{\partial x^2}\right\vert _{x,y}$ $\displaystyle =\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項を計算すればよいのである。ラプラス 方程式と同じである。同様なことを時間軸方向についても行うと

$\displaystyle \left.\frac{\partial^2 u}{\partial t^2}\right\vert _{x,t}$ $\displaystyle =\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)に代入すれば、

$\displaystyle \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次元波動方程式の差分の式である。この式を計算し易いように、もう 少し変形すると、

$\displaystyle 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{equation*}\begin{aligned}u(x,t)&=u(i\Delta x,j\Delta t) &=u_{i j} \end{aligned}\end{equation*}

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

$\displaystyle 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)

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

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

である。

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

図 3: 差分方程式の計算の様子
\includegraphics[keepaspectratio, scale=0.85]{figure/sabun.eps}

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

2.2 初期条件

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

\begin{equation*}\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}\end{equation*}

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

$\displaystyle 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)の初期条件が使える。すなわち、

$\displaystyle 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 x)$なので、とりあえず テーラー展開してみる。これを、$ t=0$の周りでテーラー展開すると、

\begin{equation*}\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}\end{equation*}

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

\begin{equation*}\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}\end{equation*}

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

$\displaystyle 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)に従い、計算すればよい。


ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成17年2月18日


no counter