% poisson2n_nh.m
% 長方形領域におけるPoisson -△u=f (非同次Neumann境界条件) を差分法で解く
% 入手 curl -O https://m-katsurada.sakura.ne.jp/program/fdm/poisson2n_nh.m
% 一応の完成をして、poisson2n_nh.m として公開した (2024/8/27)。

function poisson2n_nh(Nx,Ny)
  arguments
    Nx=50
    Ny=50
  end
  % Neumann境界条件の場合の係数行列
  % poisson2n_mat.m
  function A = poisson2n_mat(W,H,Nx,Ny)
    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;
    hx = W / Nx;
    hy = H / Ny;
    betax = 1.0 / hx^2;
    betay = 1.0 / hy^2;
    Kx=betax * (2*Ix - Jx);
    Ky=betay * (2*Iy - Jy);
    % row major (y changes first, l=j+(Ny+1)*i)
    A=kron(Kx,Iy) + kron(Ix,Ky);
    % column major (x chandes first, l=i+(Nx+1)*j)
    %  A=kron(Iy,Kx) + kron(Ky,Ix);
  end
  % 長方形領域 (a,b)×(c,d)
  a=0; b=1; c=0; d=1;
  %
  function u = exact_solution(x, y)
    u = sin(pi * x) .* sin(pi * y);
  end
  % Poisson 方程式の右辺 (exact_solution の -△)
  function val = f(x, y)
    val = 2 * pi^2 * sin(pi * x) .* sin(pi * y);
  end
  % 長方形の下の辺での b() = -∂u/∂y=-π sin πx cos π0=-π sin πx
  function val = b_b(x)
    val = - pi * sin(pi * x);
  end
  % 長方形の上の辺での b() = ∂u/∂y=π sin πx cos π・1=-π sin πx
  function val = b_t(x)
    val = - pi * sin(pi * x);
  end
  % 長方形の左の辺での b() = -∂u/∂x=-π cos π0 sin πy=-π sin πy
  function val = b_l(y)
    val = - pi * sin(pi * y);
  end
  % 長方形の右の辺での b() = ∂u/∂x=π cosπ・1 sin πy=-π sin πy
  function val = b_r(y)
    val = - pi * sin(pi * y);
  end
  % 誤差を測るためのノルム (最大値ノルム)
  function val = supnorm(a)
    [m,n]=size(a);
    val = norm(reshape(a,m*n,1),Inf);
  end
  % 格子への分割
  %Nx=100;
  %Ny=100;
  hx=(b-a)/Nx;
  hy=(d-c)/Ny;
  betax=1.0/hx^2;
  betay=1.0/hy^2;
  disp(sprintf('Nx=%d, Ny=%d, hx=%f, hy=%f, betax=%f, betay=%f', ...
	       Nx, Ny, hx, hy, betax, betay));

  % 差分方程式 A U^{n+1}=B U^n の行列
  A=poisson2n_mat(b-a,d-c,Nx,Ny);
  % 対称にするためのスケーリング用の行列
  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(er,inf)>=1
    disp('Aが対称でない。')
    pause
  else
    %disp('Aは対称である。');
  end
  clear er;
  N=(Nx+1)*(Ny+1);
  if N <= 50
    disp('A='); disp(full(A));
  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);

  % Poisson方程式の右辺の関数値 (連立1次方程式に必要)
  u=f(x,y);
  close all hidden
  meshc(x,y,u); axis([a b c d -50 50]); title('graph of f'); drawnow; shg
  % 計算に用いるため、f をベクトル U とする。
  U=reshape(u,N,1);
  if N <= 50
    disp('f='); disp(u);
  end
  % 境界値 (左の辺、右の辺)
  blvector=b_l(Y);
  brvector=b_r(Y);
  lindex=1:Ny+1;
  rindex=Nx*(Ny+1)+1:(Nx+1)*(Ny+1);
  if Ny+1 <= 20
    disp('blvector='); disp(blvector);
    disp('brvector='); disp(brvector);
    disp('lindex=');   disp(lindex);
    disp('rindex=');   disp(rindex);
  end
  % 境界値 (下の辺、上の辺)
  bbvector=b_b(X);
  btvector=b_t(X);
  bindex=1:Ny+1:(Nx+1)*(Ny+1);
  tindex=Ny+1:Ny+1:(Nx+1)*(Ny+1);
  if Nx+1 <= 20
    disp('bbvector='); disp(bbvector);
    disp('btvector='); disp(btvector);
    disp('bindex=');   disp(bindex);
    disp('tindex=');   disp(tindex);
  end

  % 境界条件により分かっている値を右辺に移項
  U(lindex)=U(lindex)+2*hx*betax*blvector;
  U(rindex)=U(rindex)+2*hx*betax*brvector;
  U(bindex)=U(bindex)+2*hy*betay*bbvector;
  U(tindex)=U(tindex)+2*hy*betay*btvector;
  % 係数行列を対称化するためのスケーリングを U にも施す
  U=S*U;
  if N <= 50
    disp('右辺のベクトル='); disp(reshape(U,Ny+1,Nx+1));
  end
  % U_{0,0}=0 として解く
  A(1,1)=1; A(1,2)=0; A(1,Ny+2)=0; A(2,1)=0; A(Ny+2,1)=0;
  U(1)=0.0;
  if N <= 20
    disp('A='); disp(full(A));
  end
  % LU分解
  [AL,AU,ap]=lu(A,'vector');
  if N<= 20
    disp('U='); disp(U);
  end
  % 連立1次方程式を解いて差分解を求める
  U=AU\(AL\U(ap,:));
  % 差分解を meshgrid に収める
  u=reshape(U,Ny+1,Nx+1);
  %
  if N <= 50
    disp('u='); disp(u);
  end
  diff=u(1,1)-exact_solution(1,1);
  %fprintf('差分解(0,0)-厳密解(0,0)=%g\n', diff);
  u=u-diff;
  % 解のグラフ
  minu=min(U);
  maxu=max(U);
  figure; meshc(x,y,u); axis([a b c d -1 1]);
  title('graph of the approximate solution');
  drawnow;
  shg;
  % 誤差を調べる
  error = supnorm(u-exact_solution(x,y));
  disp(sprintf('Nx=%d, Ny=%d, ||error||=%e, ||u||=%e',...
	       Nx, Ny, error, supnorm(u)));
  % 比較のため厳密解を表示する
  meshc(x,y,exact_solution(x,y)); axis([a b c d -1 1]);
  title('graph of the exact solution'); drawnow; shg
end