function pce_burgers ( )

%*****************************************************************************80
%
%% PCE_BURGERS applies the polynomial chaos expansion to the Burgers equation.
%
%  Discussion:
%
%    The time-dependent viscous Burgers equation to be solved is:
%
%      du/dt = - d ( u*(1/2-u)) /dx + nu d2u/dx2
%
%    with boundary conditions
%
%      u(-3.0) = 0.0, u(+3.0) = 1.0.
%
%    The viscosity nu is assumed to be an uncertain quantity with
%    normal distribution of known mean and variance.
%
%    A polynomial chaos expansion is to be used, with Hermite polynomial
%    basis functions h(i,x), 0 <= i <= n.
%
%    Because the first two Hermite polynomials are simply 1 and x, 
%    we have that 
%
%      nu = nu_mean * h(0,x) + nu_variance * h(1,x).
%
%    We replace the time derivative by an explicit Euler approximation,
%    so that the equation now describes the value of U(x,t+dt) in terms
%    of known data at time t.
%
%    Now assume that the solution U(x,t) can be approximated
%    by the truncated expansion:
%
%      U(x,t) = sum ( 0 <= i <= n ) c(i,t) * h(i,x)
%
%    In the equation, we replace U by its expansion, and then multiply
%    successively by each of the basis functions h(*,x) to get a set of
%    n+1 equations that can be used to determine the values of c(i,t+dt).
%
%    This process is repeated until the desired final time is reached.
%
%    At any time, the coefficients c(0,t) contain information definining
%    the expected value of u(x,t) at that time, while the higher order coefficients
%    can be used to deterimine higher moments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2012
%
%  Author:
%
%    The original FORTRAN90 version of this program was written by Gianluca Iaccarino.
%    This MATLAB version is by John Burkardt.
%
%  Local parameters:
%
%    Local, real DT, the timestep.
%
%    Local, real DX, the spacing between grid points.
%
%    Local, integer N, the number of intervals in the spatial domain.
%
%    Local, real NUMEAN, the mean of viscosity.
%
%    Local, real NUVARIANCE, the variance of viscosity.
%
%    Local, integer P, the order of the PC expansion.
%
%    Local, real T, the current time.
%
%    Local, real TF, the final integration time.
%
%    Local, real U1(N+1,P+1), the PCE representation at the current time.
%
%    Local, real U2(N+1,P+1), the PCE representation for the next time.
%
%    Local, real X(N+1,1), the grid points.
%
  p = 5;
  n = 32;
  nt = 2000;
  ti = 0.0;
  tf = 2.0;
  dt = ( tf - ti ) / nt;
  numean = 0.25;
  nuvariance = 0.08;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PCE_BURGERS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polynomial Chaos Expansion\n' );
  fprintf ( 1, '  1D Burgers equation\n' );
  fprintf ( 1, '  Original version by Gianluca Iaccarino\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  PCE order       = %d\n', p );
  fprintf ( 1, '  Number of cells = %d\n', n );
  fprintf ( 1, '  Time step       = %g\n', dt );
  fprintf ( 1, '  Initial time    = %g\n', ti );
  fprintf ( 1, '  Final time      = %g\n', tf );
  fprintf ( 1, '  Viscosity Mean  = %g\n', numean );
  fprintf ( 1, '  Viscosity Var   = %g\n', nuvariance );
  fprintf ( 1, '\n' );
%
%  Define some numerical parameters.
%
  dx = 6.0 / n;
  conv = dt / ( 2.0 * dx );

  visc = zeros ( p + 1, 1 );
  visc(1) = numean * dt / ( dx * dx );
  visc(2) = nuvariance * dt / ( dx * dx );
%
%  Define a uniform grid.
%
  x = ( linspace ( - 3.0, + 3.0, n + 1 ) )';
%
%  Set the initial conditions.
%
  u1 = zeros ( n + 1, p + 1 );
  u1(1:n+1,1) = 0.5 + x(1:n+1,1) / 6.0;

  u2 = zeros ( n + 1, p + 1 );
%
%  Time integration.
%
  t1 = ti;
%
%  Write the current solution.
%
  output_filename = 'burgers.history.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, '----------\n' );
  fprintf ( output_unit, 'T = %g\n', t1 );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %10g', x(i) );
    for k = 1 : p + 1
      fprintf ( output_unit, '  %10g', u1(i,k) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Time integration.
%
  for it = 1 : nt

    t2 = ( ( nt - it ) * ti   ...
         + (      it ) * tf ) ...
         / ( nt      );
%
%  Boundary conditions.
%
    u2(1,1:p+1) = 0.0;
    u2(n+1,1) = 1.0;
    u2(n+1,2:p+1) = 0.0;

    for k = 1 : p + 1

      dp = he_double_product_integral ( k - 1, k - 1 );

      for ix = 2 : n
%
%  Viscous term.
%
        term1 = visc(1) * ( u1(ix+1,k) - 2.0 * u1(ix,k) + u1(ix-1,k) );
        i = 2;
        for j = 1 : p + 1
          tp = he_triple_product_integral ( i - 1, j - 1, k - 1 );
          term1 = term1 + visc(i) * ( u1(ix+1,j) - 2.0 * u1(ix,j) + u1(ix-1,j) ) * tp / dp;
        end
%
%  Convective term.
%
        term2 = - conv * 0.5 * ( u1(ix+1,k) - u1(ix-1,k) );
        for j = 1 : p + 1
          for i = 1 : p + 1
            tp = he_triple_product_integral ( i - 1, j - 1, k - 1 );
            term2 = term2 + ( conv * u1(ix,i) * ( u1(ix+1,j) - u1(ix-1,j) ) * tp ) / dp;
          end
        end

        u2(ix,k) = u1(ix,k) + term1 + term2;

      end

    end

    t1 = t2;
    u1(1:n+1,1:p+1) = u2(1:n+1,1:p+1);
%
%  Print solution every 100 time steps.
%
    if ( mod ( it, 100 ) == 0 )
      fprintf ( output_unit, '----------\n' );
      fprintf ( output_unit, 'T = %g\n', t1 );
      for i = 1 : n + 1
        fprintf ( output_unit, '  %10g', x(i) );
        for k = 1 : p + 1
          fprintf ( output_unit, '  %10g', u1(i,k) );
        end
        fprintf ( output_unit, '\n' );
      end
    end

  end

  fclose ( output_unit );
  fprintf ( 1, '  Time history in "%s".\n', output_filename );
%
%  Compute the mean and variance.
%
  umean(1:n+1,1) = u1(1:n+1,1);

  uvariance = zeros ( n + 1, 1 );

  for i = 1 : n + 1
    for j = 2 : p + 1
      dp = he_double_product_integral ( j - 1, j - 1 );
      uvariance(i) = uvariance(i) + u1(i,j).^2 * dp;
    end
  end
%
%  Save the solution at the final time.
%
  output_filename = 'burgers.moments.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, ' X E[U] Var[U]\n' );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %18.8g  %18.8g  %18.8g\n', x(i), umean(i), uvariance(i) );
  end
  fclose ( output_unit );
  fprintf ( 1, '  Moments in "%s".\n', output_filename );

  output_filename = 'burgers.modes.txt';
  output_unit = fopen ( output_filename, 'wt' );
  fprintf ( output_unit, ' X U_0 ... U_P \n' );
  for i = 1 : n + 1
    fprintf ( output_unit, '  %10g', x(i) );
    for k = 1 : p + 1
      fprintf ( output_unit, '  %10g', u1(i,k) );
    end
    fprintf ( output_unit, '\n' );
  end
  fclose ( output_unit );
  fprintf ( 1, '  Final modes in "%s".\n', output_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PCE_BURGERS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end