3.1 Dirichlet境界値問題

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

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

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

を考える。

  $\displaystyle h_x=\frac{b-a}{N_x},\quad h_y=\frac{d-c}{N_y},$    
  $\displaystyle x_i=a+i h_x$   ( $ i=0,1,\dots,N_x$)$\displaystyle ,\quad y_j=c+j h_y$   ( $ j=0,1,\dots,N_y$)$\displaystyle ,$    
  $\displaystyle u_{ij}=u(x_i,y_j)$   ( $ i=0,1,\dots,N_x$; $ j=0,1,\dots,N_y$)    

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

  $\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)$ (1)
     ( $ 1\le i\le N_x-1$, $ 1\le j\le N_y-1$)$\displaystyle ,$    
  $\displaystyle U_{0,j}=U_{N_x,j}=0$   ( $ 0\le j\le N_y$)$\displaystyle ,$ (2)
  $\displaystyle 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] というノートでは

$\displaystyle 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次方程式の係数行列は

$\displaystyle 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) の代わりに

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

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

$\displaystyle 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 の結果
Image poisson2d

まったく別の時期に作ったプログラム。 ほとんど同じで我ながら唖然とする。 最初のプログラムが 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 の結果
Image poisson2d-v2


Subsections
桂田 祐史