next up previous contents
Next: 3.2 Neumann境界値問題 Up: 3.1 Dirichlet境界値問題 Previous: 3.1 Dirichlet境界値問題

3.1.1 非同次Neumann境界値問題

非同次の問題の解説 [3] も書いた。

ここでは、 [3] で作成したプログラムを紹介する、

poisson2d_nh.m

% poisson2d_nh.m
% 長方形領域におけるPoisson -△u=f (非同次Dirichlet境界条件) を差分法で解く
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2d_nh.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/poisson2d-nonhomo.pdf
% 一応の完成をして、poisson2d_nh.m として公開した (2024/8/25)。

function poisson2d_nh(Nx,Ny)
  arguments
    Nx=50
    Ny=50
  end
  % Dirichlet境界条件の場合の係数行列
  function A = poisson2d_mat(W,H,Nx,Ny)
    Ix=speye(Nx-1,Nx-1);
    Iy=speye(Ny-1,Ny-1);
    vx=ones(Nx-2,1);
    Jx=sparse(diag(vx,1)+diag(vx,-1));
    vy=ones(Ny-2,1);
    Jy=sparse(diag(vy,1)+diag(vy,-1));
    hx = W / Nx;
    hy = H / Ny;
    betax = 1.0 / hx^2;
    betay = 1.0 / hy^2;
    Kx=betax * (2*Ix - Jx);
    Ky=betay * (2*Iy - Jy);
    % row major (y changes first, l=j+(Ny-1)*(i-1))
    A=kron(Kx,Iy) + kron(Ix,Ky);
    % column major (x chandes first, l=i+(Nx-1)*(j-1))
    %  A=kron(Iy,Kx) + kron(Ky,Ix);
  end
  function A=poisson_coef(W, H, nx, ny)
    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));
  end
  % 長方形領域 (a,b)×(c,d)
  a=0; b=1; c=0; d=1;
  %
  function u = exact_solution(x, y)
    u = cos(pi * x) .* cos(pi * y);
  end
  % Poisson 方程式の右辺 (exact_solution の -△)
  function val = f(x, y)
    val = 2 * pi^2 * cos(pi * x) .* cos(pi * y);
  end
  % 長方形の下の辺での b()
  function val = b_b(x)
    val = exact_solution(x,c);
  end
  % 長方形の上の辺での b()
  function val = b_t(x)
    val = exact_solution(x,d);
  end
  % 長方形の左の辺での b()
  function val = b_l(y)
    val = exact_solution(a, y);
  end
  % 長方形の右の辺での b()
  function val = b_r(y)
    val = exact_solution(b, y);
  end
  % 誤差を測るためのノルム (最大値ノルム)
  function val = supnorm(a)
    [m,n]=size(a);
    val = norm(reshape(a,m*n,1),Inf);
  end
  % 格子への分割
  %Nx=100;
  %Ny=100;
  hx=(b-a)/Nx;
  hy=(d-c)/Ny;
  betax=1.0/hx^2;
  betay=1.0/hy^2;
  disp(sprintf('Nx=%d, Ny=%d, hx=%f, hy=%f, betax=%f, betay=%f', ...
	       Nx, Ny, hx, hy, betax, betay));

  % 差分方程式 A U^{n+1}=B U^n の行列
  A=poisson2d_mat(b-a,d-c,Nx,Ny);
  N=(Nx-1)*(Ny-1);
  if N <= 20
    disp('A='); full(A)
  end
  %A2=poisson_coef(b-a,d-c,Nx,Ny);
  %err=A-A2;
  %disp(sprintf('A-A2=%e', supnorm(err)));
  %clear A2 err
  % LU分解
  [AL,AU,ap]=lu(A,'vector');

  % 格子点の座標ベクトル x=(x_1,x_2,...,x_{Nx+1}), y=(y_1,y_2,...,y_{Ny+1})
  X=linspace(a,b,Nx+1)';
  Y=linspace(c,d,Ny+1)';
  % 格子点のx,y座標の配列 X={X_{ij}}, Y={Y_{ij}}
  [x,y]=meshgrid(X,Y);

  % Poisson方程式の右辺の関数値 (連立1次方程式に必要)
  u=f(x,y);
  meshc(x,y,u); axis([a b c d -50 50]); title('graph of f'); drawnow; shg
  % 計算に用いるため、f の内部格子点での値をベクトル U に格納する
  U=reshape(u(2:Ny,2:Nx),(Nx-1)*(Ny-1),1);

  % 境界値 (左の辺、右の辺)
  blvector=b_l(Y(2:Ny));
  brvector=b_r(Y(2:Ny));
  lindex=1:Ny-1;
  rindex=(Nx-2)*(Ny-1)+1:(Nx-1)*(Ny-1);
  % 境界値 (下の辺、上の辺)
  bbvector=b_b(X(2:Nx));
  btvector=b_t(X(2:Nx));
  bindex=1:Ny-1:(Nx-1)*(Ny-1);
  tindex=Ny-1:Ny-1:(Nx-1)*(Ny-1);

  % 境界条件により分かっている値を右辺に移項
  U(lindex)=U(lindex)+betax*blvector;
  U(rindex)=U(rindex)+betax*brvector;
  U(bindex)=U(bindex)+betay*bbvector;
  U(tindex)=U(tindex)+betay*btvector;

  % 連立1次方程式を解いて差分解を求める
  U=AU\(AL\U(ap,:));
  % 差分解を meshgrid に収める
  u(2:Ny,2:Nx)=reshape(U,Ny-1,Nx-1);
  % 境界条件を用いて境界上の格子点での値を求める
  u(:,1)=b_l(Y);
  u(:,Nx+1)=b_r(Y);
  u(1,:)=b_b(X);
  u(Ny+1,:)=b_t(X);
  % 解のグラフ
  figure; meshc(x,y,u); axis([a b c d -2 2]);
  title('graph of the approximate solution');
  drawnow;
  shg;
  disp('何かキーを押して下さい。');
  % 誤差を調べる
  error = supnorm(u-exact_solution(x,y));
  disp(sprintf('Nx=%d, Ny=%d, ||error||=%e, ||u||=%e',...
	       Nx, Ny, error, supnorm(u)));
  % 比較のため厳密解を表示する
  meshc(x,y,exact_solution(x,y)); axis([a b c d -2 2]);
  title('graph of the exact solution'); drawnow; shg
end

入手するには、ターミナルで 
curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2d_nh.m
これを実行するには、例えば 
cp -p poisson2d_nh.m ~/Documents/MATLAB
とコピーして (私はシンボリック・リンクをはることにしている)、 MATLAB の中で 
>> poisson2d_nh
(分割数はデフォールトの $ N_x=50$, $ N_y=50$ が採用される。)
あるいは、分割数 $ N_x$, $ N_y$ を指定して 
>> poisson2d_nh(100,100)
とする。

図 3: $ Nx=Ny=200$ 厳密解と差分解 (違いが見て分かるはずはないが…)
Image poisson2d_200_fig


誤差は次のようになった。 $ O(h_x^2+h_y^2)$ となっていると考えている。 
>> poisson2d_nh
Nx=50, Ny=50, hx=0.020000, hy=0.020000, betax=2500.000000, betay=2500.000000
Nx=50, Ny=50, ||error||=5.924640e-05, ||u||=1.000000e+00

>> poisson2d_nh(100,100)
Nx=100, Ny=100, hx=0.010000, hy=0.010000, betax=10000.000000, betay=10000.000000

Nx=100, Ny=100, ||error||=1.482114e-05, ||u||=1.000000e+00

>> poisson2d_nh(200,200)
Nx=200, Ny=200, hx=0.005000, hy=0.005000, betax=40000.000000, betay=40000.000000

Nx=200, Ny=200, ||error||=3.705884e-06, ||u||=1.000000e+00

>> poisson2d_nh(400,400)
Nx=400, Ny=400, hx=0.002500, hy=0.002500, betax=160000.000000, betay=160000.0000
00
Nx=400, Ny=400, ||error||=9.265347e-07, ||u||=1.000000e+00


桂田 祐史