5.1 Dirichlet境界条件

「長方形領域における熱方程式に対する差分法 -- MATLAB を使って数値計算」

「2次元熱方程式の非同次 Dirichlet 境界値問題を解く差分法プログラムを作る」 という解説を書いてる。

以下は、そこで説明した非同次Dirichlet境界値問題のプログラムである。

heat2d_nh.m

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

function heat2d_nh
  % Dirichlet境界条件の場合の係数行列
  function [A,B]=heat2d_mat(Nx,Ny,lambdax,lambday,theta)
    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));
    Kx=2*Ix-Jx;
    Ky=2*Iy-Jy;
    % row major (y changes first, l=j+(Ny-1)*i)
    A=kron(Ix,Iy)+theta*lambday*kron(Ix,Ky)+theta*lambdax*kron(Kx,Iy);
    B=kron(Ix,Iy)-(1-theta)*lambday*kron(Ix,Ky)-(1-theta)*lambdax*kron(Kx,Iy);
    % column major (x chandes first, l=i+(Nx-1)*j)
    % A=kron(Iy,Ix)+theta*lambdax*kron(Iy,Kx)+theta*lambday*kron(Ky,Ix);
    % B=kron(Iy,Ix)-(1-theta)*lambdax*kron(Iy,Kx)-(1-theta)*lambday*kron(Ky,Ix);
  end

  % 長方形領域 (a,b)×(c,d)
  a=0; b=1; c=0; d=1;
  % Dirichlet境界値 b() 調和関数を選ぶとこれが平衡解 (定常解)
  function val = dbv(x,y)
    val = (x-0.5).*(x-0.5)-(y-0.5).*(y-0.5);
  end
  % 初期値 (平衡解を境界値が0のもので摂動する)
  function val = iv(x,y)
    val = dbv(x,y) + sin(pi*x).*sin(2*pi*y);
  end
  % 厳密解が分かる
  function val = exact(x,y,t)
    val = dbv(x,y) + exp(- 5 * pi * pi * t) * sin(pi*x).*sin(2*pi*y);
  end
  % 長方形の下の辺での b()
  function val = bbottom(x)
    val = dbv(x,c);
  end
  % 長方形の上の辺での b()
  function val = btop(x)
    val = dbv(x,d);
  end
  % 長方形の左の辺での b()
  function val = bleft(y)
    val = dbv(a, y);
  end
  % 長方形の右の辺での b()
  function val = bright(y)
    val = dbv(b, y);
  end
  % 誤差を測るためのノルム (最大値ノルム)
  function val = mynorm(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;
  theta=0.5;
  tau=0.5/(1/hx^2+1/hy^2); % 陽解法の場合のギリギリの刻み(陰解法なので大きくてもOK)
  lambdax=tau/hx^2;
  lambday=tau/hy^2;
  disp(sprintf('Nx=%d, Ny=%d, hx=%f, hy=%f, tau=%e, λx=%f, λy=%f', ...
	       Nx, Ny, hx, hy, tau, lambdax, lambday));

  % 差分方程式 A U^{n+1}=B U^n の行列
  [A,B]=heat2d_mat(Nx,Ny,lambdax,lambday,theta);
  % 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);
  % 初期値
  u= iv(x,y);
  % 極限 (境界値 b の式が調和ならば)
  ulimit = dbv(x,y);

  % 初期値のグラフを描く
  disp('初期値のグラフ')
  meshc(x,y,u); axis([a b c d -2 2]); drawnow;

  U=reshape(u(2:Ny,2:Nx),(Nx-1)*(Ny-1),1);

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

  % どこまで計算するか
  Tmax=0.5;
  Nmax = ceil(Tmax/tau);
  % グラフを表示する時間の間隔
  dt=0.005;
  skip=dt/tau;
  %
  disp('ループに入ります。何かキーを押してください。')
  pause
  t=0;
  maxerror=0;
  for n=0:Nmax
    % 連立1次方程式の右辺の計算
    U=B*U;
    % 移項
    U(lindex)=U(lindex)+lambdax*blvector;
    U(rindex)=U(rindex)+lambdax*brvector;
    U(bindex)=U(bindex)+lambday*bbvector;
    U(tindex)=U(tindex)+lambday*btvector;
    % 連立1次方程式を解いて差分解を求める
    U=AU\(AL\U(ap,:));
    t=(n+1)*tau;
    if mod(n,skip)==0
      u(2:Ny,2:Nx)=reshape(U,Ny-1,Nx-1);
      meshc(x,y,u); axis([a b c d -2 2]); drawnow;
      error = mynorm(u-exact(x,y,t));
      if error > maxerror maxerror = error; end
      disp(sprintf('t=%f, ||error||=%e, ||u||=%e, ||u∞||=%e',...
		   t, error, mynorm(u), mynorm(ulimit)));
    end
  end
  display(sprintf('||u(・,t)-u(・,∞)||=%e', mynorm(u-ulimit)));
  display(sprintf('Nx=%d,Ny=%d,max error=%e', Nx, Ny, maxerror));
end

curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2d_nh.m
cp -p heat2d_nh.m ~/Documents/MATLAB

>> heat2d_nh

>> heat2d_nh(100,100)