CppAD: A C++ Algorithmic Differentiation Package 20110419
romberg_mul.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_ROMBERG_MUL_INCLUDED
00003 # define CPPAD_ROMBERG_MUL_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 RombergMul$$
00017 $spell
00018         cppad.hpp
00019         bool
00020         const
00021         Cpp
00022         RombergMulMul
00023 $$
00024 
00025 $section Multi-dimensional Romberg Integration$$
00026 
00027 $index integrate, multi-dimensional Romberg$$
00028 $index Romberg, multi-dimensional integrate$$
00029 $index multi, dimensional Romberg integration$$
00030 $index dimension, multi Romberg integration$$
00031 
00032 $head Syntax$$
00033 $code # include <cppad/romberg_mul.hpp>$$
00034 $pre
00035 $$
00036 $syntax%RombergMul<%Fun%, %SizeVector%, %FloatVector%, %m%> %R%$$
00037 $pre
00038 $$
00039 $syntax%%r% = %R%(%F%, %a%, %b%, %n%, %p%, %e%)%$$
00040 
00041 
00042 $head Description$$
00043 Returns the Romberg integration estimate
00044 $latex r$$ for the multi-dimensional integral
00045 $latex \[
00046 r = 
00047 \int_{a[0]}^{b[0]} \cdots \int_{a[m-1]}^{b[m-1]}
00048 \; F(x) \;
00049 {\bf d} x_0 \cdots {\bf d} x_{m-1} 
00050 \; + \; 
00051 \sum_{i=0}^{m-1} 
00052 O \left[ ( b[i] - a[i] ) / 2^{n[i]-1} \right]^{2(p[i]+1)}
00053 \] $$
00054 
00055 $head Include$$
00056 The file $code cppad/romberg_mul.hpp$$ is included by $code cppad/cppad.hpp$$
00057 but it can also be included separately with out the rest of 
00058 the $code CppAD$$ routines.
00059 
00060 $head m$$
00061 The template parameter $italic m$$ must be convertible to a $code size_t$$ 
00062 object with a value that can be determined at compile time; for example
00063 $code 2$$.
00064 It determines the dimension of the domain space for the integration.
00065 
00066 $head r$$
00067 The return value $italic r$$ has prototype
00068 $syntax%
00069         %Float% %r%
00070 %$$ 
00071 It is the estimate computed by $code RombergMul$$ for the integral above
00072 (see description of $xref/RombergMul/Float/Float/$$ below). 
00073 
00074 $head F$$
00075 The object $italic F$$ has the prototype
00076 $syntax%
00077         %Fun% &%F%
00078 %$$
00079 It must support the operation
00080 $syntax%
00081         %F%(%x%)
00082 %$$
00083 The argument $italic x$$ to $italic F$$ has prototype
00084 $syntax%
00085         const %Float% &%x%
00086 %$$
00087 The return value of $italic F$$ is a $italic Float$$ object
00088 
00089 $head a$$
00090 The argument $italic a$$ has prototype
00091 $syntax%
00092         const %FloatVector% &%a%
00093 %$$ 
00094 It specifies the lower limit for the integration
00095 (see description of $xref/RombergMul/FloatVector/FloatVector/$$ below). 
00096 
00097 $head b$$
00098 The argument $italic b$$ has prototype
00099 $syntax%
00100         const %FloatVector% &%b%
00101 %$$ 
00102 It specifies the upper limit for the integration.
00103 
00104 $head n$$
00105 The argument $italic n$$ has prototype
00106 $syntax%
00107         const %SizeVector% &%n%
00108 %$$ 
00109 A total number of $latex 2^{n[i]-1} + 1$$ 
00110 evaluations of $syntax%%F%(%x%)%$$ are used to estimate the integral
00111 with respect to $latex {\bf d} x_i$$.
00112 
00113 $head p$$
00114 The argument $italic p$$ has prototype
00115 $syntax%
00116         const %SizeVector% &%p%
00117 %$$ 
00118 For $latex i = 0 , \ldots , m-1$$,
00119 $latex n[i]$$ determines the accuracy order in the 
00120 approximation for the integral 
00121 that is returned by $code RombergMul$$.
00122 The values in $italic p$$ must be less than or equal $italic n$$; i.e.,
00123 $syntax%%p%[%i%] <= %n%[%i%]%$$.
00124 
00125 $head e$$
00126 The argument $italic e$$ has prototype
00127 $syntax%
00128         %Float% &%e%
00129 %$$ 
00130 The input value of $italic e$$ does not matter
00131 and its output value is an approximation for the absolute error in 
00132 the integral estimate.
00133 
00134 $head Float$$
00135 The type $italic Float$$ is defined as the type of the elements of
00136 $xref/RombergMul/FloatVector/FloatVector/$$.
00137 The type $italic Float$$ must satisfy the conditions
00138 for a $xref/NumericType/$$ type.
00139 The routine $xref/CheckNumericType/$$ will generate an error message
00140 if this is not the case.
00141 In addition, if $italic x$$ and $italic y$$ are $italic Float$$ objects,
00142 $syntax%
00143         %x% < %y%
00144 %$$     
00145 returns the $code bool$$ value true if $italic x$$ is less than 
00146 $italic y$$ and false otherwise.
00147 
00148 $head FloatVector$$
00149 The type $italic FloatVector$$ must be a $xref/SimpleVector/$$ class.
00150 The routine $xref/CheckSimpleVector/$$ will generate an error message
00151 if this is not the case.
00152 
00153 
00154 $children%
00155         example/romberg_mul.cpp
00156 %$$
00157 $head Example$$
00158 $comment%
00159         example/romberg_mul.cpp
00160 %$$
00161 The file
00162 $xref/RombergMul.cpp/$$
00163 contains an example and test a test of using this routine.
00164 It returns true if it succeeds and false otherwise.
00165 
00166 $head Source Code$$
00167 The source code for this routine is in the file
00168 $code cppad/romberg_mul.hpp$$.
00169 
00170 $end
00171 */
00172 
00173 # include <cppad/romberg_one.hpp>
00174 # include <cppad/check_numeric_type.hpp>
00175 # include <cppad/check_simple_vector.hpp>
00176 
00177 namespace CppAD { // BEGIN CppAD namespace
00178 
00179 template <class Fun, class FloatVector>
00180 class SliceLast {
00181         typedef typename FloatVector::value_type Float;
00182 private:
00183         Fun        *F;
00184         size_t      last;
00185         FloatVector x;
00186 public:
00187         SliceLast( Fun *F_, size_t last_, const FloatVector &x_ ) 
00188         : F(F_) , last(last_), x(last + 1)
00189         {       size_t i;
00190                 for(i = 0; i < last; i++)
00191                         x[i] = x_[i];
00192         }
00193         double operator()(const Float &xlast)
00194         {       x[last] = xlast;
00195                 return (*F)(x);
00196         }
00197 };
00198 
00199 template <class Fun, class SizeVector, class FloatVector, class Float>
00200 class IntegrateLast {
00201 private:
00202         Fun                 *F; 
00203         const size_t        last;
00204         const FloatVector   a; 
00205         const FloatVector   b; 
00206         const SizeVector    n; 
00207         const SizeVector    p; 
00208         Float               esum;
00209         size_t              ecount;
00210 
00211 public:
00212         IntegrateLast(
00213                 Fun                *F_    , 
00214                 size_t              last_ ,
00215                 const FloatVector  &a_    , 
00216                 const FloatVector  &b_    , 
00217                 const SizeVector   &n_    , 
00218                 const SizeVector   &p_    ) 
00219         : F(F_) , last(last_), a(a_) , b(b_) , n(n_) , p(p_) 
00220         { }             
00221         Float operator()(const FloatVector           &x)
00222         {       Float r, e;
00223                 SliceLast<Fun, FloatVector           > S(F, last, x);
00224                 r     = CppAD::RombergOne(
00225                         S, a[last], b[last], n[last], p[last], e
00226                 );
00227                 esum = esum + e;
00228                 ecount++;
00229                 return r;
00230         }
00231         void ClearEsum(void)
00232         {       esum   = 0.; }
00233         Float GetEsum(void)
00234         {       return esum; }
00235 
00236         void ClearEcount(void)
00237         {       ecount   = 0; }
00238         size_t GetEcount(void)
00239         {       return ecount; }
00240 };
00241 
00242 template <class Fun, class SizeVector, class FloatVector, size_t m>
00243 class RombergMul {
00244         typedef typename FloatVector::value_type Float;
00245 public:
00246         RombergMul(void)
00247         {       }
00248         Float operator() (
00249                 Fun                 &F  , 
00250                 const FloatVector   &a  ,
00251                 const FloatVector   &b  ,
00252                 const SizeVector    &n  ,
00253                 const SizeVector    &p  ,
00254                 Float               &e  )
00255         {       Float r;
00256 
00257                 typedef IntegrateLast<
00258                         Fun         , 
00259                         SizeVector  , 
00260                         FloatVector , 
00261                         Float       > IntegrateOne;
00262 
00263                 IntegrateOne Fm1(&F, m-1, a, b, n, p);
00264                 RombergMul<
00265                         IntegrateOne, 
00266                         SizeVector  ,
00267                         FloatVector ,
00268                         m-1         > RombergMulM1;
00269 
00270                 Fm1.ClearEsum();
00271                 Fm1.ClearEcount();
00272 
00273                 r  = RombergMulM1(Fm1, a, b, n, p, e);
00274 
00275                 size_t i, j;
00276                 Float prod = 1;
00277                 size_t pow2 = 1;
00278                 for(i = 0; i < m-1; i++)
00279                 {       prod *= (b[i] - a[i]);
00280                         for(j = 0; j < (n[i] - 1); j++)
00281                                 pow2 *= 2;
00282                 }
00283                 assert( Fm1.GetEcount() == (pow2+1) );
00284                 
00285                 e = e + Fm1.GetEsum() * prod / Fm1.GetEcount();
00286 
00287                 return r;
00288         }
00289 };
00290 
00291 template <class Fun, class SizeVector, class FloatVector>
00292 class RombergMul <Fun, SizeVector, FloatVector, 1> {
00293         typedef typename FloatVector::value_type Float;
00294 public:
00295         Float operator() (
00296                 Fun                 &F  , 
00297                 const FloatVector   &a  ,
00298                 const FloatVector   &b  ,
00299                 const SizeVector    &n  ,
00300                 const SizeVector    &p  ,
00301                 Float               &e  )
00302         {       Float r;
00303                 typedef IntegrateLast<
00304                         Fun         , 
00305                         SizeVector  , 
00306                         FloatVector , 
00307                         Float       > IntegrateOne;
00308 
00309                 // check simple vector class specifications
00310                 CheckSimpleVector<Float, FloatVector>();
00311 
00312                 // check numeric type specifications
00313                 CheckNumericType<Float>();
00314 
00315                 IntegrateOne F0(&F, 0, a, b, n, p);
00316 
00317                 F0.ClearEsum();
00318                 F0.ClearEcount();
00319 
00320                 r  = F0(a);
00321 
00322                 assert( F0.GetEcount() == 1 );
00323                 e = F0.GetEsum();
00324 
00325                 return r;
00326         }
00327 };
00328 
00329 } // END CppAD namespace
00330 
00331 # endif