source: trunk/test_more/rosen_34.cpp @ 1559

Last change on this file since 1559 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: rosen_34.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 OdeImplicit example now used just for valiadation testing of Rosen34
15*/
16// BEGIN PROGRAM
17
18# include <cppad/cppad.hpp>
19
20# include <iostream>
21# include <cassert>
22
23/*
24Case where
25x[0](0) = 1, x[0]'(t) = - w[0] * x[0](t)
26x[1](0) = 1, x[1]'(t) = - w[1] * x[1](t)
27x[2](0) = 0, x[2]'(t) =   w[2] * t
28
29x[0](t) = exp( - w[0] * t )
30x[1](t) = exp( - w[1] * t )
31x[2](t) = w[2] * t^2 / 2
32*/
33
34namespace {  // BEGIN Empty namespace
35        class TestFun {
36        public:
37                TestFun(const CPPAD_TEST_VECTOR< CppAD::AD<double> > &w_)
38                {       w.resize( w_.size() );
39                        w = w_;
40                }
41                void Ode(
42                        const CppAD::AD<double>                      &t, 
43                        const CPPAD_TEST_VECTOR< CppAD::AD<double> > &x, 
44                        CPPAD_TEST_VECTOR< CppAD::AD<double> >       &f) 
45                {
46                        f[0] = - w[0] * x[0];
47                        f[1] = - w[1] * x[1];
48                        f[2] =   w[2] * t;
49       
50                }
51       
52                void Ode_ind(
53                        const CppAD::AD<double>                      &t, 
54                        const CPPAD_TEST_VECTOR< CppAD::AD<double> > &x, 
55                        CPPAD_TEST_VECTOR< CppAD::AD<double> >       &f_t) 
56                {
57                        f_t[0] = 0.;
58                        f_t[1] = 0.;
59                        f_t[2] = w[2];
60       
61                }
62       
63                void Ode_dep(
64                        const CppAD::AD<double>                      &t, 
65                        const CPPAD_TEST_VECTOR< CppAD::AD<double> > &x, 
66                        CPPAD_TEST_VECTOR< CppAD::AD<double> >       &f_x) 
67                {
68                        f_x[0] = - w[0];    f_x[1] = 0.;      f_x[2] = 0.;
69                        f_x[3] = 0.;        f_x[4] = - w[1];  f_x[5] = 0.;
70                        f_x[6] = 0.;        f_x[7] = 0.;      f_x[8] = 0.;
71       
72                }
73       
74        private:
75                CPPAD_TEST_VECTOR< CppAD::AD<double> > w;
76        };
77}       // END empty namespace
78
79bool Rosen34(void)
80{       bool ok = true;
81
82        using namespace CppAD;
83
84        CPPAD_TEST_VECTOR< AD<double> > x(3);
85        CPPAD_TEST_VECTOR< AD<double> > w(3);
86        size_t         n     = 3;
87        size_t         nstep = 20;
88        AD<double>     t0    = 0.;
89        AD<double>     t1    = 1.;
90
91        // set independent variables
92        size_t i;
93        for(i = 0; i < n; i++)
94                w[i] = 100 * i + 1.;
95        Independent(w);
96
97        // construct the function object using the independent variables
98        TestFun  fun(w);
99       
100        // initial value of x
101        CPPAD_TEST_VECTOR< AD<double> > xini(3);
102        xini[0] = 1.;
103        xini[1] = 1.;
104        xini[2] = 0.;
105       
106
107        // integrate the differential equation
108        x  = Rosen34(fun, nstep, t0, t1, xini);
109
110        // create f : w -> x and vectors for evaluating derivatives
111        ADFun<double> f(w, x);
112        CPPAD_TEST_VECTOR<double> q( f.Domain() );
113        CPPAD_TEST_VECTOR<double> r( f.Range() );
114
115        // check function values
116        AD<double> x0 = exp( - w[0] * t1 );
117        ok &= NearEqual(x[0], x0, 0., 1. / (nstep * nstep) );
118
119        AD<double> x1 = exp( - w[1] * t1 );
120        ok &= NearEqual(x[1],  x1, 0., 1. / (nstep * nstep) );
121
122        AD<double> x2 = w[2] * t1 * t1 / 2.;
123        ok &= NearEqual(x[2],  x2, 1e-14, 1e-14);
124
125        // check dx[0] / dw[0]
126        for(i = 0; i < w.size(); i++)
127                q[i] = 0.;
128        q[0] = 1.;
129        r    = f.Forward(1, q);
130        ok &= NearEqual(r[0], - w[0] * x0, 0., 1. / (nstep * nstep) );
131
132        // check dx[1] / dw[1]
133        q[0] = 0.;
134        q[1] = 1.;
135        r    = f.Forward(1, q);
136        ok &= NearEqual(r[1], - w[1] * x1, 0., 1. / (nstep * nstep) );
137
138        // check dx[2] / dw[2]
139        q[1] = 0.;
140        q[2] = 1.;
141        r    = f.Forward(1, q);
142        ok &= NearEqual(r[2], x2 / w[2],  1e-14, 1e-14 );
143
144        return ok;
145}
146
147// END PROGRAM
Note: See TracBrowser for help on using the repository browser.