5.2 Neumann境界条件

「Neumann 境界条件下の熱方程式に対する差分法」

「非同次 Neumann 境界条件下の熱方程式に対する差分法 -- MATLAB を使って数値計算」

heat2n_nh.m

% heat2n_nh.m
% 長方形領域における熱方程式 u_t=△u (非同次Neumann境界条件) を差分法で解く
% 2024/8/24, written by mk.
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2n_nh.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/heat2n-nonhomo.pdf

function heat2n_nh
  function [A,B]=heat2n_mat0(Nx,Ny,lambdax,lambday,theta)
    Ix=speye(Nx+1,Nx+1);
    Iy=speye(Ny+1,Ny+1);
    vx=ones(Nx,1);
    Jx=sparse(diag(vx,1)+diag(vx,-1));
    Jx(1,2)=2; Jx(Nx+1,Nx)=2;
    vy=ones(Ny,1);
    Jy=sparse(diag(vy,1)+diag(vy,-1));
    Jy(1,2)=2; Jy(Ny+1,Ny)=2;
    Kx=2*Ix-Jx;
    Ky=2*Iy-Jy;
  % column first
    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);
  % row first
  %  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;
  % Neumann境界値 b() 調和関数を選ぶとこれが平衡解 (定常解)
  function val = harmonic(x,y)
    val = (x-0.5).*(y-0.5);
  end
  % 初期値 (平衡解を境界値が0のもので摂動する)
  function val = iv(x,y)
    val = harmonic(x,y) + cos(pi*x).*cos(pi*y);
  end
  % 厳密解が分かる
  function val = exact(x,y,t)
    val = harmonic(x,y) + exp(- 2 * pi * pi * t) * cos(pi*x).*cos(pi*y);
  end
  % 長方形の下の辺での b_b() … -∂u/∂y(x,0,t)
  function val = b_b(x)
    val = 0.5 - x;
  end
  % 長方形の上の辺での b_b() … ∂u/∂y(x,1,t)
  function val = b_t(x)
    val = x - 0.5;
  end
  % 長方形の左の辺での b_l() … -∂u/∂x(0,y,t)
  function val = b_l(y)
    val = 0.5 - y;
  end
  % 長方形の右の辺での b_r() … ∂u/∂x(1,y,t)
  function val = b_r(y)
    val = y - 0.5;
  end
  % 誤差を測るためのノルム (最大値ノルム)
  function val = mynorm(a)
    [m,n]=size(a);
    val = norm(reshape(a,m*n,1),Inf);
  end
  % 格子への分割
  Nx=50;
  Ny=50;
  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;
  display(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]=heat2n_mat0(Nx,Ny,lambdax,lambday,theta);
  % 対称にするためのスケーリング用の行列
  Sx=sparse(diag([1/2;ones(Nx-1,1);1/2]));
  Sy=sparse(diag([1/2;ones(Ny-1,1);1/2]));
  S=kron(Sy,Sx);
  % Aを対称行列に直す
  A=S*A;
  % 対称かどうかチェックする
  er=A-A';
  if norm(full(er),inf)>=1
    disp('Aが対称でない。')
    pause
  end
  % 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);
  % 極限
  ulimit = harmonic(x,y);

  % 初期値のグラフを描く
  disp('初期値のグラフ')
  meshc(x,y,u);
  axis([a b c d -4 4]);
  drawnow;
  disp('何かキーを押して下さい')
  pause

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

  % どこまで計算するか
  Tmax=0.5;
  Nmax = ceil(Tmax/tau);
  % グラフを表示する時間の間隔
  dt=0.005;
  skip=dt/tau;

  U=reshape(u,(Nx+1)*(Ny+1),1);

  disp('ループに入ります。何かキーを押してください。')
  t=0;
  maxerror=0;
  for n=0:Nmax
    % 注: これから計算するのが時刻tでの値 (最初が t=tau)
    U=B*U;
    % 移項
    U(1:Ny+1)=U(1:Ny+1)+2*hx*lambdax*blvector;
    U(end-(Ny+1)+1:end)=U(end-(Ny+1)+1:end)+2*hx*lambdax*brvector;
    U(bindex)=U(bindex)+2*hy*lambday*bbvector;
    U(tindex)=U(tindex)+2*hy*lambday*btvector;
    U=S*U;
    % 連立一次方程式を解いて差分解を求める
    U=AU\(AL\U(ap,:));
    t=(n+1)*tau;
    % 計算結果の表示
    if mod(n,skip)==0
      u=reshape(U,Ny+1,Nx+1);
      meshc(x,y,u); axis([a b c d -1 1]); drawnow;
      % 誤差の表示
      error = mynorm(u-exact(x,y,t));
      if error > maxerror maxerror = error; end
      display(sprintf("t=%f, ||error||=%e, ||u||=%e, ||u∞||=%e",...
		      t, error, mynorm(u), mynorm(ulimit)));
    end
  end
  display(sprintf("max ||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/heat2n_nh.m
cp -p heat2n_nh.m ~/Documents/MATLAB

>> heat2n_nh

>> heat2n_nh(100,100)