test/test_qtrap.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 <cmath>
#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
{

    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("trapezoidal");

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

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

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

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

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

    template<>
    template<>
    void testobject::test<6>()
    {
        int info;
        double I = NR::qtrap(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>()
    {
        int info;
        const double left = 1e-5;
        const double right = 1.0;
        double I = NR::qtrap(recip, left, 1.0, 1e-12, info);
        double expected = std::log(right) - std::log(left);
        rDebug("I - %e = %e", expected, I - expected);
        ensure(fabs(I - expected) <= tol);
    }

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

};

---