function cppspm2d
% cppspm2d: Chebyshev pseudospectral Poisson solver with penalty method in a 2D rectangular domain.
% Citation: Tzyy-Leng Horng and Chun-Hao Teng, 2012, "An error minimized pseudospectral penalty direct Poisson solver," 
% Journal of Computational Physics, 231(6):2498¡V2509 
% uxx+uyy=f(x,y); a test case with exact solution ue(x,y)=sin(4*pi*X).*sin(4*pi*Y);
% BC: amx*u(-1,y)-bmx*ux(-1,y)=gmx(y), at x=-1,
%       apx*u(1,y)+bpx*ux(1,y)=gpx(y), at x=1,
%       amy*u(x,-1)-bmy*uy(x,-1)=gmy(x), at y=-1,
%       apy*u(x,1)+bpy*uy(x,1)=gpy(x), at y=1.
% PS: This code is adapted from poisson_penalty_2d_benchmarking_final.m
close all;
  icho=input('Input 1 for less accuracte collocation derivative matrix, else for more accurate, default is 1: ');
  if isempty(icho)
      icho=1;
  end
  amx=input('Input amx, default is 1: ');
  if isempty(amx)
      amx=1;
  end
  bmx=input('Input bmx, default is 0: ');
  if isempty(bmx)
      bmx=0;
  end
  apx=input('Input apx, default is 1: ');
  if isempty(apx)
      apx=1;
  end
  bpx=input('Input bpx, default is 0: ');
  if isempty(bpx)
      bpx=0;
  end
  amy=input('Input amy, default is 1: ');
  if isempty(amy)
      amy=1;
  end
  bmy=input('Input bmy, default is 0: ');
  if isempty(bmy)
      bmy=0;
  end
  apy=input('Input apy, default is 1: ');
  if isempty(apy)
      apy=1;
  end
  bpy=input('Input bpy, default is 0: ');
  if isempty(bpy)
      bpy=0;
  end
  Nvec=input('Input N in vector, default is [16 20 24 32 40 48 64]: ');
  if isempty(Nvec)
      Nvec=[16 20 24 32 40 48 64];
  end
  errl2vec=[]; errinfvec=[];
  for N=Nvec
  Nx=N; Ny=N;
  cwx=ones(Nx+1,1); cwx([1 end])=2; 
  cwy=ones(Ny+1,1); cwy([1 end])=2; 
  cwmat=cwy*cwx.'; 
  % computing errorminized tau's in x direction:
  am=amx; bm=bmx; ap=apx; bp=bpx;
  gp=(ap/(N^2-1)-bp)/(2*N^2)+(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  gm=(ap/(N^2-1)-bp)/(2*N^2)-(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  hp=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)+(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  hm=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)-(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  taumoptx=(-1)^N*(ap^2+2*(ap+bp)^2)/((2*(ap+bp)*(am+bm)-am*ap)*gm+(ap^2+2*(ap+bp)^2)*hm);
  taupoptx=-(am^2+2*(am+bm)^2)/((am^2+2*(am+bm)^2)*gp+(2*(am+bm)*(ap+bp)-am*ap)*hp);
  % computing errorminized tau's in y direction:
  am=amy; bm=bmy; ap=apy; bp=bpy;
  gp=(ap/(N^2-1)-bp)/(2*N^2)+(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  gm=(ap/(N^2-1)-bp)/(2*N^2)-(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  hp=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)+(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  hm=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)-(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  taumopty=(-1)^N*(ap^2+2*(ap+bp)^2)/((2*(ap+bp)*(am+bm)-am*ap)*gm+(ap^2+2*(ap+bp)^2)*hm);
  taupopty=-(am^2+2*(am+bm)^2)/((am^2+2*(am+bm)^2)*gp+(2*(am+bm)*(ap+bp)-am*ap)*hp);
  taumx=taumoptx; taupx=taupoptx; taumy=taumopty; taupy=taupopty;
  if icho==1
  [Dx,x] = cheb(Nx); Dx=flipud(fliplr(Dx)); x=flipud(x); Dxx=Dx^2; eyex=eye(Nx+1);
  else
  [x,tmp] = chebdif(Nx+1,2); Dx=tmp(:,:,1); Dxx=tmp(:,:,2); eyex=eye(Nx+1);
  Dx=reshape(Dx,Nx+1,Nx+1); Dxx=reshape(Dxx,Nx+1,Nx+1); 
  Dx=flipud(fliplr(Dx)); Dxx=flipud(fliplr(Dxx)); x=flipud(x);
  end
  Dxx(1,:)=Dxx(1,:)-taumx*(amx*eyex(1,:)-bmx*Dx(1,:)); 
  Dxx(end,:)=Dxx(end,:)-taupx*(apx*eyex(end,:)+bpx*Dx(end,:));
  [Xeig,Sigma]=eig(Dxx.'); sigma=diag(Sigma); Xeiginv=inv(Xeig);
  if icho==1
  [Dy,y] = cheb(Ny); Dy=flipud(fliplr(Dy)); y=flipud(y); Dyy=Dy^2; eyey=eye(Ny+1);
  else
  [y,tmp] = chebdif(Ny+1,2); Dy=tmp(:,:,1); Dyy=tmp(:,:,2); eyey=eye(Ny+1);
  Dy=reshape(Dy,Ny+1,Ny+1); Dyy=reshape(Dyy,Ny+1,Ny+1); 
  Dy=flipud(fliplr(Dy)); Dyy=flipud(fliplr(Dyy)); y=flipud(y);
  end
  Dyy(1,:)=Dyy(1,:)-taumy*(amy*eyey(1,:)-bmy*Dy(1,:)); 
  Dyy(end,:)=Dyy(end,:)-taupy*(apy*eyey(end,:)+bpy*Dy(end,:));
  [Yeig,Lambda]=eig(Dyy); lambda=diag(Lambda); Yeiginv=inv(Yeig);
  coef=ones(Ny+1,1)*sigma.'+lambda*ones(1,Nx+1);
  [X,Y]=meshgrid(x,y); 
  ue=sin(4*pi*X).*sin(4*pi*Y);
  uex=4*pi*cos(4*pi*X).*sin(4*pi*Y);
  uey=4*pi*sin(4*pi*X).*cos(4*pi*Y);
  uexx=-16*pi^2*sin(4*pi*X).*sin(4*pi*Y);
  ueyy=-16*pi^2*sin(4*pi*X).*sin(4*pi*Y);
  lapu=uexx+ueyy; 
  lapu(:,1)=lapu(:,1)-taumx*(amx*ue(:,1)-bmx*uex(:,1)); 
  lapu(:,end)=lapu(:,end)-taupx*(apx*ue(:,end)+bpx*uex(:,end));
  lapu(1,:)=lapu(1,:)-taumy*(amy*ue(1,:)-bmy*uey(1,:)); 
  lapu(end,:)=lapu(end,:)-taupy*(apy*ue(end,:)+bpy*uey(end,:));
  bb=(Yeiginv*lapu*Xeig)./coef;
  u=Yeig*bb*Xeiginv;
  err=abs(u-ue); 
%  err([1 end],:)=0; err(:,[1 end])=0;
  maxerr=max(err(:));
  errsq=err.^2; 
  errsqocw=errsq./cwmat;
  errl2=sqrt(sum(errsqocw(:))*pi^2/(Nx*Ny));
  errl2vec=[errl2vec; errl2]; errinfvec=[errinfvec; maxerr];
  end
  format short e
  disp('L2-error   L-infinity-error'); 
  [errl2vec errinfvec]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % CHEB  compute D = differentiation matrix, x = Chebyshev grid
  function [D,x] = cheb(N)
  if N==0, D=0; x=1; return, end
  x = cos(pi*(0:N)/N)'; 
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';                  
  D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
  D  = D - diag(sum(D'));                 % diagonal entries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x, DM] = chebdif(N, M)

%  The function [x, DM] =  chebdif(N,M) computes the differentiation 
%  matrices D1, D2, ..., DM on Chebyshev nodes. 
% 
%  Input:
%  N:        Size of differentiation matrix.        
%  M:        Number of derivatives required (integer).
%  Note:     0 < M <= N-1.
%
%  Output:
%  DM:       DM(1:N,1:N,ell) contains ell-th derivative matrix, ell=1..M.
%
%  The code implements two strategies for enhanced 
%  accuracy suggested by W. Don and S. Solomonoff in 
%  SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).
%  The two strategies are (a) the use of trigonometric 
%  identities to avoid the computation of differences 
%  x(k)-x(j) and (b) the use of the "flipping trick"
%  which is necessary since sin t can be computed to high
%  relative precision when t is small whereas sin (pi-t) cannot.
%  Note added May 2003:  It may, in fact, be slightly better not to
%  implement the strategies (a) and (b).   Please consult the following
%  paper for details:   "Spectral Differencing with a Twist", by
%  R. Baltensperger and M.R. Trummer, to appear in SIAM J. Sci. Comp. 

%  J.A.C. Weideman, S.C. Reddy 1998.  Help notes modified by 
%  JACW, May 2003.

     I = eye(N);                          % Identity matrix.     
     L = logical(I);                      % Logical identity matrix.

    n1 = floor(N/2); n2  = ceil(N/2);     % Indices used for flipping trick.

     k = [0:N-1]';                        % Compute theta vector.
    th = k*pi/(N-1);

     x = sin(pi*[N-1:-2:1-N]'/(2*(N-1))); % Compute Chebyshev points.

     T = repmat(th/2,1,N);                
    DX = 2*sin(T'+T).*sin(T'-T);          % Trigonometric identity. 
    DX = [DX(1:n1,:); -flipud(fliplr(DX(1:n2,:)))];   % Flipping trick. 
 DX(L) = ones(N,1);                       % Put 1's on the main diagonal of DX.

     C = toeplitz((-1).^k);               % C is the matrix with 
C(1,:) = C(1,:)*2; C(N,:) = C(N,:)*2;     % entries c(k)/c(j)
C(:,1) = C(:,1)/2; C(:,N) = C(:,N)/2;

     Z = 1./DX;                           % Z contains entries 1/(x(k)-x(j))  
  Z(L) = zeros(N,1);                      % with zeros on the diagonal.

     D = eye(N);                          % D contains diff. matrices.
                                          
for ell = 1:M
          D = ell*Z.*(C.*repmat(diag(D),1,N) - D); % Off-diagonals
       D(L) = -sum(D');                            % Correct main diagonal of D
DM(:,:,ell) = D;                                   % Store current D in DM
end