CppAD: A C++ Algorithmic Differentiation Package  20130102
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-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse 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 $icode%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 $icode r$$ has prototype
00051 $codei%
00052      %Float% %r%
00053 %$$ 
00054 It is the estimate computed by $code RombergOne$$ for the integral above.
00055 
00056 $head F$$
00057 The object $icode F$$ can be of any type, but it must support 
00058 the operation
00059 $codei%
00060      %F%(%x%)
00061 %$$
00062 The argument $icode x$$ to $icode F$$ has prototype
00063 $codei%
00064      const %Float% &%x%
00065 %$$
00066 The return value of $icode F$$ is a $icode Float$$ object
00067 (see description of $cref/Float/RombergOne/Float/$$ below). 
00068 
00069 $head a$$
00070 The argument $icode a$$ has prototype
00071 $codei%
00072      const %Float% &%a%
00073 %$$ 
00074 It specifies the lower limit for the integration.
00075 
00076 $head b$$
00077 The argument $icode b$$ has prototype
00078 $codei%
00079      const %Float% &%b%
00080 %$$ 
00081 It specifies the upper limit for the integration.
00082 
00083 $head n$$
00084 The argument $icode n$$ has prototype
00085 $codei%
00086      size_t %n%
00087 %$$ 
00088 A total number of $latex 2^{n-1} + 1$$ evaluations of $icode%F%(%x%)%$$
00089 are used to estimate the integral.
00090 
00091 $head p$$
00092 The argument $icode p$$ has prototype
00093 $codei%
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 $icode e$$ has prototype
00107 $codei%
00108      %Float% &%e%
00109 %$$ 
00110 The input value of $icode 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 $icode Float$$ must satisfy the conditions
00119 for a $cref NumericType$$ type.
00120 The routine $cref CheckNumericType$$ will generate an error message
00121 if this is not the case.
00122 In addition, if $icode x$$ and $icode y$$ are $icode Float$$ objects,
00123 $codei%
00124      %x% < %y%
00125 %$$     
00126 returns the $code bool$$ value true if $icode x$$ is less than 
00127 $icode 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 $cref romberg_one.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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines