source: trunk/test_more/runge_45.cpp @ 2354

Last change on this file since 2354 was 1370, checked in by bradbell, 11 years ago

trunk: Fix svn_add_id.sh and use it set Id property for some missed files.

  • Property svn:keywords set to Id
File size: 3.5 KB
Line 
1/* $Id: runge_45.cpp 1370 2009-05-31 05:31:50Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14Old OdeRunge example now used just for valiadation testing of Runge45
15*/
16
17# include <cppad/cppad.hpp>
18# include <iostream>
19# include <cassert>
20
21namespace { // BEGIN Empty namespace
22        class TestFun {
23        public:
24                TestFun(const CPPAD_TEST_VECTOR< CppAD::AD<double> > &w_)
25                {       w.resize( w_.size() );
26                        w = w_;
27                }
28                void Ode(
29                        const CppAD::AD<double>                      &t, 
30                        const CPPAD_TEST_VECTOR< CppAD::AD<double> > &x, 
31                        CPPAD_TEST_VECTOR< CppAD::AD<double> >       &f) 
32                {
33                        using CppAD::exp;
34       
35                        size_t n = x.size();
36       
37                        size_t i;
38                        f[0]  = 0.;
39                        for(i = 1; i < n-1; i++)
40                                f[i] = w[i] * x[i-1];
41       
42                        f[n-1] = x[0] * x[1];
43                }
44        private:
45                CPPAD_TEST_VECTOR< CppAD::AD<double> > w;
46        };
47} // END Empty namespace
48
49bool Runge45(void)
50{       bool ok = true;
51
52        using namespace CppAD;
53
54        size_t i;
55        size_t j;
56        size_t k;
57
58        size_t n = 6;
59        size_t m = n - 1;
60
61        CPPAD_TEST_VECTOR< AD<double> > x(n);
62        AD<double>                t0    = 0.;
63        AD<double>                t1    = 2.;
64        size_t                    nstep = 2;
65
66        // vector of independent variables
67        CPPAD_TEST_VECTOR< AD<double> > w(m);
68        for(i = 0; i < m; i++)
69                w[i] = double(i);
70        Independent(w);
71
72        // construct function object using independent variables
73        TestFun fun(w);
74
75        // initial value of x
76        CPPAD_TEST_VECTOR< AD<double> > x0(n);
77        for(i = 0; i < n; i++)
78                x0[i] = 0.;
79        x0[0] = exp( w[0] );
80
81        // solve the differential equation
82        x = Runge45(fun, nstep, t0, t1, x0);
83
84        // create f : w -> x and vectors for evaluating derivatives
85        ADFun<double> f(w, x);
86        CPPAD_TEST_VECTOR<double> q( f.Domain() );
87        CPPAD_TEST_VECTOR<double> r( f.Range() );
88
89        // for i < n-1,
90        // x[i](2) = exp( w[0] ) * (w[1] / 1) * ... * (w[i] / i) * 2^i
91        AD<double> xi2 = exp(w[0]);
92        for(i = 0; i < n-1; i++)
93        {       ok &= NearEqual(x[i],  xi2, 1e-14, 1e-14);
94                if( i < n-2 )
95                        xi2 *= w[i+1] * 2. / double(i+1);
96        }
97
98        // x[n-1](2) = exp(2 * w[0]) * w[1] * 2^2 / 2
99        xi2 = exp(2. * w[0]) * w[1] * 2.;
100        ok &= NearEqual(x[n-1], xi2, 1e-14, 1e-14);
101
102        // the partial of x[i](2) with respect to w[j] is
103        //      x[i](2) / w[j] if 0 < j <= i < n-1
104        //      x[i](2)        if j == 0 and i < n-1
105        //      2*x[i](2)      if j == 0 and i = n-1
106        //      x[i](2) / w[j] if j == 1 and i = n-1
107        //      zero           otherwise
108
109        for(i = 0; i < n-1; i++)
110        {       // compute partials of x[i]
111                for(k = 0; k < n; k++)
112                        r[k] = 0.;
113                r[i] = 1.;
114                q    = f.Reverse(1,r);
115
116                for(j = 0; j < m; j++)
117                {       // check partial of x[i] w.r.t w[j]
118                        if (j == 0 )
119                                ok &= NearEqual(q[j], x[i], 1e-14, 1e-14);
120                        else if( j <= i  ) 
121                                ok &= NearEqual(
122                                        q[j], x[i]/w[j], 1e-14, 1e-14);
123                        else    ok &= NearEqual(q[j], 0., 1e-14, 1e-14);
124                }
125        }
126
127        // compute partials of x[n-1]
128        i = n-1;
129        for(k = 0; k < n; k++)
130                r[k] = 0.;
131        r[i] = 1.;
132        q    = f.Reverse(1,r);
133
134        for(j = 0; j < m; j++)
135        {       // check partial of x[n-1] w.r.t w[j]
136                if (j == 0 )
137                        ok &= NearEqual(q[j], 2.*x[i], 1e-14, 1e-14);
138                else if( j == 1  ) 
139                        ok &= NearEqual(
140                                q[j], x[i]/w[1], 1e-14, 1e-14);
141                else    ok &= NearEqual(q[j], 0., 1e-14, 1e-14);
142        }
143
144        return ok;
145}
Note: See TracBrowser for help on using the repository browser.