% Solves numerically the one-dimensional wave equation,
%
%          u_t = u_{xx}      (x,t) \in [0,1]^2
%       u(x,0) = \sin 2\pi x  x    \in [0,1]
%     u_t(x,0) = 0            x    \in [0,1]
%       u(0,t) = 0              t  \in [0,1]            
%       u(1,t) = 0              t  \in [0,1]             
%
% whose exact solution is
%
%       u(x,t) = \sin 2\pi x \cos 2\pi t),
%
% using explicit and implicit methods. 

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

%% EXPLICIT METHOD
% U(I,J) = u((I-1)*Dx,(J-1)*Dt), I = 1,...,Nx+1, J = 1,...,Nt+1
% (U(I,J+1) - 2*U(I,J) + U(I,J-1)/Dt^2 = (U(I+1,J) - 2*U(I,J) + U(I-1,J))/Dx^2
% U(I,1) = sin(2*pi*(I-1)*Dx), U(I,2) = U(I,2), U(1,J) = 0, U(Nx+1,J) = 0 
% Order of local truncation error: Dt^2 + Dx^2
% Conditionally stable: Dx^2/Dt^2 > 1 (Courant-Friedrichs-Lewy)

clear

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

Nx = 50; Dx = 1/Nx; x = (0:Nx)'*Dx;
Dt = .25/Nx; Nt = fix(1/Dt); t = (0:Nt)'*Dt;
s = Dt^2/Dx^2;

S = repmat(s,Nx-1,1);
T = spdiags([S 2-2*S S],-1:1,Nx-1,Nx-1);

U = zeros(Nx-1,Nt+1);
U(:,1) = u(x(2:end-1),0);
U(:,2) = U(:,1);
for J = 2:Nt
   U(:,J+1) = T*U(:,J) - U(:,J-1);
end
U = [zeros(1,Nt+1); U; zeros(1,Nt+1)];

figure(1), clf

subplot(221)
plot(x,U(:,[1 fix(rand(1,15)*Nt)+1]),'b')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x)$','FontS',14,'Interp','latex')
text(1.25,1.2,['Explicit Method ($N_x = '...
   num2str(Nx) '$, $N_t = ' num2str(Nt) '$)'],...
   'Hor','cen','FontS',14,'Interp','latex')

subplot(222)
surf(t,x,U), shading flat
set(gca,'YDir','reverse')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$x$','FontS',14,'Interp','latex')
zlabel('$u(x,t)$','FontS',14,'Interp','latex')

subplot(223)
[X T] = meshgrid(x,t);
Ux025 = interp2(X,T,U',repmat(.25,size(t)),t);
plot(t,Ux025,'bo'), hold on
h = plot(t,u(.25,t),'r');
h =legend(h,'exact'); set(h,'FontSize',14,'Interp','latex')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$u(x=0.25,t)$','FontS',14,'Interp','latex')

subplot(224)
Ut05 = interp2(X,T,U',x,repmat(.5,size(x)));
plot(x,Ut05,'bo',x,u(x,.5),'r')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x,t=0.5)$','FontS',14,'Interp','latex')

%% IMPLICIT METHOD
% U(I,J) = u((I-1)*Dx,(J-1)*Dt), I = 1,...,Nx+1, J = 1,...,Nt+1
% (U(I,J+1) - 2*U(I,J) + U(I,J-1)/Dt^2 = (U(I+1,J+1) - 2*U(I,J+1) + U(I-1,J+1))/(2*Dx^2) +
%                                        (U(I+1,J-1) - 2*U(I,J-1) + U(I-1,J-1))/(2*Dx^2)
% U(I,1) = sin(2*pi*(I-1)*Dx), U(I,2) = U(I,2), U(1,J) = 0, U(Nx+1,J) = 0 
% Order of local truncation error: Dt^2 + Dx^2
% Unconditionally stable.

clear

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

Nx = 50; Dx = 1/Nx; x = (0:Nx)'*Dx;
Nt = 50; Dt = 1/Nt; t = (0:Nt)'*Dt;
s = Dt^2/Dx^2;

S = repmat(s,Nx-1,1);
T = spdiags([-.5*S  1+S -.5*S],-1:1,Nx-1,Nx-1);

U = zeros(Nx-1,Nt+1);
U(:,1) = u(x(2:end-1),0);
U(:,2) = U(:,1);
for J = 2:Nt
   U(:,J+1) = T\(2*U(:,J)-T*U(:,J-1));
end
U = [zeros(1,Nt+1); U; zeros(1,Nt+1)];

figure(2), clf

subplot(221)
plot(x,U(:,[1 fix(rand(1,15)*Nt)+1]),'b')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x)$','FontS',14,'Interp','latex')
text(1.25,1.2,['Implicit Method ($N_x = '...
   num2str(Nx) '$, $N_t = ' num2str(Nt) '$)'],...
   'Hor','cen','FontS',14,'Interp','latex')

subplot(222)
surf(t,x,U), shading flat
set(gca,'YDir','reverse')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$x$','FontS',14,'Interp','latex')
zlabel('$u(x,t)$','FontS',14,'Interp','latex')

subplot(223)
[X T] = meshgrid(x,t);
Ux025 = interp2(X,T,U',repmat(.25,size(t)),t);
plot(t,Ux025,'bo'), hold on
h = plot(t,u(.25,t),'r');
h =legend(h,'exact'); set(h,'FontSize',14,'Interp','latex')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$u(x=0.25,t)$','FontS',14,'Interp','latex')

subplot(224)
Ut05 = interp2(X,T,U',x,repmat(.5,size(x)));
plot(x,Ut05,'bo',x,u(x,.5),'r')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x,t=0.5)$','FontS',14,'Interp','latex'