\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

romberg_mul.cpp

View page source

Multi-Dimensional Romberg Integration: Example and Test

# include <cppad/utility/romberg_mul.hpp>
# include <cppad/utility/vector.hpp>
# include <cppad/utility/near_equal.hpp>


namespace {

    class TestFun {
    private:
        const CppAD::vector<size_t> deg;
    public:
        // constructor
        TestFun(const CppAD::vector<size_t> deg_)
        : deg(deg_)
        { }

        // function F(x) = x[0]^deg[0] * x[1]^deg[1]
        double operator () (const CppAD::vector<double> &x)
        {   size_t i;
            double   f = 1;
            for(i = 0; i < deg[0]; i++)
                f *= x[0];
            for(i = 0; i < deg[1]; i++)
                f *= x[1];
            return f;
        }
    };

}

bool RombergMul(void)
{   bool ok = true;
    size_t i;
    size_t k;

    CppAD::vector<size_t> deg(2);
    deg[0] = 5;
    deg[1] = 3;
    TestFun F(deg);

    CppAD::RombergMul<
        TestFun              ,
        CppAD::vector<size_t>,
        CppAD::vector<double>,
        2                    > RombergMulTest;

    // arguments to RombergMul
    CppAD::vector<double> a(2);
    CppAD::vector<double> b(2);
    CppAD::vector<size_t> n(2);
    CppAD::vector<size_t> p(2);
    for(i = 0; i < 2; i++)
    {   a[i] = 0.;
        b[i] = 1.;
    }
    n[0] = 4;
    n[1] = 3;
    double r, e;

    // int_a1^b1 dx1 int_a0^b0 F(x0,x1) dx0
    //    = [ b0^(deg[0]+1) - a0^(deg[0]+1) ] / (deg[0]+1)
    //    * [ b1^(deg[1]+1) - a1^(deg[1]+1) ] / (deg[1]+1)
    double bpow = 1.;
    double apow = 1.;
    for(i = 0; i <= deg[0]; i++)
    {   bpow *= b[0];
        apow *= a[0];
    }
    double check = (bpow - apow) / double(deg[0]+1);
    bpow = 1.;
    apow = 1.;
    for(i = 0; i <= deg[1]; i++)
    {   bpow *= b[1];
        apow *= a[1];
    }
    check *= (bpow - apow) / double(deg[1]+1);

    double step = (b[1] - a[1]) / exp(log(2.)*double(n[1]-1));
    double spow = 1;
    for(k = 0; k <= n[1]; k++)
    {   spow = spow * step * step;
        double bnd = 3 * double(deg[1] + 1) * spow;

        for(i = 0; i < 2; i++)
            p[i] = k;
        r    = RombergMulTest(F, a, b, n, p, e);

        ok  &= e < bnd;
        ok  &= CppAD::NearEqual(check, r, 0., e);

    }

    return ok;
}