B..3.3 MATLABプログラム

MATLAB は定評のあるシミュレーション・ソフトウェアである (桂田 [6])。

明治大学ではライセンス契約をしているので、 学生は申請すれば使えるようになる (MATLAB TAHライセンス)。

(もちろん研究テーマによるが) 卒業研究等で有効に利用できる可能性があり、学ぶ価値は高い。 この授業で MATLAB の解説をする余裕はないが、 ここでは既に MATLAB を知っている人向けにサンプル・プログラムを紹介しておく。 (なお、熱方程式を解くための MATLAB プログラムについては、 桂田 [7], [8] を見よ。)

(B.16) の 係数行列 $A$ は次のプログラムで計算出来る。

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

$$W=3,\quad H=2,\quad f\equiv 1,\quad N_x=30, \quad N_y=20$$

の場合、次のプログラムで計算出来る。

poisson2d_f1.m

% 長方形領域 (0,W)×(0,H) で Poisson 方程式の同次Dirichlet境界値問題を解く
W=3.0;
H=2.0;
nx=30;
ny=20;
m=nx-1;
n=ny-1;
% 連立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); % これについては meshgrid() の説明を見よ。
% f≡1 の場合
F=ones(m*n,1); % ここを頑張ると一般の非同次項の問題が解ける。
%
u=zeros(ny+1,nx+1);
u(2:ny,2:nx)=reshape(A\F,ny-1,nx-1);
%
% グラフの鳥瞰図
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 で分かる

以下の図は、 $f(x,y)=-2(x^2-3x+y^2-2y)$ の場合 (厳密解は $u(x,y)=x(3-x)y(2-y)$) の解を可視化したものである。

図: $-\Laplacian u=-2(x^2-3x+y^2-2y)$