% -------------------------------------------------------------------------
% 
% The Mathematics behind biological invasions:
%  
%  - University of Utah,
%    Department of Mathematics
%    June 2 - 13, 2003
%  
%  - Fred Adler
%    Mark Lewis
%    Mike Neubert
%    Nancy Sundell-Turner
% 
% Forward Euler solution to a Reaction-Diffusion equation.
% 
% -------------------------------------------------------------------------

%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Pre - Initialization:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clf;
F_P = [370 240 900 700];
Fig = gcf;
set(Fig, 'position', F_P);



%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Initialization:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Equation Constants:
K   = 2.0;                    % Diffusion constant


%%%%% Numerical Constant:
Nu  = 0.4;                    % Grid constant 
                              % ( < 0.5 for numerical stability of Forward-Euler.)

%%%%% Spatial Discretization:
Xo  =  0;                     % Initial X-coordinate
Xf  =  50.0;                  % Final X-coordinate
Dx  =  0.1;                   % Space step
X   = [Xo : Dx : Xf];         % Vector of X indices
Nx  = size(X,2);              % Number of grid points


%%%%% Temporal Discretization:
Tf  = 10;                      % Time of simulation
Dt  = Nu * (Dx)^2/K;          % Set the time step for numerical stability


%%%%% Plot Discretization:
Plot_Count = 10;
Plot_Ndx   = 0;
Plot_Base  = floor((Tf/Dt)/Plot_Count);


%%%%% Wave Front:
% Population threshold for the wave speed calculation.
wave_front = 0.01;



%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Difference Matrix Initialization:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Initialize the Forward-Euler difference Matrix:
Iv(1:Nx) = 1;                 % Vector of 1s
A        = (1 - 2*Nu)*eye(Nx) + ...
           Nu*diag(Iv(1:Nx-1),1) + ...
           Nu*diag(Iv(1:Nx-1),-1);


%%%%% Modify the difference matrix for the boundary conditions:
% Neumann boundary conditions:
A(1,1)     =  0;
A(1,2)     =  0;
A(Nx,Nx)   =  0;
A(Nx,Nx-1) =  0;


%%%%% Make the difference matrix sparse for computational efficiency:
A = sparse(A);



%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Initial Conditions:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Compact initial conditions at the origin:
U       = zeros( Nx,1);
U(1:50)  = 1.0;
ndx(1)  = 0;                       % Plot Variable


%%%%% Plot the initial population:
subplot(2,1,1);
plot( X, U, 'k-')
xlabel( 'x');
ylabel( 'u(x,t)');
title(  'Population Density');



%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 
%%%%% Solve the equation:
%%%%% 
%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Compute the solution for T=0 to T=Tf:
for k = 1 : floor(Tf/Dt)
  
% Compute the reaction term:
% Logistic Growth:
  Ur   = Dt*(U - U.*U);

% Exponential Growth:
%  Ur   = Dt*U;
  
% Compute the diffusion term and add the reaction term:
  U = A*U + Ur;
  
% Modify the solution for the boundary conditions:
% Neumann conditions
  U(1) = U(2);
  U(Nx) = U(Nx-1);
  
% Plot the solution and compute the wave front location:
  if mod(k,Plot_Base) == 0
    Plot_Ndx = Plot_Ndx + 1
    
%   Plot the density:
    subplot(2,1,1);
    hold on
    plot( X, U, 'k-')
    hold off
    pause( 0.5);
    
%   Compute the location of the wave front:
    T_1 = U - wave_front;
    T_2 = sign( T_1);
    [Y,ndx(Plot_Ndx+1)]  = min(T_2);
    
%   Plot the location of the wave front:
    subplot(2,1,2);
    plot( [0:Plot_Base:Plot_Ndx*Plot_Base]*Dt, ndx*Dx, 'ko-');
    axis( [0, Tf, 0, 1.1*max(ndx*Dx)] );
    xlabel( 'Time');
    ylabel( 'X');
    title(  'Wave front location');
    
  end
  
end