root/lang/cplusplus/boost-supplement/branches/numeric-ode/libs/numeric/ode/test_runge_kutta_rk4.cpp @ 6875

/* $Id$ */

#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>

#include <boost_supplement/numeric/ode/runge_kutta.hpp>

#include <boost/test/floating_point_comparison.hpp>
#include <boost/tuple/tuple.hpp>

typedef boost::numeric::ublas::vector<double> dvector;
typedef boost::numeric::ublas::matrix<double> dmatrix;

namespace bs = boost_supplement;

namespace {

void
nr_rk4(double y[], double dydx[], int n, double x, double h, double yout[],
       void (*derivs)(double, double[], double[]))
{
  int i;
  double xh, hh, h6, *dym, *dyt, *yt;
  dym = new double[n+1];
  dyt = new double[n+1];
  yt  = new double[n+1];
  hh = h*0.5;
  h6 = h/6.0;
  xh = x+hh;
  for (i = 1; i <= n; i++) yt[i] = y[i]+hh*dydx[i];
  (*derivs)(xh, yt, dyt);
  for (i = 1; i <= n; i++) yt[i] = y[i]+hh*dyt[i];
  (*derivs)(xh, yt, dym);
  for (i = 1; i <= n; i++) {
    yt[i] = y[i] + h*dym[i];
    dym[i] += dyt[i];
  }
  (*derivs)(x+h, yt, dyt);
  for (i = 1; i <= n; i++)
    yout[i] = y[i] + h6*(dydx[i]+dyt[i]+2.0*dym[i]);
  delete[] yt;
  delete[] dyt;
  delete[] dym;
}

void
nr_rkdumb(double vstart[], int nvar, double x1, double x2, int nstep,
          double* xx, double* y,
          void (*derivs)(double, double[], double[]))
{
  int i, k;
  double x, h;
  double *v, *vout, *dv;

  v = new double[nvar+1];
  vout = new double[nvar+1];
  dv = new double[nvar+1];
  for (i = 1; i <= nvar; i++) {
    v[i] = vstart[i];
    y[i*(nstep+2) + 1] = v[i];
  }
  xx[1] = x1;
  x = x1;
  h = (x2 - x1)/nstep;
  for (k = 1; k <= nstep; k++) {
    (*derivs)(x, v, dv);
    nr_rk4(v, dv, nvar, x, h, vout, derivs);
    if ((double)(x+h) == x)
      throw "Step size too small in routine rkdumb";
    x += h;
    xx[k+1] = x;
    for (i = 1; i <= nvar; i++) {
      v[i] = vout[i];
      y[i*(nstep+2) + k+1] = v[i];
    }
  }

  delete[] dv;
  delete[] vout;
  delete[] v;
}

void
nr_derivative(double x, double y[], double dy[])
{
  dy[1] = y[2];
  dy[2] = y[3];
  dy[3] = y[1];
}

}

dvector
test_derivative(double x, dvector y)
{
  BOOST_ASSERT(3 == y.size());

  dvector dy(3);
  dy[0] = y[1];
  dy[1] = y[2];
  dy[2] = y[0];

  return dy;
}

BOOST_AUTO_TEST_CASE(test_rk4)
{
  dvector y0(3);
  y0[0] = y0[1] = y0[2] = 1.0;
  dvector dydx0(test_derivative(0.0, y0));

  double nr_y0[4] = { 0.0, 1.0, 1.0, 1.0 };
  double nr_dydx0[4];
  nr_derivative(0.0, nr_y0, nr_dydx0);

  BOOST_CHECK_CLOSE(nr_dydx0[1], dydx0[0], 1e-13);
  BOOST_CHECK_CLOSE(nr_dydx0[2], dydx0[1], 1e-13);
  BOOST_CHECK_CLOSE(nr_dydx0[3], dydx0[2], 1e-13);

  BOOST_CHECK_CLOSE(1.0, y0[0], 1e-13);
  dvector y(bs::numeric::ode::rk4(y0, dydx0, 0.0, 0.1, test_derivative));

  double nr_y[4];
  nr_rk4(nr_y0, nr_dydx0, 3, 0.0, 0.1, nr_y, nr_derivative);

  BOOST_CHECK_CLOSE(1.105170833333333, y[0], 1e-13);
  BOOST_CHECK_CLOSE(nr_y[1], y[0], 1e-13);
}

BOOST_AUTO_TEST_CASE(test_rkdumb)
{
  dvector y0(3);
  y0[0] = y0[1] = y0[2] = 1.0;

  double nr_y0[4] = { 0.0, 1.0, 1.0, 1.0 };

  dvector x;
  dmatrix y;
  boost::tie(x, y) = bs::numeric::ode::rkdumb(y0, 0.0, 1.0, 10, test_derivative);

  double nr_y[4*12];
  double nr_xx[11];
  nr_rkdumb(nr_y0, 3, 0.0, 1.0, 10, nr_xx, nr_y, nr_derivative);

  BOOST_CHECK_EQUAL(11, y.size1());
  BOOST_CHECK_EQUAL(3, y.size2());

  BOOST_CHECK_CLOSE(nr_y[1*12 + 11], y(10, 0), 1e-13);
}

/*
 * Local Variables:
 * coding: utf-8
 * mode: c++
 * c-basic-offset: 2
 * indent-tabs-mode: nil
 * End:
 */