/home/coin/SVN-release/CoinAll-1.1.0/cppad/cppad/romberg_one.hpp

Go to the documentation of this file.

# ifndef CPPAD_ROMBERG_ONE_INCLUDED
# define CPPAD_ROMBERG_ONE_INCLUDED

/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 Bradley M. Bell

CppAD is distributed under multiple licenses. This distribution is under
the terms of the 
                    Common Public License Version 1.0.

A copy of this license is included in the COPYING file of this distribution.
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
-------------------------------------------------------------------------- */
/*
$begin RombergOne$$
$spell
        cppad.hpp
        bool
        const
        Cpp
        RombergOne
$$

$section One DimensionalRomberg Integration$$

$index integrate, Romberg$$
$index Romberg, Integrate$$

$head Syntax$$
$code # include <cppad/romberg_one.hpp>$$
$pre
$$
$syntax%%r% = RombergOne(%F%, %a%, %b%, %n%, %e%)%$$


$head Description$$
Returns the Romberg integration estimate
$latex r$$ for a one dimensional integral
$latex \[
r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)}
\] $$

$head Include$$
The file $code cppad/romberg_one.hpp$$ is included by $code cppad/cppad.hpp$$
but it can also be included separately with out the rest of 
the $code CppAD$$ routines.

$head r$$
The return value $italic r$$ has prototype
$syntax%
        %Float% %r%
%$$ 
It is the estimate computed by $code RombergOne$$ for the integral above.

$head F$$
The object $italic F$$ can be of any type, but it must support 
the operation
$syntax%
        %F%(%x%)
%$$
The argument $italic x$$ to $italic F$$ has prototype
$syntax%
        const %Float% &%x%
%$$
The return value of $italic F$$ is a $italic Float$$ object
(see description of $xref/RombergOne/Float/Float/$$ below). 

$head a$$
The argument $italic a$$ has prototype
$syntax%
        const %Float% &%a%
%$$ 
It specifies the lower limit for the integration.

$head b$$
The argument $italic b$$ has prototype
$syntax%
        const %Float% &%b%
%$$ 
It specifies the upper limit for the integration.

$head n$$
The argument $italic n$$ has prototype
$syntax%
        size_t %n%
%$$ 
A total number of $latex 2^{n-1} + 1$$ evaluations of $syntax%%F%(%x%)%$$
are used to estimate the integral.

$head p$$
The argument $italic p$$ has prototype
$syntax%
        size_t %p%
%$$ 
It must be less than or equal $latex n$$
and determines the accuracy order in the approximation for the integral
that is returned by $code RombergOne$$. 
To be specific
$latex \[
r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)}
\] $$


$head e$$
The argument $italic e$$ has prototype
$syntax%
        %Float% &%e%
%$$ 
The input value of $italic e$$ does not matter
and its output value is an approximation for the error in 
the integral estimates; i.e.,
$latex \[
        e \approx \left| r - \int_a^b F(x) {\bf d} x \right|
\] $$

$head Float$$
The type $italic Float$$ must satisfy the conditions
for a $xref/NumericType/$$ type.
The routine $xref/CheckNumericType/$$ will generate an error message
if this is not the case.
In addition, if $italic x$$ and $italic y$$ are $italic Float$$ objects,
$syntax%
        %x% < %y%
%$$     
returns the $code bool$$ value true if $italic x$$ is less than 
$italic y$$ and false otherwise.

$children%
        example/romberg_one.cpp
%$$
$head Example$$
$comment%
        example/romberg_one.cpp
%$$
The file
$xref/RombergOne.cpp/$$
contains an example and test a test of using this routine.
It returns true if it succeeds and false otherwise.

$head Source Code$$
The source code for this routine is in the file
$code cppad/romberg_one.hpp$$.

$end
*/

# include <cppad/check_numeric_type.hpp>
# include <cppad/local/cppad_assert.hpp>
# include <cppad/vector.hpp>

namespace CppAD { // BEGIN CppAD namespace

template <class Fun, class Float>
Float RombergOne(
        Fun           &F , 
        const Float   &a , 
        const Float   &b , 
        size_t         n , 
        size_t         p ,
        Float         &e )
{
        size_t ipow2 = 1;
        size_t k, i;
        Float pow2, sum, x; 

        Float  zero  = Float(0);
        Float  two   = Float(2);

        // check specifications for a NumericType
        CheckNumericType<Float>();

        CPPAD_ASSERT_KNOWN(
                n >= 2,
                "RombergOne: n must be greater than or equal 2"
        );
        CppAD::vector<Float> r(n);

        //  set r[i] = trapazoidal rule with 2^i intervals in [a, b]
        r[0]  = ( F(a) + F(b) ) * (b - a) / two; 
        for(i = 1; i < n; i++)
        {       ipow2 *= 2;
                // there must be a conversion from int to any numeric type
                pow2   = Float(int(ipow2)); 
                sum    = zero;
                for(k = 1; k < ipow2; k += 2)
                {       // start = a + (b-a)/pow2, increment = 2*(b-a)/pow2
                        x    = ( (pow2 - Float(int(k))) * a + k * b ) / pow2;
                        sum  = sum + F(x);
                }
                // combine function evaluations in sum with those in T[i-1]
                r[i] = r[i-1] / two + sum * (b - a) / pow2;
        }

        // now compute the higher order estimates
        size_t ipow4    = 1;   // order of accuract for previous estimate
        Float pow4, pow4minus;
        for(i = 0; i < p; i++)
        {       // compute estimate accurate to O[ step^(2*(i+1)) ]
                // put resutls in r[n-1], r[n-2], ... , r[n-i+1]
                ipow4    *= 4;
                pow4      = Float(int(ipow4));
                pow4minus = Float(ipow4-1);
                for(k = n-1; k > i; k--)
                        r[k] = ( pow4 * r[k] - r[k-1] ) / pow4minus;
        }

        // error estimate for r[n]
        e = r[n-1] - r[n-2];
        if( e < zero )
                e = - e;
        return r[n-1];
}

} // END CppAD namespace

# endif

---