% heat2d_mat.m
% 長方形領域における熱方程式 u_t=△u (Dirichlet境界条件) を解くための差分方程式
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/heat2d_mat.m
% 解説 https://m-katsurada.sakura.ne.jp/labo/text/heat2d.pdf
%   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);
  
