% Solves numerically a two-dimensional wave equation, using explicit or
% implicit methods.

% Francisco J. Beron-Vera, 2006/15/05

clear

disp('Solves numerically the two-dimensional wave equation:') 
disp(' ')
disp('     u_{tt} = u_{xx} +u_{yy}          (x,y,t) \in [0,1] \times [0,5],')
disp('   u(x,y,0) = \sin 2\pi x \sin 2\pi y (x,y)   \in [0,1]^2,')
disp(' u_t(x,y,0) = 0                       (x,y)   \in [0,1]^2,')
disp('   u(0,y,t) = 0                       (y,t)   \in [0,1] \times [0,5],')            
disp('   u(1,y,t) = 0                       (y,t)   \in [0,1] \times [0,5],')            
disp('   u(x,0,t) = 0                       (x,t)   \in [0,1] \times [0,5],')             
disp('   u(x,1,t) = 0                       (x,t)   \in [0,1] \times [0,5],') 
disp(' ')
disp('whose exact solution is')
disp(' ')
disp('   u(x,y,t) = \sin 2\pi x \sin 2\pi y \cos(2\pi t).')
disp(' ')
disp('Available methods:') 
disp(' ')
disp('(1) Explicit: U(t+Dt) = T*U(t) - U(t-Dt) (N.B. (Dx^2 + Dy^2)/Dt^2 > 1)')
disp('(2) Implicit: T*U(t+Dt) = 2*U(t) - T*U(t-Dt)')
disp(' ')
method = input('Your method choice: ');
Nx = input('Number of grid points in x direction: ') - 1;
Ny = input('Number of grid points in y direction: ') - 1;
Nt = input('Number of time steps = ');

u = @(x,y,t) sin(2*pi*x).*sin(2*pi*y)*cos(2*pi*t);

Dx = 1/Nx; x = (0:Nx)'*Dx;
Dy = 1/Ny; y = (0:Ny)'*Dy; 
Dt = 5/Nt; t = (0:Nt)'*Dt;

sx = Dt^2/Dx^2;
sy = Dt^2/Dy^2;

nx = Nx-1;  
ny = Ny-1;
Sx = repmat(sx,nx,1);
Sy = repmat(sy,ny,1);
E  = speye(nx*ny);
Ex = speye(nx); 
Ey = speye(ny);

d2dx2  = spdiags([Sx -2*Sx Sx],-1:1,nx,nx);      
d2dy2  = spdiags([Sy -2*Sy Sy],-1:1,ny,ny);
nabla2 = kron(d2dy2,Ex) + kron(Ey,d2dx2);   

if method == 1
   T = 2*E + nabla2/2;
else
   T = E - nabla2/4;
end

% Initializing U...
[x y] = ndgrid(x,y);
UdU = u(x,y,0);

subplot(221)
surf(x,y,UdU), shading flat
axis([0 1 0 1 -1 1])
xlabel('x'), ylabel('y'), zlabel('u'), title('NUMERIC')
text(2.05,1,.5,'t = 0','Hor','cen')

subplot(222)
surf(x,y,UdU), shading flat
axis([0 1 0 1 -1 1])
xlabel('x'), ylabel('y'), zlabel('u'), title('EXACT')

drawnow

Uold = UdU;
Uold([1 end],:) = [];
Uold(:,[1 end]) = []; 
Uold = Uold(:);
U = Uold;

% Time stepping U...
for J = 2:Nt+1
   
   if method == 1
      Unew = T*U - Uold; 
   else
       Unew = T\(2*U - T*Uold); 
   end
   Uold = U; U = Unew;
   
   subplot(221)
	UdU(2:end-1,2:end-1) = reshape(Unew,[nx ny]);
   surf(x,y,UdU), shading flat
   axis([0 1 0 1 -1 1]) 
   xlabel('x'), ylabel('y'), zlabel('u'),title('NUMERIC')
   text(2.05,1,.5,['t = ' num2str(t(J))],'Hor','cen')
   
   subplot(222)
   surf(x,y,u(x,y,t(J))), shading flat
   axis([0 1 0 1 -1 1]) 
   xlabel('x'), ylabel('y'), zlabel('u'), title('EXACT')
   
   drawnow
   
en