test/test_romberg.cpp

Go to the documentation of this file.

//
// C++ Implementation: test_romberg
//
// Description: unit tests for romberg integration
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <math.h>
#include <tut.h>
#include "romberg.h"

static double zero(double x) {
    return 0.0;
}

static double one(double x) {
    return 1.0;
}

static double ramp(double x) {
    return x;
}

static double sawtooth(double x) {
    return x-floor(x);
}

static double recip(double x) {
    return 1.0/x;
}

namespace
{
    // Example used in Gander Gautschi paper
    struct gg_example {
        gg_example() : num_eval(0)
        {
        }
        double operator()(double x) {
            num_eval ++;
            if (x < 1.0)
                return x + 1.0;
            else if (x < 3.0)
                return 3.0 - x;
            else
                return 2.0;
        }
        int num_eval;
    };

    struct romberg_data
    {
        romberg_data()
            : pi(::acos(-1.0)),
              tol(1e-8)
        {
        }
        const double pi;
        const double tol;
    };

    typedef tut::test_group<romberg_data> testgroup;
    typedef testgroup::object testobject;
    testgroup romberg_group("romberg");

    template<>
    template<>
    void testobject::test<1>()
    {
        int info;
        double I = NR::qromb(zero, 0.0, 1.0, 1e-12, info);
        ensure(I == 0.0);
    }

    template<>
    template<>
    void testobject::test<2>()
    {
        int info;
        double I = NR::qromb(one, 0.0, 1.0, 1e-12, info);
        ensure(I == 1.0);
    }

    template<>
    template<>
    void testobject::test<3>()
    {
        int info;
        double I = NR::qromb(ramp, 0.0, 1.0, 1e-12, info);
        ensure(I == 0.5);
    }

    template<>
    template<>
    void testobject::test<4>()
    {
        int info;
        double I = NR::qromb(::sin, 0.0, pi, 1e-12, info);
        ensure(fabs(I - 2.0) <= tol);
    }

    template<>
    template<>
    void testobject::test<5>()
    {
        int info;
        double I = NR::qromb(::sin, 0.0, 2.0*pi, 1e-12, info);
        ensure(fabs(I) <= tol);
    }

    template<>
    template<>
    void testobject::test<6>()
    {
        // discontinous function
        int info;
        double I = NR::qromb(sawtooth, 0.0, 3.0, 1e-12, info);
        rDebug("I - 1.5 = %e", I - 1.5);
        ensure(fabs(I - 1.5) <= tol);
    }

    template<>
    template<>
    void testobject::test<7>()
    {
        // function with very big derivative
        int info;
        const double left = 1e-5;
        const double right = 1.0;
        double I = NR::qromb(recip, left, right, 1e-12, info);
        double expected = ::log(right) - ::log(left);
        rDebug("I - %e = %e", expected, I - expected);
        ensure(fabs(I - expected) <= tol);
    }

    template<>
    template<>
    void testobject::test<8>()
    {
        int info;
        gg_example f;
        double I = NR::qromb(f, 0.0, 5.0, 1e-12, info);
        double expected = 7.5;
        rDebug("I - %e = %e, num_eval=%d", expected, I - expected, f.num_eval);
        ensure(fabs((I - expected)/expected) <= 1e-5);
    }

    template<>
    template<>
    void testobject::test<9>()
    {
        int info;
        double I = NR::qromb(sin, pi, 0.0, 1e-12, info);
        ensure(fabs(I + 2.0) <= tol);
    }

};

---