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

2.1 差分方程式

2次元のラプラス方程式を数値計で解くことを考える。まずは、いつものよう に、解$ \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$ で区切り計算 を進めるが、等間隔である必要はない。多少、式(7) は異なる が同じような計算は可能である。これまでの説明が理解できていれば、xとy方 向の間隔が異なっても、式(7)に対応する差分の式が作れるはず である。

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

数値計算をする場合、$ \phi(x,y)$ $ \phi(x+h,y)$の形は不便なので、形式を 改める。図5の左下の座標を(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-h}-4U_{i\,j}=0$ (9)

となり、数値計算し易い形になる。各格子点の様子を図 5に示す。

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

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

2.2 境界条件(外周)

実際に、連立方程式(6)を計算する場合、困った問題が生じる。 このままだと、式の数と未知数の数が合わないのである。たとえば、図 [*]に示す境界を考える。すると、境界が式 6$ (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年8月21日


no counter