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

3.1 ラプラス方程式

2次元のラプラス方程式

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

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

$\displaystyle \phi(x+h,y)$ $\displaystyle =\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$ (13)
$\displaystyle \phi(x-h,y)$ $\displaystyle =\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$ (14)

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

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

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

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

が得られる.

これらの式(15)と(16)を元の2次元ラプラス 方程式(12)に代入すれば,

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

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

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

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

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

とする.このようにすると,式(17)は

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

となり,数値計算し易い形になる.

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

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

3.2 波動方程式

弦の長さが1,そこを伝わる波の速度を1として,弦の横波の様子を数値計算で解くことを 考える.1次元波動方程式をコンピューターを用いて近似計算するのである.計算に移る 前に,解くべき方程式と条件をきちんと書いておく.解くべき方程式と条件は,

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

となる.この最初の式は波動方程式,2番目を初期条件,3番目を境界条件と言う.

波動方程式を解くためには,初期条件と境界条件が必要である.ある時刻の位置と速度が 決まれば,それ以降を力学的状態は決まってしまう--ということに対応している.振動 の場合は,これに加えて更に,振動の境界条件を決める必要がある.これらが決まって初 めて,振動の状態--ある時刻の変位と速度--が決まるわけである.図3に 初期条件と境界条件の様子を示す.

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

まずは,波動方程式を差分方程式に書き直すことからはじめる.これも,いつものように, 解$ 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*} (21)

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

$\displaystyle \left.\frac{\partial^2 u}{\partial x^2}\right\vert _{x,t}$ $\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)$ (22)

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

が得られる.

これらの式(22)と(23)を元の波動 方程式(20)に代入すれば,

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

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

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

実際に式(25)を数値計算する場合,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*} (26)

と表現を改める.このようにすると,式(25)は

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

となり,数値計算し易い形になる.ただし,

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

である.
ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成20年1月31日


no counter