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

2.1 差分方程式

二次元のラプラス方程式を数値計算で近似解を求めることを考える.まずは,いつものよう に,解$ \phi(x,y)$をテイラー展開する.$ x$$ y$方向に微小変位$ \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$ (3)
$\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$ (4)

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

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

が得られる.このことから,2階の偏導関数の値は微小変位$ h$の場所の関数の値を用いて, $ h^2$の精度で近似計算ができることが分かる.すなわち,式(5) の右辺の第1項を計算すればよいのである.同じことを$ 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)$ (6)

が得られる.

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

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

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

実際にこの式を数値計算する場合,例えば図1のポテンシャルを求め る時には,図5のように格子状2に区切り,その交点での値を求めることになる.ここでは, $ x$および$ y$方向には等間隔$ h$ で区切り計算を進めるが,等間隔である必要はない.そ のようにしても差分の方程式は複雑になるが,同じような計算は可能である.これまでの説明が理解で きていれば,$ x$$ y$方向の間隔が異なっても,式(7)に対応する差分の式が作 れるはずである.

図 4: 解くべき領域を格子に分割
\includegraphics[keepaspectratio, scale=1.0]{figure/lattice.eps}

数値計算をする場合,$ \phi(x,y)$ $ \phi(x+h,y)$の形は不便なので,形式を 改める.図4の左下の座標を(0,0)として,格子点で のポテンシャルを

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

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

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

となり,数値計算し易い形になる.このようにした場合の各格子点の様子を図 5に示す.

次の節で述べる境界条件を考えないとすると,ラプラス方程式は式(9)の 連立方程式を解くだけである.格子に領域を分割することにより,難しげな偏微分方程式 が連立方程式に還元されたわけである.

図 5: 差分の格子
\includegraphics[keepaspectratio, scale=1.0]{figure/sabun.eps}

2.2 境界条件(外周)

実際に,連立方程式(9)を計算する場合,困った問題が生じる.このまま だと,式の数と未知数の数が合わないのである.たとえば,図6 に示す境界を考える.すると,境界が式9$ (i,\,j)$になる場合,式が 作れないのである.すると,図5の中の電極が無い場合,可能 な連立方程式の数は, $ (N_x-1)\times(N_y-1)$個である.未知数の数は, $ (N_x+1)\times(N_y+1)$個である.未知数の数の方が, $ 2(N_x+N_y)$個多いのである.そ のため,予め,この余分の未知数となっている値を決めなくてはならない.実際これは, 偏微分方程式の境界条件を決めることに相当する.
図 6: 境界付近の差分方程式の可能性
\includegraphics[keepaspectratio, scale=1.0]{figure/eq_boundary.eps}

そこで,境界上の格子点 $ 2(N_x+N_y)$個の値を予め決める.こうすれば,式の数を減らさ ないで,未知数の数を減らすことができる.要するに偏微分方程式を解くときの境界条件 を決めるのと同じ.驚いたことに,外周部(境界条件)の格子点の数が,ちょうど,不足し ている方程式の数と同じなのである.自然は,都合良くできているのである.

懸命な諸君であれば,予め決める値は外周の境界上の格子点でなくても良いと考えるだろ う.しかし,内部の点の値を決めてしまうと,連立方程式が1個減ってしまうので,未知 数と式の数の差は変わらない.これについては,良い説明が思い浮かばなかったので,そ ういうものだと思ってください.

2.3 境界条件(内部の電極)

先ほどの説明通り内部の格子点のポテンシャルを決めてしまうと,その数だけ方程式が減 少する.したがって,必要なだけ内部のポテンシャルを決めても,式と未知数の数は同 じで,連立方程式は解ける.そのような理由で,図5 の電極内部の格子点のポテンシャル$ U_{i\,j}$ は,問題で与えられたとおり,30Vと-20V と固定しても,連立方程式の数は変わらない.
ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年1月15日


no counter