6 正方形領域での Poisson 方程式の Dirichlet 境界値問題

正方形領域 $\Omega=(0,1)\times(0,1)$ における Poisson 方程式の Dirichlet 境界値問題

$$-\Laplacian u=f$$   $$\mbox{in $\Omega$}$$$$,\quad u=0$$   $$\mbox{on $\rd\Omega$}$$

を差分法で解くには、.... 差分法でどういう連立1次方程式が導かれるかについては、 『Poisson方程式に対する差分解』 を見よ。

$f\equiv 1$ の場合に解くには、次のようなプログラムを用いれば良い。

report_poisson.m

% report_poisson.m
%   正方形領域で -△u=f, 同次Dirichlet境界値問題を解く
%   正方形の各辺を n 等分
function u = report_poisson(n)
  h=1/n;
  J=diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
  Id=eye(n-1,n-1);
  A=-4*kron(Id,Id)+kron(J,Id)+kron(Id,J);
  % b_{ij}=-h^2 f(x_i,y_j)
  b=-h*h*ones((n-1)*(n-1),1);
  U=A\b;
  u=zeros(n+1,n+1);
  for j=1:n-1
    u(2:n,j+1)=U((j-1)*(n-1)+1:j*(n-1));
  end
%endfunction

このプログラムは、効率について無頓着な書き方をされているので、 $n=50$ 程度しか解けない ($n=100$ とすると、未知数の個数は約 $100$ 万で、 小さいコンピューターで解くのは困難である)。

oyabun% octave
GNU Octave, version 2.0.16 (i386-unknown-freebsd3.4).
Copyright (C) 1996, 1997, 1998, 1999, 2000 John W. Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.

octave:1> n=10; u=report_poisson(n);
octave:2> x=0:1/n:1; y=x; contour(u,10,x,y);
octave:3> mesh(x,y,u);
octave:4> n=20; tic; u=report_poisson(n); toc
ans = 1.1213
octave:5> x=0:1/n:1; y=x; contour(u,10,x,y);
octave:6> mesh(x,y,u);
octave:7> n=40; tic; u=report_poisson(n); toc
ans = 16.870
octave:8> x=0:1/n:1; y=x; contour(u,10,x,y);
octave:9> quit
oyabun% 

係数行列の表現を、少し一般化して

$$I_n\otimes\left[\frac{1}{h_x^2}(2I_m-J_m)\right] + \left[\frac{1}{h_y^2}(2I_n-J_n)\right]\otimes I_m,\quad m:=N_x-1,\quad n:=N_y-2$$

と変更した (やはり http://nalab.mind.meiji.ac.jp/~mk/labo/text/poisson.pdf を見よ)。 それを元に書いたのが次の solve_poisson.m である。

solve_poisson.m

% solve_poisson.m --- 正方形領域における Poisson方程式 (2010/2/11)
function solve_poisson(nx,ny)
     hx=1/nx;
     hy=1/ny;
     m=nx-1;
     n=ny-1;
     Im=eye(m,m);
     In=eye(n,n);
     Jm=diag(ones(m-1,1),1)+diag(ones(m-1,1),-1);
     Jn=diag(ones(n-1,1),1)+diag(ones(n-1,1),-1);
     Lx=(2*Im-Jm)/(hx*hx);
     Ly=(2*In-Jn)/(hy*hy);
     A=kron(Ly,Im)+kron(In,Lx);
     f=ones(m*n,1);
     U=zeros(m,n);
     U(:)=A\f;
%
     u=zeros(nx+1,ny+1);
     u(2:nx,2:ny)=U;
%
     x=0:1/nx:1;
     y=0:1/ny:1;
%
     colormap hsv
     subplot(1,2,1);
     mesh(x,y,u);
     colorbar
%
     subplot(1,2,2);
     contour(x,y,u);
%
     disp('saving graphs');
     print -depsc2 solve_poisson.eps

図: 正方形領域で $-\Laplacian u=1$ を解く