% wave_01.m
% Adapted from the code written initially by Dr. D. Pelinovsky
% (Assignment 6, Problem 3);
% Solving the wave equation:
%    u_tt = u_xx for t >= 0 and x in [0,1]
%
% with the initial conditions
%    u(x,0) = sin(pi x),   u_t(x,0) = cos(pi x)
% 
% and the boundary conditions
%    u_x(0,t) = 1,     u_x(1,t) = 0
%

close all;
clear all;

% initializing constants;
h = 0.1; x = 0 : h : 1; nX = length(x); % grid in x
tau = 0.05; r = tau/h; t = 0 : tau : 1; nT = length(t); % grid in t
U1 = zeros(nX,nT); % matrix for solution surface
U1(:,1) = sin(pi*x'); % initial value

% Assembling the matrix for the explicit approach
A = 2*(1-r^2)*diag(ones(1,nX))+r^2*(diag(ones(1,nX-1),1)+diag(ones(1,nX-1),-1));  
A(1,2) = 2*r^2; A(nX,nX-1) = 2*r^2; % due to derivative boundary conditions
b = zeros(nX,1); b(1) = -2*h*r^2; % derivative boundary conditions

% solving the system with the explicit approach and plotting the solution;
U1(:,2) = 0.5*A*U1(:,1)+0.5*b+tau*cos(pi*x'); % first step of the explicit method
for k = 3 : nT
        U1(:,k) = A*U1(:,k-1)-U1(:,k-2)+b; % further steps of the explicit method
end
[X,T] = meshgrid(x,t);
figure(1); surf(X,T,U1');
xlabel('x'); ylabel('t'); zlabel('u'); % solution to (a)


U2 = zeros(nX,nT); % matrix for solution surface
U2(:,1) = sin(pi*x'); % initial value
% Assembling the matrix for the implicit approach;
A = (1+r^2)*diag(ones(1,nX))-0.5*r^2*(diag(ones(1,nX-1),1)+diag(ones(1,nX-1),-1)); 
A(1,2) = -r^2; A(nX,nX-1) = -r^2; % due to derivative boundary conditions

% solving the system with the implicit approach and plotting the solution
U2(:,2) = A \ (U2(:,1)+0.5*b+tau*A*cos(pi*x')); % first step of the implicit method
for k = 3 : nT
        U2(:,k) = A \ (2*U2(:,k-1)-A*U2(:,k-2)+b); % further steps of the implicit method
end
figure(2); surf(X,T,U2');
xlabel('x'); ylabel('t'); zlabel('u'); % solution to (b)

% Calculating the distance between the solutions as a function of time;
E = max(abs(U1 - U2));
figure(3); 
plot(t,E,'*',t,E); 
xlabel('t'); ylabel('maximum distance between two solutions'); % solution to (c)


% Repeating solution with the explicit method, but for the value of
% the time step that results in an instability
tau = 0.15; r = tau/h; t = 0 : tau : 0.5; nT = length(t); % grid in t
U1 = zeros(nX,nT); % matrix for solution surface
U1(:,1) = sin(pi*x'); % initial value
A = 2*(1-r^2)*diag(ones(1,nX))+r^2*(diag(ones(1,nX-1),1)+diag(ones(1,nX-1),-1));  
A(1,2) = 2*r^2; A(nX,nX-1) = 2*r^2; % due to derivative boundary conditions
U1(:,2) = 0.5*A*U1(:,1)+0.5*b+tau*cos(pi*x'); % first step of the explicit method
for k = 3 : nT
        U1(:,k) = A*U1(:,k-1)-U1(:,k-2)+b; % further steps of the explicit method
end
[X,T] = meshgrid(x,t);
figure(4); surf(X,T,U1');
xlabel('x'); ylabel('t'); zlabel('u'); % solution to (d)
disp('Explicit method is unstable for r > 1'); 
