testIntegrator.cc

Go to the documentation of this file.

#include <cmath>
#include <iostream>
#include <utl/Integrator.h>

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

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


class SinFunctor {
public:
  SinFunctor(const double frequency, const double amplitude)
    : fFrequency(frequency), fAmplitude(amplitude) { }

  double operator()(const double x) const
  { return fAmplitude * sin(fFrequency * x); }

private:
  const double fFrequency;
  const double fAmplitude;
};


class CountedLog {
public:
  CountedLog() : fCounter(0) { }

  double operator()(const double x)
  { ++fCounter; return log(x); }

  int GetCounter() const { return fCounter; }

  void Reset() { fCounter = 0; }

private:
  int fCounter;
};


double
MyLog(const double x)
{
  return log(x);
}


class IntegratorTest : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(IntegratorTest);
  CPPUNIT_TEST(testFactory);
  CPPUNIT_TEST(testRomberg);
  CPPUNIT_TEST(testTrapezoidal);
  CPPUNIT_TEST(testSimpson);
  CPPUNIT_TEST(testCancellation);
  CPPUNIT_TEST(testZero);
  CPPUNIT_TEST(testFunction);
  CPPUNIT_TEST(testSpeed);
  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() { }

  void tearDown() { }

  void
  testFactory()
  {
    SinFunctor sinf(M_PI, 2);
    const double res =
      MakeIntegrator(sinf).GetIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res, 4/M_PI));
  }

  void
  testRomberg()
  {
    SinFunctor sinf(M_PI, 2);
    const double res1 =
      Integrator<SinFunctor>(sinf, 1e-5).GetRombergIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res1, 4/M_PI, 1e-11));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res1, 4/M_PI, 1e-12));
    const double res2 =
      Integrator<SinFunctor>(sinf, 1e-6).GetRombergIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res2, 4/M_PI, 1e-11));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res2, 4/M_PI, 1e-12));
    const double res3 =
      Integrator<SinFunctor>(sinf, 1e-5).GetRombergIntegral(1, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res3, -4/M_PI, 1e-11));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res3, -4/M_PI, 1e-12));
  }

  void
  testTrapezoidal()
  {
    SinFunctor sinf(M_PI, 2);
    const double res1 =
      Integrator<SinFunctor>(sinf, 1e-5).GetTrapezoidalIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res1, 4/M_PI, 1e-5));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res1, 4/M_PI, 1e-6));
    const double res2 =
      Integrator<SinFunctor>(sinf, 1e-6).GetTrapezoidalIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res2, 4/M_PI, 1e-6));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res2, 4/M_PI, 1e-7));
    const double res3 =
      Integrator<SinFunctor>(sinf, 1e-5).GetTrapezoidalIntegral(1, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res3, -4/M_PI, 1e-5));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res3, -4/M_PI, 1e-6));
  }

  void
  testSimpson()
  {
    SinFunctor sinf(M_PI, 2);
    const double res1 =
      Integrator<SinFunctor>(sinf, 1e-5).GetSimpsonIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res1, 4/M_PI, 1e-5));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res1, 4/M_PI, 1e-7));
    const double res2 =
      Integrator<SinFunctor>(sinf, 1e-6).GetSimpsonIntegral(0, 1);
    CPPUNIT_ASSERT(Verify<CloseTo>(res2, 4/M_PI, 1e-6));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res2, 4/M_PI, 1e-8));
    const double res3 =
      Integrator<SinFunctor>(sinf, 1e-5).GetSimpsonIntegral(1, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res3, -4/M_PI, 1e-5));
    CPPUNIT_ASSERT(Verify<Not<CloseTo> >(res1, -4/M_PI, 1e-7));
  }

  void
  testCancellation()
  {
    SinFunctor sinf(M_PI, 10);
    const double res1 =
      MakeIntegrator(sinf).GetIntegral(0, 3);
    CPPUNIT_ASSERT(Verify<CloseTo>(res1, 20/M_PI));
    const double res2 =
      MakeIntegrator(sinf).GetIntegral(1, 4);
    CPPUNIT_ASSERT(Verify<CloseTo>(res2, -20/M_PI));
  }

  void
  testZero()
  {
    SinFunctor sinf(M_PI, 1);
    const double res1 =
      MakeIntegrator(sinf).GetRombergIntegral(0, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res1, 0.));
    const double res2 =
      MakeIntegrator(sinf).GetTrapezoidalIntegral(0, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res2, 0.));
    const double res3 =
      MakeIntegrator(sinf).GetSimpsonIntegral(0, 2);
    CPPUNIT_ASSERT(Verify<CloseTo>(res3, 0.));

    const double resa =
      MakeIntegrator(sinf).GetRombergIntegral(1, 3);
    CPPUNIT_ASSERT(Verify<CloseTo>(resa, 0.));
    const double resb =
      MakeIntegrator(sinf).GetTrapezoidalIntegral(1, 3);
    CPPUNIT_ASSERT(Verify<CloseTo>(resb, 0.));
    const double resc =
      MakeIntegrator(sinf).GetSimpsonIntegral(1, 3);
    CPPUNIT_ASSERT(Verify<CloseTo>(resc, 0.));
  }

  void
  testFunction()
  {
    const double res = MakeIntegrator(MyLog).GetIntegral(1, M_E);
    CPPUNIT_ASSERT(Verify<CloseTo>(res, 1.));
  }

  void
  testSpeed()
  {
    const double eps = 1e-9;
    CountedLog clog;
    Integrator<CountedLog> clogint(clog, eps);
    const double ress = clogint.GetSimpsonIntegral(1, M_E);
    CPPUNIT_ASSERT(Verify<CloseTo>(ress, 1.));
    const int ns = clog.GetCounter();
    clog.Reset();
    const double rest = clogint.GetTrapezoidalIntegral(1, M_E);
    CPPUNIT_ASSERT(Verify<CloseTo>(rest, 1.));
    const int nt = clog.GetCounter();
    clog.Reset();
    const double resr = clogint.GetRombergIntegral(1, M_E);
    CPPUNIT_ASSERT(Verify<CloseTo>(resr, 1.));
    const int nr = clog.GetCounter();
    cerr << "\nNumber of function calls:\n"
            "    Romberg " << nr << "\n"
            "    Simpson " << ns << "\n"
            "Trapezoidal " << nt
         << endl;
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(IntegratorTest);


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: