5 差分法による偏微分方程式の数値計算

5.1 ラプラス方程式

2次元のラプラス方程式

$$\frac{\partial^2 \phi}{\partial x^2}+ \frac{\partial^2 \phi}{\partial y^2}=0$$ (22)

を数値計で解くことを考える。まずは、いつものように、解$\phi(x,y)$をテイラー展開 する。xおよび、y方向に微小変位$\pm h$があった場合、

$$\phi(x+h,y) =\phi(x,y) +\frac{\partial\phi}{\partial x}h +\frac{1}{2!}\frac{\... ...i}{\partial x^3}h^3 +\frac{1}{4!}\frac{\partial^4\phi}{\partial x^4}h^4 +\cdots$$ (23)

$$\phi(x-h,y) =\phi(x,y) -\frac{\partial\phi}{\partial x}h +\frac{1}{2!}\frac{\... ...i}{\partial x^3}h^3 +\frac{1}{4!}\frac{\partial^4\phi}{\partial x^4}h^4 -\cdots$$ (24)

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

$$\left.\frac{\partial^2\phi}{\partial x^2}\right\vert _{x,y} =\frac{1}{h^2}\left[ \phi(x+h,y)-2\phi(x,y)+\phi(x-h,y)\right]-O(h^2)$$ (25)

が得られる。このことから、2階の偏導関数の値は微小変位$h$の場所の関数の値を用いて、 $h^2$の精度で近似計算ができることが分かる。すなわち、式(27) の$O(h^2)$を除いた右辺を計算すればよいのである。同じことをy方向についても行うと

$$\left.\frac{\partial^2\phi}{\partial y^2}\right\vert _{x,y} =\frac{1}{h^2}\left[ \phi(x,y+h)-2\phi(x,y)+\phi(x,y-h)\right]-O(h^2)$$ (26)

が得られる。

これらの式(27)と(28)を元の2次元ラプラス 方程式(24)に代入すれば、

$$\phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)-4\phi(x,y)=0$$ (27)

となる。これが、2次元ラプラス方程式の差分の式である。この式を眺めると、座標 $(x,y)$のポテンシャルの値$\phi(x,y)$は、周りの値の平均であることがわかる。

実際にこの式を数値計算する場合、計算領域を間隔$h$で格子状³に区切り、その交点での値を求めることになる。 ここでは、xおよびy方向には等間隔$h$ で区切り計算を進めるが、等間隔である必要はな い。多少、式(29) は異なるが同じような計算は可能である。これまでの説 明が理解できていれば、xとy方向の間隔が異なっても、式(29)に対応する差 分の式が作れるはずである。

実際、数値計算をする場合、$\phi(x,y)$や $\phi(x+h,y)$の形は不便なので、形式を改め る。各格子点でのポテンシャルを

$$\begin{aligned}\phi(x,y)&=\phi(ih,jh)\ &=U_{i j} \end{aligned}$$

とする。このようにすると、式(29)は

$$U_{i+1 j}+U_{i-1 j}+U_{i j+1}+U_{i j-1}-4U_{i j}=0$$ (29)

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

ラプラス方程式は式(31)の連立方程式を解くだけである。格子に領域を分 割することにより、難しげな偏微分方程式が連立方程式に還元されたわけである。

連立方程式を解くわけであるが、このままでは、式の数と未知数の数が異なる。格子点で のポテンシャルの値を求めるためには、境界条件を設定する必要がある。それにより、式 の数と未知数の数が同一になり、格子点でのポテンシャルを求めることができる。

5.2 波動方程式

弦の長さが1、そこを伝わる波の速度を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番目を境界条件と言う。

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

図 4: 時刻$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)$$ (32)

が得られる。このことから、2階の偏導関数の値は微小変位$\Delta x$の場所の関数の 値を用いて、 $(\Delta x)^2$の精度で近似計算ができることが分かる。すなわち、式( 34)の右辺の第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)$$ (33)

が得られる。

これらの式(34)と(35)を元の波動 方程式(32)に代入すれば、

$$\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]$$ (34)

となる。これが、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,y) \right]$$ (35)

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

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

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

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

$$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)$$ (37)

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

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

である。
ホームページ: Yamamoto's laboratory
著者: 山本昌志