B..2.2 2次元の場合

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

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

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

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

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

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

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

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

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

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

とおく。

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

$$\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$ と書くとき

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

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

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

とおくと

$$\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$ で割った商である。ゆえに

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

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

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

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

$$A\bm{U}=\bm{f},$$ (8)

$$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)

$$\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) でなく

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

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

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

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

$$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);

に置き換えれば良い。