CppAD: A C++ Algorithmic Differentiation Package 20110419
romberg_one.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ROMBERG_ONE_INCLUDED
00003 # define CPPAD_ROMBERG_ONE_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 /*
00016 $begin RombergOne$$
00017 $spell
00018         cppad.hpp
00019         bool
00020         const
00021         Cpp
00022         RombergOne
00023 $$
00024 
00025 $section One DimensionalRomberg Integration$$
00026 
00027 $index integrate, Romberg$$
00028 $index Romberg, Integrate$$
00029 
00030 $head Syntax$$
00031 $code # include <cppad/romberg_one.hpp>$$
00032 $pre
00033 $$
00034 $syntax%%r% = RombergOne(%F%, %a%, %b%, %n%, %e%)%$$
00035 
00036 
00037 $head Description$$
00038 Returns the Romberg integration estimate
00039 $latex r$$ for a one dimensional integral
00040 $latex \[
00041 r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)}
00042 \] $$
00043 
00044 $head Include$$
00045 The file $code cppad/romberg_one.hpp$$ is included by $code cppad/cppad.hpp$$
00046 but it can also be included separately with out the rest of 
00047 the $code CppAD$$ routines.
00048 
00049 $head r$$
00050 The return value $italic r$$ has prototype
00051 $syntax%
00052         %Float% %r%
00053 %$$ 
00054 It is the estimate computed by $code RombergOne$$ for the integral above.
00055 
00056 $head F$$
00057 The object $italic F$$ can be of any type, but it must support 
00058 the operation
00059 $syntax%
00060         %F%(%x%)
00061 %$$
00062 The argument $italic x$$ to $italic F$$ has prototype
00063 $syntax%
00064         const %Float% &%x%
00065 %$$
00066 The return value of $italic F$$ is a $italic Float$$ object
00067 (see description of $xref/RombergOne/Float/Float/$$ below). 
00068 
00069 $head a$$
00070 The argument $italic a$$ has prototype
00071 $syntax%
00072         const %Float% &%a%
00073 %$$ 
00074 It specifies the lower limit for the integration.
00075 
00076 $head b$$
00077 The argument $italic b$$ has prototype
00078 $syntax%
00079         const %Float% &%b%
00080 %$$ 
00081 It specifies the upper limit for the integration.
00082 
00083 $head n$$
00084 The argument $italic n$$ has prototype
00085 $syntax%
00086         size_t %n%
00087 %$$ 
00088 A total number of $latex 2^{n-1} + 1$$ evaluations of $syntax%%F%(%x%)%$$
00089 are used to estimate the integral.
00090 
00091 $head p$$
00092 The argument $italic p$$ has prototype
00093 $syntax%
00094         size_t %p%
00095 %$$ 
00096 It must be less than or equal $latex n$$
00097 and determines the accuracy order in the approximation for the integral
00098 that is returned by $code RombergOne$$. 
00099 To be specific
00100 $latex \[
00101 r = \int_a^b F(x) {\bf d} x + O \left[ (b - a) / 2^{n-1} \right]^{2(p+1)}
00102 \] $$
00103 
00104 
00105 $head e$$
00106 The argument $italic e$$ has prototype
00107 $syntax%
00108         %Float% &%e%
00109 %$$ 
00110 The input value of $italic e$$ does not matter
00111 and its output value is an approximation for the error in 
00112 the integral estimates; i.e.,
00113 $latex \[
00114         e \approx \left| r - \int_a^b F(x) {\bf d} x \right|
00115 \] $$
00116 
00117 $head Float$$
00118 The type $italic Float$$ must satisfy the conditions
00119 for a $xref/NumericType/$$ type.
00120 The routine $xref/CheckNumericType/$$ will generate an error message
00121 if this is not the case.
00122 In addition, if $italic x$$ and $italic y$$ are $italic Float$$ objects,
00123 $syntax%
00124         %x% < %y%
00125 %$$     
00126 returns the $code bool$$ value true if $italic x$$ is less than 
00127 $italic y$$ and false otherwise.
00128 
00129 $children%
00130         example/romberg_one.cpp
00131 %$$
00132 $head Example$$
00133 $comment%
00134         example/romberg_one.cpp
00135 %$$
00136 The file
00137 $xref/RombergOne.cpp/$$
00138 contains an example and test a test of using this routine.
00139 It returns true if it succeeds and false otherwise.
00140 
00141 $head Source Code$$
00142 The source code for this routine is in the file
00143 $code cppad/romberg_one.hpp$$.
00144 
00145 $end
00146 */
00147 
00148 # include <cppad/check_numeric_type.hpp>
00149 # include <cppad/local/cppad_assert.hpp>
00150 # include <cppad/vector.hpp>
00151 
00152 namespace CppAD { // BEGIN CppAD namespace
00153 
00154 template <class Fun, class Float>
00155 Float RombergOne(
00156         Fun           &F , 
00157         const Float   &a , 
00158         const Float   &b , 
00159         size_t         n , 
00160         size_t         p ,
00161         Float         &e )
00162 {
00163         size_t ipow2 = 1;
00164         size_t k, i;
00165         Float pow2, sum, x; 
00166 
00167         Float  zero  = Float(0);
00168         Float  two   = Float(2);
00169 
00170         // check specifications for a NumericType
00171         CheckNumericType<Float>();
00172 
00173         CPPAD_ASSERT_KNOWN(
00174                 n >= 2,
00175                 "RombergOne: n must be greater than or equal 2"
00176         );
00177         CppAD::vector<Float> r(n);
00178 
00179         //  set r[i] = trapazoidal rule with 2^i intervals in [a, b]
00180         r[0]  = ( F(a) + F(b) ) * (b - a) / two; 
00181         for(i = 1; i < n; i++)
00182         {       ipow2 *= 2;
00183                 // there must be a conversion from int to any numeric type
00184                 pow2   = Float(int(ipow2)); 
00185                 sum    = zero;
00186                 for(k = 1; k < ipow2; k += 2)
00187                 {       // start = a + (b-a)/pow2, increment = 2*(b-a)/pow2
00188                         x    = ( (pow2 - Float(int(k))) * a + k * b ) / pow2;
00189                         sum  = sum + F(x);
00190                 }
00191                 // combine function evaluations in sum with those in T[i-1]
00192                 r[i] = r[i-1] / two + sum * (b - a) / pow2;
00193         }
00194 
00195         // now compute the higher order estimates
00196         size_t ipow4    = 1;   // order of accuract for previous estimate
00197         Float pow4, pow4minus;
00198         for(i = 0; i < p; i++)
00199         {       // compute estimate accurate to O[ step^(2*(i+1)) ]
00200                 // put resutls in r[n-1], r[n-2], ... , r[n-i+1]
00201                 ipow4    *= 4;
00202                 pow4      = Float(int(ipow4));
00203                 pow4minus = Float(ipow4-1);
00204                 for(k = n-1; k > i; k--)
00205                         r[k] = ( pow4 * r[k] - r[k-1] ) / pow4minus;
00206         }
00207 
00208         // error estimate for r[n]
00209         e = r[n-1] - r[n-2];
00210         if( e < zero )
00211                 e = - e;
00212         return r[n-1];
00213 }
00214 
00215 } // END CppAD namespace
00216 
00217 # endif