Prev Next Index-> contents reference index search external Up-> CppAD utility Runge45 runge45_2.cpp CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix utility-> ErrorHandler NearEqual speed_test SpeedTest time_test test_boolofvoid NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union sparse_rc sparse_rcv Runge45-> runge45_1.cpp runge45_2.cpp runge45_2.cpp Headings

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
Runge45: Example and Test
Define $X : \B{R} \times \B{R} \rightarrow \B{R}^n$ by $$X_j (b, t) = b \left( \sum_{k=0}^j t^k / k ! \right)$$ for $j = 0 , \ldots , n-1$. It follows that $$\begin{array}{rcl} X_j (b, 0) & = & b \\ \partial_t X_j (b, t) & = & b \left( \sum_{k=0}^{j-1} t^k / k ! \right) \\ \partial_t X_j (b, t) & = & \left\{ \begin{array}{ll} 0 & {\rm if} \; j = 0 \\ X_{j-1} (b, t) & {\rm otherwise} \end{array} \right. \end{array}$$ For a fixed $t_f$, we can use Runge45 to define $f : \B{R} \rightarrow \B{R}^n$ as an approximation for $f(b) = X(b, t_f )$. We can then compute $f^{(1)} (b)$ which is an approximation for $$\partial_b X(b, t_f ) = \sum_{k=0}^j t_f^k / k !$$  # include <cstddef> // for size_t # include <limits> // for machine epsilon # include <cppad/cppad.hpp> // for all of CppAD namespace { template <class Scalar> class Fun { public: // constructor Fun(void) { } // set return value to X'(t) void Ode( const Scalar &t, const CPPAD_TESTVECTOR(Scalar) &x, CPPAD_TESTVECTOR(Scalar) &f) { size_t n = x.size(); f[0] = 0.; for(size_t k = 1; k < n; k++) f[k] = x[k-1]; } }; } bool runge_45_2(void) { typedef CppAD::AD<double> Scalar; using CppAD::NearEqual; bool ok = true; // initial return value size_t j; // temporary indices size_t n = 5; // number components in X(t) and order of method size_t M = 2; // number of Runge45 steps in [ti, tf] Scalar ad_ti = 0.; // initial time Scalar ad_tf = 2.; // final time // value of independent variable at which to record operations CPPAD_TESTVECTOR(Scalar) ad_b(1); ad_b[0] = 1.; // declare b to be the independent variable Independent(ad_b); // object to evaluate ODE Fun<Scalar> ad_F; // xi = X(0) CPPAD_TESTVECTOR(Scalar) ad_xi(n); for(j = 0; j < n; j++) ad_xi[j] = ad_b[0]; // compute Runge45 approximation for X(tf) CPPAD_TESTVECTOR(Scalar) ad_xf(n), ad_e(n); ad_xf = CppAD::Runge45(ad_F, M, ad_ti, ad_tf, ad_xi, ad_e); // stop recording and use it to create f : b -> xf CppAD::ADFun<double> f(ad_b, ad_xf); // evaluate f(b) CPPAD_TESTVECTOR(double) b(1); CPPAD_TESTVECTOR(double) xf(n); b[0] = 1.; xf = f.Forward(0, b); // check that f(b) = X(b, tf) double tf = Value(ad_tf); double term = 1; double sum = 0; double eps = 10. * CppAD::numeric_limits<double>::epsilon(); for(j = 0; j < n; j++) { sum += term; ok &= NearEqual(xf[j], b[0] * sum, eps, eps); term *= tf; term /= double(j+1); } // evalute f'(b) CPPAD_TESTVECTOR(double) d_xf(n); d_xf = f.Jacobian(b); // check that f'(b) = partial of X(b, tf) w.r.t b term = 1; sum = 0; for(j = 0; j < n; j++) { sum += term; ok &= NearEqual(d_xf[j], sum, eps, eps); term *= tf; term /= double(j+1); } return ok; } 
Input File: example/general/runge45_2.cpp