5.2 Neumann境界条件

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

heat2n.m

% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
%   A U^{n+1}=B U^n
% の行列 A, B を求める。
% 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]=heat2n_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);
  if Nx<=5 && Ny<=5
    x
    y
    X
    Y
    u
  end
%
% 初期値のグラフを描く
  disp('初期値')
  mesh(x,y,u);
  [AL,AU,AP]=lu(A);
  Tmax=1;
  t=tau;
  disp('繰り返し')
  k=0;
  dt=0.005;
  skip=dt/tau;
  v=u;
  v(1,:)=v(1,:)/sqrt(2); v(Ny+1,:)=v(Ny+1,:)/sqrt(2);
  v(:,1)=v(:,1)/sqrt(2); v(:,Nx+1)=v(:,Nx+1)/sqrt(2);
  V=reshape(v,(Nx+1)*(Ny+1),1);
  while t<=Tmax
    V=AU\(AL\(AP*(B*V)));
    if mod(k,skip)==0
      u(1:Ny+1,1:Nx+1)=reshape(V,Ny+1,Nx+1);
      u(1,:)=u(1,:)*sqrt(2); u(Ny+1,:)=u(Ny+1,:)*sqrt(2);
      u(:,1)=u(:,1)*sqrt(2); u(:,Nx+1)=u(:,Nx+1)*sqrt(2);
      meshc(x,y,u);
      axis([a b c d -1 1]);
      drawnow;
    end
    t=t+tau
    k=k+1;
  end

heat2n_mat.m

% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
%   A U^{n+1}=B U^n
% の行列 A, B を求める。
% 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;
% heat2n_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=[sqrt(2); ones(Nx-2,1); sqrt(2)];
  Jx=sparse(diag(vx,1)+diag(vx,-1));
  vy=[sqrt(2); ones(Ny-2,1); sqrt(2)];
  Jy=sparse(diag(vy,1)+diag(vy,-1));
  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);