2009年久しぶりの見直し: 疎行列対応万歳

何年か放置していた間に、 Octave にも sparse() などの疎行列関係の関数が用意されるようになっ たので、少し勉強してプログラムをちょっと改良した。

sparse_poisson.m

% sparse_poisson.m --- 正方形領域における Poisson 方程式 (2009/12/29)
function [x,y,u]=sparse_poisson(n)
  h=1/n;
  J=sparse(diag(ones(n-2,1),1)+diag(ones(n-2,1),-1));
  I=speye(n-1,n-1);
  A=-4*kron(I,I)+kron(J,I)+kron(I,J);
  b=-h*h*ones((n-1)*(n-1),1);
% 2次元化を少し工夫
  U=zeros(n-1,n-1);
  U(:)=A\b;
  u=zeros(n+1,n+1);
  u(2:n,2:n)=U;
  
  x=0:1/n:1;
  y=x;
% まず鳥瞰図
  subplot(1,2,1);
  colormap hsv
  mesh(x,y,u);
  colorbar
% 等高線
  subplot(1,2,2);
  contour(x,y,u);
% PostScript を出力
  disp('saving graphs');
  print -depsc2 sparsepoisson.eps

計算内容は同じであるが、$n=100$ くらいスイスイ解けるようになった。 すばらしい。

• ここで subplot() を使って、グラフの鳥瞰図と等高線を並べて描いてある。 subplot(m,n,k) はウィンドウを m 行 n 列の領域に分けて、 その k 番目の領域に描画することを意味する。
• print -depsc2 ファイル名 あるいは print('-depsc2', ファイル名を表す文字列) で、 EPS 形式で出力する (出力サイズが大きい？！)。
• print -dpdf ファイル名 あるいは print('-dpdf', ファイル名を表す文字列) で、 PDF 形式で出力する。
• colorbar でカラーバーを出している。
• subplot(1,2,2); contour(x,y,u); のところは次のようにすべきだった (2012年加筆)。

    right=subplot(1,2,2);
    contour(x,y,u);
    pbaspect(right,[1 1 1]);
  

次に示すのは、solve_poisson.m の疎行列対応版である。

sparse_poisson3.m (2010/2/14)

% sparse_poisson3.m (2010/2/14) かなり効率が高い。描画の方に時間がかかる。
function sparse_poisson3(nx,ny)
% 係数行列の作成 (nx=ny=100で0.031秒)
     hx=1/nx;
     hy=1/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(Ly,speye(m,m))+kron(speye(n,n),Lx);
% 連立方程式を作成して解く (結果は2次元配列に格納, nx=ny=100で0.3秒)
     f=ones(m*n,1);
     U=zeros(m,n);
     U(:)=A\f;
% 境界値0をつける
     u=zeros(nx+1,ny+1);
     u(2:nx,2:ny)=U;
%
     x=0:1/nx:1;
     y=0:1/ny:1;
%
     clf
     colormap hsv
     subplot(1,2,1);
     mesh(y,x,u);
     colorbar
%
     subplot(1,2,2);
     contour(y,x,u);
%
     disp('saving graphs');
     print -depsc2 sparse_poisson3.eps