5.1 Dirichlet境界条件

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

heat2d.m

% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
%   A U^{n+1}=B U^n
% の行列 A, B を求める。 (2015/5/1 作成, 2015/5/30 修正)
% Nx=3; Ny=3; hx=1/Nx; hy=1/Ny; theta=0.5; tau=0.5/(1/hx^2+1/hy^2); lamx=tau/hx^2; lamy=tau/hy^2;
% heat2d(Nx,Ny,lamx,lamy,theta)
% A =
%    1.5000   -0.1250   -0.1250         0
%   -0.1250    1.5000         0   -0.1250
%   -0.1250         0    1.5000   -0.1250
%         0   -0.1250   -0.1250    1.5000
% B =
%    0.5000    0.1250    0.1250         0
%    0.1250    0.5000         0    0.1250
%    0.1250         0    0.5000    0.1250
%         0    0.1250    0.1250    0.5000
function heat2d
  a=0; b=2; c=0; d=1;
  Nx=100;
  Ny=50;
  hx=(b-a)/Nx;
  hy=(d-c)/Ny;
  theta=0.5;
  tau=0.5/(1/hx^2+1/hy^2);
  lambdax=tau/hx^2;
  lambday=tau/hy^2;
  
% 差分方程式 A U^{n+1}=B U^n の行列
  [A,B]=heat2d_mat(Nx,Ny,lambdax,lambday,theta);
  if ((Nx <= 5) && (Ny <= 5))
    A
    B
  end
% 格子点の座標ベクトル 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);
% 初期値 sin(pi x) sin(pi y)
  u=sin(pi*x) .* sin(pi*y);
%
% 初期値のグラフを描く
  disp('初期値')
  mesh(x,y,u);
  [AL,AU,ap]=lu(A,'vector');
  Tmax=1;
  t=tau;
  disp('繰り返し')
  k=0;
  dt=0.005;
  skip=dt/tau;
  U=reshape(u(2:Ny,2:Nx),(Nx-1)*(Ny-1),1);
  while t<=Tmax
    U=B*U;
    U=AU\(AL\U(ap,:));
    if mod(k,skip)==0
      u(2:Ny,2:Nx)=reshape(U,Ny-1,Nx-1);
      meshc(x,y,u);
      axis([a b c d -1 1]);
      drawnow;
    end
    t=t+tau
    k=k+1;
  end

heat2d_mat.m

% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
%   A U^{n+1}=B U^n
% の行列 A, B を求める。 (2015/5/1 作成, 2015/5/31 コメント修正)
% Nx=3; Ny=3; hx=1/Nx; hy=1/Ny; theta=0.5; tau=0.5/(1/hx^2+1/hy^2); lamx=tau/hx^2; lamy=tau/hy^2;
% heat2d_mat(Nx,Ny,lamx,lamy,theta)
% A =
%    1.5000   -0.1250   -0.1250         0
%   -0.1250    1.5000         0   -0.1250
%   -0.1250         0    1.5000   -0.1250
%         0   -0.1250   -0.1250    1.5000
% B =
%    0.5000    0.1250    0.1250         0
%    0.1250    0.5000         0    0.1250
%    0.1250         0    0.5000    0.1250
%         0    0.1250    0.1250    0.5000
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;
% x 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);
% y 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);