CppAD: A C++ Algorithmic Differentiation Package 20110419
|
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