# include <math.h>
# include <stdlib.h>
# include <stdio.h>

# include "pce_ode_hermite.h"

/******************************************************************************/

double he_double_product_integral ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    he_double_product_integral(): integral of He(i,x)*He(j,x)*e^(-x^2/2).

  Discussion:

    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x) exp(-x^2/2) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 March 2012

  Author:

    John Burkardt

  Reference:

    Dongbin Xiu,
    Numerical Methods for Stochastic Computations: A Spectral Method Approach,
    Princeton, 2010,
    ISBN13: 978-0-691-14212-8,
    LC: QA274.23.X58.

  Input:

    int I, J, the polynomial indices.

  Output:

    double HE_DOUBLE_PRODUCT_INTEGRAL, the value of the integral.
*/
{
  double value;

  if ( i == j )
  {
    value = r8_factorial ( i );
  }
  else
  {
    value = 0.0;
  }
  return value;
}
/******************************************************************************/

double he_triple_product_integral ( int i, int j, int k )

/******************************************************************************/
/*
  Purpose:

    he_triple_product_integral(): integral of He(i,x)*He(j,x)*He(k,x)*e^(-x^2/2).

  Discussion:

    VALUE = integral ( -oo < x < +oo ) He(i,x)*He(j,x)*He(k,x) exp(-x^2/2) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    18 March 2012

  Author:

    John Burkardt

  Reference:

    Dongbin Xiu,
    Numerical Methods for Stochastic Computations: A Spectral Method Approach,
    Princeton, 2010,
    ISBN13: 978-0-691-14212-8,
    LC: QA274.23.X58.

  Input:

    int I, J, K, the polynomial indices.

  Output:

    double HE_TRIPLE_PRODUCT_INTEGRAL, the value of the integral.
*/
{
  int s;
  double value;

  s = ( i + j + k ) / 2;

  if ( s < i || s < j || s < k )
  {
    value = 0.0;
  }
  else if ( ( ( i + j + k ) % 2 ) != 0 )
  {
    value = 0.0;
  }
  else
  {
    value = r8_factorial ( i ) / r8_factorial ( s - i ) 
          * r8_factorial ( j ) / r8_factorial ( s - j ) 
          * r8_factorial ( k ) / r8_factorial ( s - k );
  }

  return value;
}
/******************************************************************************/

void pce_ode_hermite ( double ti, double tf, int nt, double ui, int np, 
  double alpha_mu, double alpha_sigma, double t[], double u[] )

/******************************************************************************/
/*
  Purpose:

    pce_ode_hermite() applies the polynomial chaos expansion to a scalar ODE.

  Discussion:

    The deterministic equation is

      du/dt = - alpha * u,
      u(0) = u0

    In the stochastic version, it is assumed that the decay coefficient
    ALPHA is a Gaussian random variable with mean value ALPHA_MU and variance
    ALPHA_SIGMA^2.

    The exact expected value of the stochastic equation will be

      u(t) = u0 * exp ( t^2/2)

    This should be matched by the first component of the polynomial chaos
    expansion.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 September 2012

  Author:

    John Burkardt

  Input:

    double TI, TF, the initial and final times.

    int  NT, the number of output points.

    double UI, the initial condition.

    int NP, the degree of the expansion.  Polynomials 
    of degree 0 through NP will be used.

    double ALPHA_MU, ALPHA_SIGMA, the mean and standard 
    deviation of the decay coefficient.

  Output:

    double T[NT+1], U[(NT+1)*(NP+1)], the times and the PCE 
    coefficients at the successive time steps.
*/
{
  double dp;
  double dt;
  int i;
  int it;
  int j;
  int k;
  double t1;
  double t2;
  double term;
  double tp;
  double *u1;
  double *u2;

  u1 = ( double * ) malloc ( ( np + 1 ) * sizeof ( double ) );
  u2 = ( double * ) malloc ( ( np + 1 ) * sizeof ( double ) );

  dt = ( tf - ti ) / ( double ) ( nt );
/*
  Set the PCE coefficients for the initial time.
*/
  t1 = ti;

  u1[0] = ui;
  for ( j = 1; j <= np; j++ )
  {
    u1[j] = 0.0;
  }
/*
  Copy into the output arrays.
*/
  t[0] = t1;
  for ( j = 0; j <= np; j++ )
  {
    u[0+j*(nt+1)] = u1[j];
  }
/*
  Time integration.
*/
  for ( it = 1; it <= nt; it++ )
  {
    t2 = ( ( double ) ( nt - it ) * ti   
         + ( double ) (      it ) * tf ) 
         / ( double ) ( nt      );

    for ( k = 0; k <= np; k++ )
    {
      dp = he_double_product_integral ( k, k );

      term = - alpha_mu * u1[k];

      i = 1;
      for ( j = 0; j <= np; j++)
      {
        tp = he_triple_product_integral ( i, j, k );
        term = term - alpha_sigma * u1[j] * tp / dp;
      }
      u2[k] = u1[k] + dt * term;
    }
/*
  Prepare for next step.
*/
    t1 = t2;
    for ( j = 0; j <= np; j++ )
    {
      u1[j] = u2[j];
    }
/*
  Copy into the output arrays.
*/
    t[it] = t1;
    for ( j = 0; j <= np; j++ )
    {
      u[it+j*(nt+1)] = u1[j];
    }
  }

  free ( u1 );
  free ( u2 );

  return;
}
/******************************************************************************/

double r8_factorial ( int n )

/******************************************************************************/
/*
  Purpose:

    r8_factorial() computes the factorial of N.

  Discussion:

    factorial ( N ) = product ( 1 <= I <= N ) I

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 January 1999

  Author:

    John Burkardt

  Input:

    int N, the argument of the factorial function.
    If N is less than 1, the function value is returned as 1.

  Output:

    double R8_FACTORIAL, the factorial of N.
*/
{
  int i;
  double value;

  value = 1.0;

  for ( i = 1; i <= n; i++ )
  {
    value = value * ( double ) ( i );
  }

  return value;
}