% heat2d.m
% 長方形領域における熱方程式 u_t=△u (同次Dirichlet境界条件) を差分法で解く
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2d.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/heat2d.pdf
function heat2d
  a=0; b=2; c=0; d=1;
  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;
  
% 差分方程式 A U^{n+1}=B U^n の行列
  [A,B]=heat2d_mat(Nx,Ny,lambdax,lambday,theta);
% 格子点の座標ベクトル 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);
				%
  msg = sprintf("Nx=%d, Ny=%d, τ=%e, θ=%f, λx=%f, λy=%f",...
		Nx, Ny, tau, theta, lambdax, lambday);
  disp(msg);
% 初期値のグラフを描く
  mesh(x,y,u);
  disp('初期値を書きました。何かキーを押してください。')
  pause
% 径数行列のLU分解
  [AL,AU,ap]=lu(A,'vector');

  Tmax=1;
  disp('繰り返し')
  dt=0.005;
  skip=dt/tau;
  U=reshape(u(2:Ny,2:Nx),(Nx-1)*(Ny-1),1);
  Nmax = ceil(Tmax / tau);
  for n=1:Nmax
    t=n*tau;
    U=B*U;
    U=AU\(AL\U(ap,:));
    if mod(n,skip)==0
      u(2:Ny,2:Nx)=reshape(U,Ny-1,Nx-1);
      msg=sprintf("t=%f", t);
      disp(msg);
      meshc(x,y,u);
      axis([a b c d -1 1]);
      drawnow;
    end
  end
