% Solves numerically the one-dimensional diffusion equation,
%
%          u_t = u_{xx}      (x,t) \in [0,1]^2
%       u(x,0) = \sin \pi x   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 \pi x \exp(-\pi^2 t),
%
% using forward-time difference, backward-time difference, and 
% Crank-Nicolson methods. 

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

%% FORWARD-TIME DIFFERENCE (EXPLICIT)
%
% U(I,J) = u((I-1)*Dx,(J-1)*Dt), I = 1,...,Nx+1, J = 1,...,Nt+1
% (U(I,J+1) - U(I,J))/Dt = (U(I+1,J) - 2*U(I,J) + U(I-1,J))/Dx^2
% U(I,1) = sin(pi*(I-1)*Dx), U(1,J) = 0, U(Nx+1,J) = 0 
% s = Dt/Dx^2, Nx = 6:
% |U(2,J+1)|   |1-2*s   s     0     0     0  | |U(2,J)|   |s*U(1,J)=0|
% |U(3,J+1)|   |  s   1-2*s   s     0     0  | |U(3,J)|   |     0    |
% |U(4,J+1)| = |  0     s   1-2*s   s     0  | |U(4,J)| + |     0    |
% |U(5,J+1)|   |  0     0     s   1-2*s   s  | |U(5,J)|   |     0    |
% |U(6,J+1)|   |  0     0     0     s   1-2*s| |U(6,J)|   |s*U(7,J)=0|
% Order of local truncation error: Dt + Dx^2
% Conditionally stable (s < 1/2).

clear

Nx = 50; Dx = 1/Nx; x = (0:Nx)'*Dx;
Dt = .49/Nx^2; Nt = fix(1/Dt); t = (0:Nt)'*Dt;
s = Dt/Dx^2;
T = diag(repmat(s,Nx-2,1),1) + (1-2*s)*eye(Nx-1) + diag(repmat(s,Nx-2,1),-1);
U = zeros(Nx-1,Nt+1);
U(:,1) = sin(pi*x(2:end-1));
for J = 1:Nt-1
   U(:,J+1) = T*U(:,J);
end
U = [zeros(1,Nt+1); U; zeros(1,Nt+1)];

figure(1), clf

subplot(221)
plot(x,U(:,[1 fix(rand(1,25)*Nt)+1]),'b')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x)$','FontS',14,'Interp','latex')
text(1.25,1.15,['Forward-Time Difference ($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);
Ux05 = interp2(X,T,U',repmat(.5,size(t)),t);
plot(t,Ux05,'bo'), hold on
h = plot(t,sin(pi*.5)*exp(-pi^2*t),'r');
h =legend(h,'exact'); set(h,'FontSize',14,'Interp','latex')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$u(x=0.5,t)$','FontS',14,'Interp','latex')

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

%% BACKWARD-TIME DIFFERENCE (IMPLICIT)
%
% U(I,J) = u((I-1)*Dx,(J-1)*Dt), I = 1,...,Nx+1, J = 1,...,Nt+1
% (U(I,J) - U(I,J-1))/Dt = (U(I+1,J) - 2*U(I,J) + U(I-1,J))/Dx^2
% U(I,1) = sin(pi*(I-1)*Dx), U(1,J) = 0, U(Nx+1,J) = 0 
% s = Dt/Dx^2, Nx = 6:
% |1+2*s -s    0    0    0  | |U(2,J)|   |U(2,J-1)|   |s*U(1,J)=0|
% | -s  1+2*s -s    0    0  | |U(3,J)|   |U(3,J-1)|   |    0     |
% |  0   -s  1+2*s -s    0  | |U(4,J)| = |U(4,J-1)| + |    0     |
% |  0    0   -s  1+2*s -s  | |U(5,J)|   |U(5,J-1)|   |    0     |
% |  0    0    0   -s  1+2*s| |U(6,J)|   |U(6,J-1)|   |s*U(7,J)=0|
% Order of local truncation error: Dt + Dx^2
% Unconditionally stable.

clear

Nx = 50; Dx = 1/Nx; x = (0:Nx)'*Dx;
Nt = 1500; Dt = 1/Nt; t = (0:Nt)'*Dt;
s = Dt/Dx^2;
T = diag(repmat(-s,Nx-2,1),1) + (1+2*s)*eye(Nx-1) + diag(repmat(-s,Nx-2,1),-1);
U = zeros(Nx-1,Nt+1);
U(:,1) = sin(pi*x(2:end-1));
for J = 2:Nt+1
   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,25)*Nt)+1]),'b')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x)$','FontS',14,'Interp','latex')
text(1.25,1.15,['Backward-Time Difference ($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);
Ux05 = interp2(X,T,U',repmat(.5,size(t)),t);
plot(t,Ux05,'bo'), hold on
h = plot(t,sin(pi*.5)*exp(-pi^2*t),'r');
h =legend(h,'exact'); set(h,'FontSize',14,'Interp','latex')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$u(x=0.5,t)$','FontS',14,'Interp','latex')

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

%% CRANK-NICOLSON (IMPLICIT)
%
% U(I,J) = u((I-1)*Dx,(J-1)*Dt), I = 1,...,Nx+1, J = 1,...,Nt+1
% (U(I,J) - U(I,J-1))/Dt = (U(I+1,J) - 2*U(I,J) + U(I-1,J) + U(I+1,J-1) - 2*U(I,J-1) + U(I-1,J-1))/(2*Dx^2)
% U(I,1) = sin(pi*(I-1)*Dx), U(1,J) = 0, U(Nx+1,J) = 0 
% s = Dt/Dx^2, Nx = 6:
% |2+2*s -s    0    0    0  | |U(2,J)|   |2-2*s  s    0    0    0  | |U(2,J-1)|   |s*(U(1,J-1)+U(1,J))=0|
% | -s  2+2*s -s    0    0  | |U(3,J)|   |  s  2-2*s  s    0    0  | |U(3,J-1)|   |         0           |
% |  0   -s  2+2*s -s    0  | |U(4,J)| = |  0    s  2-2*s -s    0  | |U(4,J-1)| + |         0           |
% |  0    0   -s  2+2*s -s  | |U(5,J)|   |  0    0    s  2-2*s  s  | |U(5,J-1)|   |         0           |
% |  0    0    0   -s  2+2*s| |U(6,J)|   |  0    0    0    s  2-2*s| |U(6,J-1)|   |s*(U(1,J-1)+U(1,J))=0|
% Order of local truncation error: Dt^2 + Dx^2
% Unconditionally stable.

clear

Nx = 50; Dx = 1/Nx; x = (0:Nx)'*Dx;
Nt = 100; Dt = 1/Nt; t = (0:Nt)'*Dt;
s = Dt/Dx^2;
T1 = diag(repmat(-s,Nx-2,1),1) + (2+2*s)*eye(Nx-1) + diag(repmat(-s,Nx-2,1),-1);
T2 = diag(repmat(s,Nx-2,1),1) + (2-2*s)*eye(Nx-1) + diag(repmat(s,Nx-2,1),-1);
U = zeros(Nx-1,Nt+1);
U(:,1) = sin(pi*x(2:end-1));
for J = 2:Nt+1
   U(:,J) = T1\(T2*U(:,J-1));
end
U = [zeros(1,Nt+1); U; zeros(1,Nt+1)];

figure(3), clf
subplot(221)
plot(x,U(:,[1 fix(rand(1,25)*Nt)+1]),'b')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x)$','FontS',14,'Interp','latex')
text(1.25,1.15,['Crank--Nicolson ($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);
Ux05 = interp2(X,T,U',repmat(.5,size(t)),t);
plot(t,Ux05,'bo'), hold on
h = plot(t,sin(pi*.5)*exp(-pi^2*t),'r');
h =legend(h,'exact'); set(h,'FontSize',14,'Interp','latex')
xlabel('$t$','FontS',14,'Interp','latex')
ylabel('$u(x=0.5,t)$','FontS',14,'Interp','latex')
subplot(224)
Ut05 = interp2(X,T,U',x,repmat(.5,size(x)));
plot(x,Ut05,'bo',x,sin(pi*x)*exp(-pi^2*.5),'r')
xlabel('$x$','FontS',14,'Interp','latex')
ylabel('$u(x,t=0.5)$','FontS',14,'Interp','latex'