testODEIntegrator.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>

#include <utl/RK4ODEIntegrator.h>
#include <utl/RK5ODEIntegrator.h>
#include <utl/SVector.h>

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>

using namespace std;
using namespace utl;
using namespace tst;

#define ASSERT_CLOSE(x, y) CPPUNIT_ASSERT(Verify<CloseTo>(x, y))
#define ASSERT_CLOSE_EPS(x, y, eps) CPPUNIT_ASSERT(Verify<CloseTo>(x, y, eps))
#define ASSERT_EQUAL(x, y) CPPUNIT_ASSERT(Verify<Equal>(x, y))


class HarmonicOscillator {
public:
  HarmonicOscillator(const double m, const double k) : fM(m), fK(k) { }

  template<typename Vector>
  bool
  operator()(const double /*x*/, const Vector& y, Vector& dYdX)
    const
  {
    // m \ddot{x} = -k x
    // y_0 = x
    // y_1 = p
    dYdX[0] = y[1] / fM;
    dYdX[1] = -fK * y[0];
    return true;
  }

  operator unsigned int() const { return 2; }

private:
  const double fM;
  const double fK;
};


class TestODEIntegrator : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestODEIntegrator);
  CPPUNIT_TEST(TestRK4);
  CPPUNIT_TEST(TestRK4SVector);
  CPPUNIT_TEST(TestRK5);
  CPPUNIT_TEST(TestRK5SVector);
  CPPUNIT_TEST(TestAdaptiveRK5);
  CPPUNIT_TEST_SUITE_END();

public:
  void
  TestRK4()
  {
    typedef double Vector[2];

    const double x0 = 0;
    const double x1 = 100;
    const double dx = 0.1;
    const Vector y0 = { 1, 0 };

    const double m = 1;
    const double k = 1;
    HarmonicOscillator ho(m, k);

    RK4ODEIntegrator<HarmonicOscillator> rk4(ho);

    RK4Iterator<HarmonicOscillator, Vector> it = rk4.Begin(x0, y0);

    const double omega = sqrt(k/m);

    const double eps = 1e-4;

    double x = x0;

    for ( ; it.GetX() < x1; it += dx) {
      ASSERT_CLOSE(it.GetX(), x);
      ASSERT_CLOSE_EPS(it.GetY()[0],  y0[0] * cos(omega * it.GetX()), eps);
      ASSERT_CLOSE_EPS(it.GetY()[1], -y0[0] * sin(omega * it.GetX()), eps);

      x += dx;
    }
  }

  void
  TestRK4SVector()
  {
    const double x0 = 0;
    const double x1 = 100;
    const double dx = 0.1;
    const double t[2] = { 1, 0 };
    const SVector<2> y0(t);

    const double m = 1;
    const double k = 1;
    HarmonicOscillator ho(m, k);

    RK4ODEIntegrator<HarmonicOscillator> rk4(ho);

    RK4Iterator<HarmonicOscillator, SVector<2> > it = rk4.Begin(x0, y0);

    const double omega = sqrt(k/m);

    const double eps = 1e-4;

    double x = x0;

    for ( ; it.GetX() < x1; it += dx) {
      ASSERT_CLOSE(it.GetX(), x);
      ASSERT_CLOSE_EPS(it.GetY()[0],  y0[0] * cos(omega * it.GetX()), eps);
      ASSERT_CLOSE_EPS(it.GetY()[1], -y0[0] * sin(omega * it.GetX()), eps);

      x += dx;
    }
  }

  void
  TestRK5()
  {
    typedef double Vector[2];

    const double x0 = 0;
    const double x1 = 100;
    const double dx = 0.1;
    const Vector y0 = { 1, 0 };

    const double m = 1;
    const double k = 1;
    HarmonicOscillator ho(m, k);

    RK5ODEIntegrator<HarmonicOscillator> rk5(ho);

    RK5Iterator<HarmonicOscillator, Vector> it = rk5.Begin(x0, y0);

    const double omega = sqrt(k/m);

    const double eps = 1e-6;

    double x = x0;

    for ( ; it.GetX() < x1; it += dx) {
      ASSERT_CLOSE(it.GetX(), x);
      ASSERT_CLOSE_EPS(it.GetY()[0],  y0[0] * cos(omega * it.GetX()), eps);
      ASSERT_CLOSE_EPS(it.GetY()[1], -y0[0] * sin(omega * it.GetX()), eps);

      x += dx;
    }
  }

  void
  TestRK5SVector()
  {
    const double x0 = 0;
    const double x1 = 100;
    const double dx = 0.1;
    const double t[2] = { 1, 0 };
    const SVector<2> y0(t);

    const double m = 1;
    const double k = 1;
    HarmonicOscillator ho(m, k);

    RK5ODEIntegrator<HarmonicOscillator> rk5(ho);

    RK5Iterator<HarmonicOscillator, SVector<2>> it = rk5.Begin(x0, y0);

    const double omega = sqrt(k/m);

    const double eps = 1e-6;

    double x = x0;

    for ( ; it.GetX() < x1; it += dx) {
      ASSERT_CLOSE(it.GetX(), x);
      ASSERT_CLOSE_EPS(it.GetY()[0],  y0[0] * cos(omega * it.GetX()), eps);
      ASSERT_CLOSE_EPS(it.GetY()[1], -y0[0] * sin(omega * it.GetX()), eps);

      x += dx;
    }
  }

  void
  TestAdaptiveRK5()
  {
    typedef double Vector[2];

    const double x0 = 0;
    const double x1 = 100;
    const double dx = 0.1;
    const Vector y0 = { 1, 0 };

    const double m = 1;
    const double k = 1;
    HarmonicOscillator ho(m, k);

    RK5ODEIntegrator<HarmonicOscillator> rk5(ho);

    const double accuracy = 1e-5;

    AdaptiveRK5Iterator<HarmonicOscillator, Vector> it =
      rk5.AdaptiveBegin(x0, dx, y0, accuracy);

    const double omega = sqrt(k/m);

    const double eps = 1e-3;

    // try first with step dx
    // if step is too large, a smaller one gets calcualted
    // but never larger
    int count = 0;
    for ( ; it.GetX() < x1; it += dx) {
      const double x = it.GetX();
      ASSERT_CLOSE_EPS(it.GetY()[0],  y0[0] * cos(omega * x), eps);
      ASSERT_CLOSE_EPS(it.GetY()[1], -y0[0] * sin(omega * x), eps);

      ++count;
    }

    cout << "\nRK5 +=: " << count << endl;

    AdaptiveRK5Iterator<HarmonicOscillator, Vector> it2 =
      rk5.AdaptiveBegin(x0, dx, y0, accuracy);

    // run all the time with optimal step which can become
    // larger than the initial dx
    count = 0;
    for ( ; it2.GetX() < x1; ++it2) {
      const double x = it2.GetX();
      ASSERT_CLOSE_EPS(it2.GetY()[0],  y0[0] * cos(omega * x), eps);
      ASSERT_CLOSE_EPS(it2.GetY()[1], -y0[0] * sin(omega * x), eps);

      ++count;
    }

    cout << "RK5 ++: " << count << endl;
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(TestODEIntegrator);