B..2.2 2次元の場合

長方形領域 $ \Omega=(0,W)\times (0,H)$ を横方向に $ N_x$ 等分、 縦方向に $ N_y$ 等分する。

$\displaystyle h_x=\frac{W}{N_x},\quad h_y=\frac{H}{N_y},
$

$\displaystyle x_i=i h_x$   ( $ i=0,1,\cdots,N_x$)$\displaystyle ,$

$\displaystyle y_j=j h_y$   ( $ j=0,1,\cdots,N_y$)$\displaystyle ,$

$\displaystyle u_{ij}=u(x_i,y_j)$   ( $ i=0,1,\cdots,N_x$; $ j=0,1,\cdots,N_y$)$\displaystyle .
$

$ u_{ij}$ の近似値 $ U_{ij} $ を求めることが目標となることが多い。

例えば Poisson 方程式のDirichlet 境界値問題

$\displaystyle -\Laplacian u=f$   in $ \Omega$$\displaystyle ,\quad u=\phi$   on $ \rd\Omega$

を解く場合、 $ U_{ij} $ ( $ 1\le i\le N_x-1$, $ 1\le j\le N_y-1$) が未知数となる。

$\displaystyle m:=N_x-1,\quad n:=N_y-1
$

とおく。

連立1次方程式を行列とベクトルで表示するためには、 $ U_{ij} $ を1次元的に並べたベクトルを作る必要がある。 例えば

$\displaystyle \bm{U}=
(U_{1,1},U_{2,1},\cdots,U_{m,1},
U_{1,2},U_{2,2},\cdots...
...U_{1,3},U_{2,3},\cdots,U_{m,3},
\cdots,
U_{1,n},U_{2,n},\cdots,U_{m,n}
)^T.
$

$ \bm{U}$ の成分を $ U_\ell$ と書くとき

$\displaystyle U_{i,j}=U_{i+n(j-1)}
$

という式が成り立つことが分かる。

$\displaystyle \ell=i+n(j-1)$ (6)

とおくと

$\displaystyle \ell-1=(i-1)+n(j-1),\quad 0\le i-1\le n-1
$

が成り立つので、$ i-1$$ \ell-1$$ n$ で割った余り、 $ j-1$$ \ell-1$$ n$ で割った商である。ゆえに

$\displaystyle U_\ell=U_{i,j},\quad j=\lfloor (\ell-1)/n \rfloor+1,\quad i=\mathrm{mod}(\ell-1,n)+1.$ (7)


\begin{yodan}
例えば C 言語のプログラムならば、$\ell$\ から $i$...
...j の式;
}
}\end{verbatim}
\end{screen}とすれば良い。 \qed
\end{yodan}

Laplacian $ \Laplacian=\frac{\rd^2}{\rd x^2}+\frac{\rd^2}{\rd y^2}$ には、 2階微分が現れるが、それを2階中心差分近似すると、 $ -\Laplacian u=f$

$\displaystyle -\left[
\frac{U_{i+1,j}-2U_{i,j}+U_{i-1,j}}{h_x^2}
+\frac{U_{i,j+1}-2U_{i,j}+U_{i,j-1}}{h_y^2}
\right]
=f(x_i,y_j)
$

という差分方程式で置き換えられる。

  $\displaystyle A\bm{U}=\bm{f},$ (8)
  $\displaystyle A= I_n\otimes\left[\frac{1}{h_x^2}\left(2I_m-J_m\right)\right] +\left[\frac{1}{h_y^2}\left(2I_n-J_n\right)\right]\otimes I_m,$ (9)
  $\displaystyle \bm{f}=(f_\ell),\quad f_\ell=f(x_i,y_j).$ (10)

% W, H, Nx, Ny に値が適当に記憶されているとして
% 例: W=2; H=1; Nx=8; Ny=4;
hx=W/Nx;
hy=H/Ny;
m=Nx-1;
n=Ny-1;
ex=ones(m,1);
ax=spdiags([-ex 2*ex -ex],-1:1,m,m)/(hx*hx);
ey=ones(n,1);
ay=spdiags([-ey 2*ey -ey],-1:1,n,n)/(hy*hy);
Im=speye(m,m);
In=speye(n,n);
a=kron(ay,Im)+kron(In,ax);

なお、$ U_{ij} $ を並べる際に、(6) でなく

$\displaystyle U_\ell=U_{i,j},\quad \ell=m(i-1)+j$ (11)

とする場合もある。むしろそちらの方が多いかもしれない。 MATLAB で meshgrid() などを使うには、こちらの方が都合が良い (本文の例を参照せよ)。 この場合は

$\displaystyle j=\mathrm{mod}(\ell-1,m)+1,\quad i=\lfloor(\ell-1)/m\rfloor+1.$ (12)

また連立1次方程式の係数行列は

$\displaystyle A=I_m\otimes\left[\frac{1}{h_y^2}(2I_n-J_n)\right] +\left[\frac{1}{h_x^2}(2I_m-J_m)\right]\otimes I_n.$ (13)

コードの方は最後の a=kron(ay,Im)+kron(In,ax);
a=kron(Im,ay)+kron(ax,In);
に置き換えれば良い。

桂田 祐史