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

5.1 ラプラス方程式

2次元のラプラス方程式

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

を数値計で解くことを考える．まずは，いつものように，解$\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$$ (26)

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

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

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

が得られる．このことから，2階の偏導関数の値は微小変位$h$ の場所の関数の値を用いて， $h^2$ の精度で近似計算ができることが分かる．すなわち，式(28) の$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)$$ (29)

が得られる．

これらの式(28)と(29)を元の2次元ラプラス 方程式(25)に代入すれば，

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

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

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

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

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

とする．このようにすると，式(30)は

$$U_{i+1\,j}+U_{i-1\,j}+U_{i\,j+1}+U_{i,j-1}-4U_{i\,j}=0$$ (32)

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

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

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

5.2 波動方程式

弦の長さが1，そこを伝わる波の速度を1として，弦の横波の様子を数値計算で解くことを 考える．1次元波動方程式を数値計で解くことを考える．計算に移る前に，解くべき方程 式と条件をきちんと書いておく．解くべき方程式と条件は，

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

となる．弦を伝わる波の速度は1，弦の長さも1としている．この最初の式は波動方程式で あるが，2番目を初期条件，3番目を境界条件と言う．

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

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

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

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

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

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

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

が得られる．

これらの式(35)と(36)を元の波動 方程式(33)に代入すれば，

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

となる．これが，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]$$ (38)

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

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

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

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

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

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

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

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