function advection_ftcs_pde ( )

%*****************************************************************************80
%
%% advection_ftcs_pde solves the advection PDE using the FTCS method.
%
%  Discussion:
%
%    The FTCS method (Forward Time, Center Space) is unstable for the 
%    advection problem.
%
%    Given a smooth initial condition, successive FTCS approximations will
%    exhibit erroneous oscillations of increasing magnitude.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent advection-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'advection_ftcs_pde:\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the constant-velocity advection PDE in 1D,\n' );
  fprintf ( 1, '    du/dt + c du/dx = 0\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with periodic boundary conditions:\n' );
  fprintf ( 1, '    u(0) = u(1)\n' );
  fprintf ( 1, '  and advection velocity\n' );
  fprintf ( 1, '    c = constant\n' );
  fprintf ( 1, '  and initial condition\n' );
  fprintf ( 1, '    u(0,x) = (10x-6)^2 (8-10x)^2 for 0.6 <= x <= 0.8\n' );
  fprintf ( 1, '           = 0 elsewhere.\n' );
  fprintf ( 1, '  and NX equally spaced nodes in X,\n' );
  fprintf ( 1, '  and NT equally spaced points in T,\n' );
  fprintf ( 1, '  using the FTCS method:\n' );
  fprintf ( 1, '    FT: Forward Time  : du/dt = (u(t+dt,x)-u(t,x))/dt\n' );
  fprintf ( 1, '    CS: Centered Space: du/dx = (u(t,x+dx)-u(t,x-dx))/2/dx\n' );

  nx = 101;
  xmin = 0.0;
  xmax = 1.0;
  dx = ( xmax - xmin ) / ( nx - 1 );
  x = linspace ( xmin, xmax, nx );

  nt = 201;
  tmin = 0.0;
  tmax = 1.0;
  dt = ( tmax - tmin ) / ( nt - 1 );
  t = linspace ( tmin, tmax, nt );

  h = zeros ( nt, 1 );

  c = 0.5;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  Space ranges from %g to %g\n', xmin, xmax );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
  fprintf ( 1, '  Time ranges from %g to %g\n', tmin, tmax );
  fprintf ( 1, '  Constant velocity C = %g\n', c );
%
%  Setting up im1, i, ip1 allows us to define the centered spatial difference
%  with the periodic boundary condition.
%
  im1 = [ nx, 1:nx-2, nx-1 ];
  i   = [ 1,  2:nx-1, nx   ];
  ip1 = [ 2,  3:nx,   1    ];
%
%  Initialize the movie.
%
  colormap ( 'jet' );
  moviename = 'advection_ftcs_pde.avi';
  Video = VideoWriter ( moviename, 'Motion JPEG AVI' );
  Video.FrameRate = 10;
  open ( Video );
%
%  Compute the solution at each time.
%
  for j = 1 : nt
%
%  Apply the initial condition to get solution at initial time.
%
    if ( j == 1 )
      u = advection_pde_initial ( x );
      umax = max ( u ) * ones ( 1, nx );
      umin = min ( u ) * ones ( 1, nx );
%
%  At later times, rearrange
%
%    u(x(i),t(j)) - u(x(i),t(j-1))       u(x(i+1),t(j-1)) - u(x(i-1),t(j-1))
%    ----------------------------- + c * ----------------------------------- = 0
%                dt                                    2 dx
%
%  to get
%
%    u(x(i),t(j)) = u(x(i),t(j-1)) - c * dt / ( 2 dx ) * ( u(x(i+1),t(j-1)) - u(x(i-1),t(j-1)) )
%
    else
      u(i) = u(i) - c * dt / dx / 2.0 * ( u(ip1) - u(im1) );
    end

    h(j) = advection_pde_conserved ( x, u );
%
%  Display the solution.
%
    plot ( x, umin, 'r-', ...
           x, umax, 'r-', ...
           x, u, 'b-', ...
           'linewidth', 3 );

    grid ( 'on' );
    xlabel ( '<--- x --->' );
    ylabel ( '<--- u(x,t) --->' );
    axis ( [ 0.0, 1.0, -2.0, +2.0 ] )
    title_string = sprintf ( 'Advection PDE with FTCS, Time = %6.2f\n', t(j) );
    title ( title_string );

    if ( j == 1 )
      filename = 'advection_ftcs_pde_initial.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( j == nt )
      filename = 'advection_ftcs_pde_final.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    end

    frame = getframe ( 1 );
    writeVideo ( Video, frame );

    pause ( 0.125 );

  end
%
%  Close the movie.
%
  close ( Video );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );
%
%  Plot the conservation quantity.
%
  plot ( t, h, 'r-', ...
         [ tmin, tmax ], [ h(1), h(1) ], 'b--', ...
         [ tmin, tmax ], [ 0, 0 ], 'k-', ...
          'linewidth', 3 );

  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- h(t) --->' );
  legend ( 'H computed', 'H initial', 'X axis', 'Location', 'South' );
  title ( 'Conservation for advection pde using FTCS' );

  filename = 'advection_ftcs_pde_conserved.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'advection_ftcs_pde\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function h = advection_pde_conserved ( x, u )

%*****************************************************************************80
%
%% advection_pde_conserved evaluates the conserved quantity for advection PDE.
%
%  Discussion:
%
%    The integral of U over the region should be conserved.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(1,NX), the coordinates of the nodes.
%
%    real U(1,NX), the current solution.
%
%  Output:
%
%    real H, the value of the conserved quantity.
%
  nx = size ( x, 2 );
  dx = ( max ( x ) - min ( x ) ) / ( nx - 1 );
  h = sum ( u ) * dx;

  return
end
function u = advection_pde_initial ( x )

%*****************************************************************************80
%
%% advection_pde_initial sets the initial condition for the model advection PDE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(1,NX), the coordinates of the nodes.
%
%  Output:
%
%    U(1,NX), the initial condition for the solution.
%
  u = zeros ( size ( x ) );
  k = find ( 0.6 <= x & x <= 0.8 );
  u(k) = ( 10.0 * x(k) - 6.0 ).^2 .* ( 8.0 - 10.0 * x(k) ).^2;

  return
end
