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