有限差分法 - 数理生物学演習

有限差分法では，微分を差分商で近似することで微分方程式を数値的に解く．ここでは，1階と2階微分の近似としての各種差分商の導出と，それを利用した微分方程式の離散化の例をポアソン方程式を用いて示す．

離散化の手法¶

離散化とは連続な関数を数値的に扱うために，有限で非連続（離散）な値に分割することをいう．離散化の方法は様々だが，目的や条件をにより使い分ける．例えば，空間方向の離散化には

有限差分法
有限体積法
有限要素法
粒子法

などがある．

これらは空間的に連続関数の離散化に用いられる手法であり，時間的に連続な関数の離散化には差分法が用いられる．ここでは有限差分法について解説する．

有限差分法 finite difference method¶

有限差分法(FDM)は，物体を格子点の集まりとみなし，テーラー展開に基づき微分を差分商で近似して，差分方程式を作り，微分方程式を解く．

差分 difference¶

関数が2つの変数値に対してとる値の有限な差を差分(difference)という．この差分を変数値の差(刻み幅)で割ったものを差分商(difference quotient)という．

ある2次元空間上の関数 $u(x,y)$ の離散化を考える． 2次元空間を格子で区切り，各格子点の座標を $(x_i,y_j)$ ，格子間隔はそれぞれ $\Delta x, \Delta y$ とする．格子番号を用いて， $u_{i,j}=u(x_i,y_j)$ とすれば， $u_{i+1,j}, u_{i-1,j}$ は $(x_i, y_j)$ のまわりで

u_{i+1,j}=u(x_i+\Delta x, y_j)=u_{i,j}+\left( \frac{\partial u}{\partial x}\right)_i \Delta x +\frac{1}{2!} \left( \frac{\partial^2 u}{\partial x^2}\right)_i \Delta x^2 +\frac{1}{3!} \left( \frac{\partial^3 u}{\partial x^3}\right)_i \Delta x^3 +\dots

u_{i-1,j}=u(x_i-\Delta x, y_j)=u_{i,j}-\left( \frac{\partial u}{\partial x}\right)_i \Delta x +\frac{1}{2!} \left( \frac{\partial^2 u}{\partial x^2}\right)_i \Delta x^2 -\frac{1}{3!} \left( \frac{\partial^3 u}{\partial x^3}\right)_i \Delta x^3 +\dots

とそれぞれテイラー展開できる．

それぞれ2次以上の項を無視すれば，

u_{i+1,j}=u_{i,j}+\left( \frac{\partial u}{\partial x}\right)_i \Delta x+O(\Delta x^2) \approx u_{i,j}+\left( \frac{\partial u}{\partial x}\right)_i \Delta x \\ u_{i-1,j}=u_{i,j}-\left( \frac{\partial u}{\partial x}\right)_i \Delta x+O(\Delta x^2) \approx u_{i,j}-\left( \frac{\partial u}{\partial x}\right)_i \Delta x

なので，

u_{i+1,j}-u_{i,j}=\left( \frac{\partial u}{\partial x}\right)_i \Delta x \\ u_{i,j}-u_{i-1,j}=\left( \frac{\partial u}{\partial x}\right)_i \Delta x

これを，それぞれ前進差分(forward difference)，後退差分(backward difference)と呼ぶ．

これを用いれば，

\begin{align*} \left( \frac{\partial u}{\partial x}\right)_i&=&\frac{u_{i+1,j}-u_{i,j}}{\Delta x} \\ \left( \frac{\partial u}{\partial x}\right)_i&=&\frac{u_{i,j}-u_{i-1,j}}{\Delta x} \end{align*}

と微分を差分商により $\Delta x$ について1次の精度で近似できる．

また，(1) $-$ (2)より，

\begin{split} u_{i+1,j}-u_{i-1,j}&=2\left( \frac{\partial u}{\partial x}\right)_i \Delta x+O(\Delta x^3)\\ &\approx 2\left( \frac{\partial u}{\partial x}\right)_i \Delta x \end{split}

\left( \frac{\partial u}{\partial x}\right)_i=\frac{u_{i+1,j}-u_{i-1,j}}{2 \Delta x}

と中心(中央)差分(central difference)及びその差分商による $\Delta x$ について2次の精度での近似が求まる．

さらに，(1) $+$ (2)より，

\begin{split} u_{i+1,j}+u_{i-1,j}&=2u_{i,j}+\left( \frac{\partial^2 u}{\partial x^2}\right)_i \Delta x^2+O(\Delta x^4)\\ &\approx 2u_{i,j}+\left( \frac{\partial^2 u}{\partial x^2}\right)_i \Delta x^2 \end{split}

\left( \frac{\partial^2 u}{\partial x^2}\right)_i=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}

と2階(中心)差分(2nd difference)が求まり，2階微分の係数を $\Delta x$ について2次の精度で近似できる．ここでは $x$ 方向についてのみ考えたが， $y$ 方向についても同様の議論が成り立つ．

ポアソン方程式 Poisson’s equation¶

ポアソン方程式とは

\begin{split} \bigtriangleup u&=f\\ \nabla^2 u&=f\\ \frac{\partial^2 u}{\partial x_1^2}+\frac{\partial^2 u}{\partial x_2^2}+\cdots+\frac{\partial^2 u}{\partial x_n^2}&=f(x_1,x_2, \cdots, x_n) \end{split}

と表現される，2階の楕円型偏微分方程式．境界値問題となる．

例えば，静電ポテンシャル $\phi$ は電荷の分布 $\rho$ に対し，真空での誘電率を $\epsilon_0$ とすれば

\bigtriangleup \phi=-\frac{\rho}{\epsilon_0}

となり，ポアソン方程式になる．

ポアソン方程式を例に有限差分法を使ってみる．今回は2次元だけを考える．とすれば，

\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=f(x,y)

空間を格子で区切り各格子点を $(x_i,y_j)$ ，格子幅を $x,y$ 両方向とも $h$ としよう．すると，

\begin{align*} u_{xx}&=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}+O(h^2)\\ u_{yy}&=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}+O(h^2)\\ u_{xx}+u_{yy}&=-\frac{4u_{i,j}-u_{i-1,j}-u_{i,j+1}-u_{i+1,j}-u_{i,j-1}}{h^2}+O(h^2) \end{align*}

よって

\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}= -\frac{4u_{i,j}-u_{i-1,j}-u_{i,j+1}-u_{i+1,j}-u_{i,j-1}}{h^2}+O(h^2)

なので $h$ が十分小さければ $h^2$ の以上の項は無視できるので，

\frac{4u_{i,j}-u_{i-1,j}-u_{i,j+1}-u_{i+1,j}-u_{i,j-1}}{h^2}=-f_{i,j}

と離散化される． $f_{i,j}$ は既知量である．よって，得られた $i \times j$ 個の連立一次方程式をとけば良い．ただし，境界領域での近似の仕方を考える必要がある．