3.1 Dirichlet境界値問題

まずは同次Dirichlet境界値問題のプログラムを示す。

長方形領域 $\Omega=(a,b)\times(c,d)$ における Poisson 方程式の境界値問題

$$-\Laplacian u=f$$   in $\Omega$$$,\quad u=0$$   (on $\rd\Omega$)

を考える。

$$h_x=\frac{b-a}{N_x},\quad h_y=\frac{d-c}{N_y},$$

$$x_i=a+i h_x i=0,1,\dots,N_x ,\quad y_j=c+j h_y j=0,1,\dots,N_y ,$$

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

とおき、$u_{ij}$ の近似値を $U_{ij}$を求めるための差分方程式は

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

$$1\le i\le N_x-1 1\le j\le N_y-1 ,$$

$$U_{0,j}=U_{N_x,j}=0 0\le j\le N_y ,$$ (2)

$$U_{i,0}=U_{0,N_y}=0 0\le i\le N_x$$ (3)

Dirichlet 境界値問題については、 $U_{ij}$ ( $1\le i\le N_x$, $1\le j\le N_y$) が未知数になる。 桂田 [2] というノートでは

$$U_\ell=U_{i,j},\quad \ell= i+j(N_x-1) 1\le i\le N_x-1 1\le j\le N_y-1$$ (4)

とおき、 $\bm{U}=\begin{pmatrix}U_1\\ U_2\\ \vdots \\ U_{N}\end{pmatrix}$ を未知ベクトルとする連立1次方程式を考えた。 連立1次方程式の係数行列は

$$A=\frac{1}{h_y^2}(2I_{N_y-1}-J_{N_y-1})\otimes I_{N_x-1} +I_{N_y-1}\otimes \frac{1}{h_x^2}(2I_{N_x-1}-J_{N_x-1})$$

であった。

以下のプログラムでは (4) の代わりに

$$U_\ell=U_{i,j},\quad \ell= j+i(N_y-1) 1\le i\le N_x-1 1\le j\le N_y-1$$ (5)

とおく。この場合の係数行列は

$$A=I_{N_x-1}\otimes\frac{1}{h_y^2}(2I_{N_y-1}-J_{N_y-1}) +\frac{1}{h_x^2}(2I_{N_x-1}-J_{N_x-1})\otimes I_{N_y-1}.$$

poisson_coef.m

function A=poisson_coef(W, H, nx, ny)
% 長方形領域 (0,W)×(0,H) における Poisson 方程式の Dirichlet 境界値問題
% Laplacian を差分近似した行列を求める。
% 長方形を nx×ny 個の格子に分割して差分近似する。
% MATLABでは
%  (1) 行列は Fotran と同様の column first であり、
%  (2) mesh(), contour() による「行列描画」は Z(j,i) と添字の順が普通と逆なので、
% l=i+(j-1)*(nx-1) と row first となるように1次元的番号付けする
  hx=W/nx;
  hy=H/ny;
  m=nx-1;
  n=ny-1;
  ex=ones(nx,1);
  ey=ones(ny,1);
  Lx=spdiags([-ex,2*ex,-ex],-1:1,m,m)/(hx*hx);
  Ly=spdiags([-ey,2*ey,-ey],-1:1,n,n)/(hy*hy);
  A=kron(speye(m,m),Ly)+kron(Lx,speye(n,n));

poisson2d.m

% 長方形領域 (0,W)×(0,H) で Poisson 方程式の同次Dirichlet境界値問題を解く
W=3.0;
H=2.0;
nx=30;
ny=20;
m=nx-1;
n=ny-1;
% 連立方程式を作成して解く
% MATLABの行列は Fotran と同様の column first であり、
%「行列描画」は Z(j,i) と添字の順が普通と逆なので、
% l=i+(j-1)*(nx-1) と row first となるように1次元的番号付けする
A=poisson_coef(W, H, nx, ny);
%
x=linspace(0,W,nx+1); % x=[x_0,x_1,...,x_nx]
y=linspace(0,H,ny+1); % y=[y_0,y_1,...,y_ny]
[X,Y]=meshgrid(x,y);
% f≡1 の場合
%F=ones(m*n,1);
f=-2*(X.^2-3*X+Y.^2-2*Y);
f=f(2:ny,2:nx);
F=f(:);
%
U=zeros(n,m);
U(:)=A\F;
% 境界値0をつける
u=zeros(ny+1,nx+1);
u(2:ny,2:nx)=U;
%
% グラフの鳥瞰図
clf
colormap hsv
subplot(1,2,1);
mesh(X,Y,u);
colorbar
% 等高線
subplot(1,2,2);
contour(X,Y,u);
%
disp('図を保存する');
print -dpdf poisson2d.pdf % 利用できるフォーマットは doc print で分かる
print -dpng poisson2d.png % 利用できるフォーマットは doc print で分かる
print -deps poisson2d.eps % 利用できるフォーマットは doc print で分かる

図 1: poisson2d.m の結果

まったく別の時期に作ったプログラム。 ほとんど同じで我ながら唖然とする。 最初のプログラムが mesh() と contour() で、 鳥瞰図と等高線を別々に描いたが、 こちらは meshc() で同時に描いている。

poisson2d_v2.m

% poisson2d.m --- Poisson equation -△u=sin(π x)sin(2π y) (0<x<3, 0<y<2), u=0
% [0,3]×[0,2]を 30×20 に分割する
a=0; b=3; c=0; d=2;
nx=30; ny=20;
%nx=6; ny=4;
X=linspace(a,b,nx+1);
Y=linspace(c,d,ny+1);
[x,y]=meshgrid(X,Y);
% f(x,y)=sin(x)sin(2y)
F=sin(pi*x).*sin(2*pi*y);
% 係数行列
hx=(b-a)/nx; hy=(d-c)/ny;
e=ones(nx-1,1);
ax=spdiags([-e 2*e -e],-1:1,nx-1,nx-1)/(hx*hx);
e=ones(ny-1,1);
ay=spdiags([-e 2*e -e],-1:1,ny-1,ny-1)/(hy*hy);
a=kron(speye(nx-1),ay)+kron(ax,speye(ny-1));
% F の境界部分の値を削除して、1次元化
f=F(2:end-1,2:end-1);
f=f(:);
%
u=zeros(ny+1,nx+1);
u(2:end-1,2:end-1)=reshape(a\f,ny-1,nx-1);
%
figure('Name','Poisson equation -△u=sin(π x)sin(2π y) (0<x<3, 0<y<2)')
meshc(x,y,u)
figure(gcf)

図 2: poisson2d_v2.m の結果