source: trunk/test_more/ode_err_control.cpp @ 982

Last change on this file since 982 was 982, checked in by bradbell, 13 years ago

trunk: Add nan to return specifications for ode eqautions.

nan.hpp: routines for returning and detecting nan.
example/ode_err_control.cpp: a case that uses nan in ODE return.
text_more/ode_err_control.cpp: move old ode_err_control example to here.
makefile.am: add nan.hpp to distribution.
library.omh: add nan routines to library.
whats_new_07.omh: user's view of the changes.
fun_construct.hpp: fix call to Dependent in documentation.
track_new_del.hpp: improve TrackDelVec? error message.
ode_err_control.hpp: add use of nan for backing up and using smaller steps.
runge_45.hpp: use nan to communicate when step is too large.

File size: 5.0 KB
Line 
1/* --------------------------------------------------------------------------
2CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
3
4CppAD is distributed under multiple licenses. This distribution is under
5the terms of the
6                    Common Public License Version 1.0.
7
8A copy of this license is included in the COPYING file of this distribution.
9Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
10-------------------------------------------------------------------------- */
11
12# include <cstddef>                    // for size_t
13# include <cmath>                      // for exp
14# include <cppad/ode_err_control.hpp>  // CppAD::OdeErrControl
15# include <cppad/near_equal.hpp>       // CppAD::NearEqual
16# include <cppad/vector.hpp>           // CppAD::vector
17# include <cppad/runge_45.hpp>         // CppAD::Runge45
18
19/* -------------------------------------------------------------------------
20Test relative error with zero initial conditions.
21(Uses minimum step size to integrate).
22*/
23
24namespace {
25        // --------------------------------------------------------------
26        class Fun_one {
27        private:
28                 size_t n;   // dimension of the state space
29        public:
30                // constructor
31                Fun_one(size_t n_) : n(n_)
32                { } 
33
34                // given x(0) = 0
35                // solution is x_i (t) = t^(i+1)
36                void Ode(
37                        const double                &t, 
38                        const CppAD::vector<double> &x, 
39                        CppAD::vector<double>       &f)
40                {       size_t i;
41                        f[0] = 1.;
42                        for(i = 1; i < n; i++)
43                                f[i] = (i+1) * x[i-1];
44                }
45        };
46
47        // --------------------------------------------------------------
48        class Method_one {
49        private:
50                Fun_one F;
51        public:
52                // constructor
53                Method_one(size_t n_) : F(n_)
54                { }
55                void step(
56                        double ta, 
57                        double tb, 
58                        CppAD::vector<double> &xa ,
59                        CppAD::vector<double> &xb ,
60                        CppAD::vector<double> &eb )
61                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
62                }
63                size_t order(void)
64                {       return 4; }
65        };
66}
67
68bool OdeErrControl_one(void)
69{       bool   ok = true;     // initial return value
70
71        // Runge45 should yield exact results for x_i (t) = t^(i+1), i < 4
72        size_t  n = 6;       
73       
74        // construct method for n component solution
75        Method_one method(n);
76
77        // inputs to OdeErrControl
78
79        double ti   = 0.;
80        double tf   = .9;
81        double smin = 1e-2;
82        double smax = 1.;
83        double scur = .5;
84        double erel = 1e-7;
85
86        CppAD::vector<double> xi(n);
87        CppAD::vector<double> eabs(n);
88        size_t i;
89        for(i = 0; i < n; i++)
90        {       xi[i]   = 0.;
91                eabs[i] = 0.;
92        }
93
94        // outputs from OdeErrControl
95
96        CppAD::vector<double> ef(n);
97        CppAD::vector<double> xf(n);
98       
99        xf = OdeErrControl(method,
100                ti, tf, xi, smin, smax, scur, eabs, erel, ef);
101
102        double check = 1.;
103        for(i = 0; i < n; i++)
104        {       check *= tf;
105                ok &= CppAD::NearEqual(check, xf[i], erel, 0.);
106        }
107
108        return ok;
109}
110
111/*
112Old example new just a test
113Define
114$latex X : \R \rightarrow \R^2$$ by
115$latex \[
116\begin{array}{rcl}
117X_0 (t) & = &  - \exp ( - w_0 t )  \\
118X_1 (t) & = & \frac{w_0}{w_1 - w_0} [ \exp ( - w_0 t ) - \exp( - w_1 t )]
119\end{array}
120\] $$
121It follows that $latex X_0 (0) = 1$$, $latex X_1 (0) = 0$$ and
122$latex \[
123\begin{array}{rcl}
124        X_0^{(1)} (t) & = & - w_0 X_0 (t)  \\
125        X_1^{(1)} (t) & = & + w_0 X_0 (t) - w_1 X_1 (t)
126\end{array}
127\] $$
128*/
129
130namespace {
131        // --------------------------------------------------------------
132        class Fun_two {
133        private:
134                 CppAD::vector<double> w;
135        public:
136                // constructor
137                Fun_two(const CppAD::vector<double> &w_) : w(w_)
138                { } 
139
140                // set f = x'(t)
141                void Ode(
142                        const double                &t, 
143                        const CppAD::vector<double> &x, 
144                        CppAD::vector<double>       &f)
145                {       f[0] = - w[0] * x[0];
146                        f[1] = + w[0] * x[0] - w[1] * x[1];     
147                }
148        };
149
150        // --------------------------------------------------------------
151        class Method_two {
152        private:
153                Fun_two F;
154        public:
155                // constructor
156                Method_two(const CppAD::vector<double> &w_) : F(w_)
157                { }
158                void step(
159                        double ta, 
160                        double tb, 
161                        CppAD::vector<double> &xa ,
162                        CppAD::vector<double> &xb ,
163                        CppAD::vector<double> &eb )
164                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
165                }
166                size_t order(void)
167                {       return 4; }
168        };
169}
170
171bool OdeErrControl_two(void)
172{       bool ok = true;     // initial return value
173
174        CppAD::vector<double> w(2);
175        w[0] = 10.;
176        w[1] = 1.;
177        Method_two method(w);
178
179        CppAD::vector<double> xi(2);
180        xi[0] = 1.;
181        xi[1] = 0.;
182
183        CppAD::vector<double> eabs(2);
184        eabs[0] = 1e-4;
185        eabs[1] = 1e-4;
186
187        // inputs
188        double ti   = 0.;
189        double tf   = 1.;
190        double smin = 1e-4;
191        double smax = 1.;
192        double scur = .5;
193        double erel = 0.;
194
195        // outputs
196        CppAD::vector<double> ef(2);
197        CppAD::vector<double> xf(2);
198        CppAD::vector<double> maxabs(2);
199        size_t nstep;
200
201       
202        xf = OdeErrControl(method,
203                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
204
205        double x0 = exp(-w[0]*tf);
206        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
207        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
208
209        double x1 = w[0] * (exp(-w[0]*tf) - exp(-w[1]*tf))/(w[1] - w[0]);
210        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
211        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
212
213        return ok;
214}
215
216bool OdeErrControl(void)
217{       bool ok = true;
218        ok     &= OdeErrControl_one();
219        ok     &= OdeErrControl_two();
220        return ok;
221}
222
Note: See TracBrowser for help on using the repository browser.