source: trunk/test_more/ode_err_control.cpp @ 1368

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

trunk: use svn_add_id.sh to add Id keyword and property to most files.

File size: 10.4 KB
Line 
1/* $Id$ */
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# include <cstddef>                    // for size_t
14# include <cmath>                      // for exp
15# include <cppad/ode_err_control.hpp>  // CppAD::OdeErrControl
16# include <cppad/near_equal.hpp>       // CppAD::NearEqual
17# include <cppad/vector.hpp>           // CppAD::vector
18# include <cppad/runge_45.hpp>         // CppAD::Runge45
19
20/* -------------------------------------------------------------------------
21Test relative error with zero initial conditions.
22(Uses minimum step size to integrate).
23*/
24
25namespace {
26        // --------------------------------------------------------------
27        class Fun_one {
28        private:
29                 size_t n;   // dimension of the state space
30        public:
31                // constructor
32                Fun_one(size_t n_) : n(n_)
33                { } 
34
35                // given x(0) = 0
36                // solution is x_i (t) = t^(i+1)
37                void Ode(
38                        const double                &t, 
39                        const CppAD::vector<double> &x, 
40                        CppAD::vector<double>       &f)
41                {       size_t i;
42                        f[0] = 1.;
43                        for(i = 1; i < n; i++)
44                                f[i] = (i+1) * x[i-1];
45                }
46        };
47
48        // --------------------------------------------------------------
49        class Method_one {
50        private:
51                Fun_one F;
52        public:
53                // constructor
54                Method_one(size_t n_) : F(n_)
55                { }
56                void step(
57                        double ta, 
58                        double tb, 
59                        CppAD::vector<double> &xa ,
60                        CppAD::vector<double> &xb ,
61                        CppAD::vector<double> &eb )
62                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
63                }
64                size_t order(void)
65                {       return 4; }
66        };
67}
68
69bool OdeErrControl_one(void)
70{       bool   ok = true;     // initial return value
71
72        // Runge45 should yield exact results for x_i (t) = t^(i+1), i < 4
73        size_t  n = 6;       
74       
75        // construct method for n component solution
76        Method_one method(n);
77
78        // inputs to OdeErrControl
79
80        double ti   = 0.;
81        double tf   = .9;
82        double smin = 1e-2;
83        double smax = 1.;
84        double scur = .5;
85        double erel = 1e-7;
86
87        CppAD::vector<double> xi(n);
88        CppAD::vector<double> eabs(n);
89        size_t i;
90        for(i = 0; i < n; i++)
91        {       xi[i]   = 0.;
92                eabs[i] = 0.;
93        }
94
95        // outputs from OdeErrControl
96
97        CppAD::vector<double> ef(n);
98        CppAD::vector<double> xf(n);
99       
100        xf = OdeErrControl(method,
101                ti, tf, xi, smin, smax, scur, eabs, erel, ef);
102
103        double check = 1.;
104        for(i = 0; i < n; i++)
105        {       check *= tf;
106                ok &= CppAD::NearEqual(check, xf[i], erel, 0.);
107        }
108
109        return ok;
110}
111
112/*
113Old example now just a test
114Define
115$latex X : \R \rightarrow \R^2$$ by
116$latex \[
117\begin{array}{rcl}
118X_0 (t) & = &  - \exp ( - w_0 t )  \\
119X_1 (t) & = & \frac{w_0}{w_1 - w_0} [ \exp ( - w_0 t ) - \exp( - w_1 t )]
120\end{array}
121\] $$
122It follows that $latex X_0 (0) = 1$$, $latex X_1 (0) = 0$$ and
123$latex \[
124\begin{array}{rcl}
125        X_0^{(1)} (t) & = & - w_0 X_0 (t)  \\
126        X_1^{(1)} (t) & = & + w_0 X_0 (t) - w_1 X_1 (t)
127\end{array}
128\] $$
129*/
130
131namespace {
132        // --------------------------------------------------------------
133        class Fun_two {
134        private:
135                 CppAD::vector<double> w;
136        public:
137                // constructor
138                Fun_two(const CppAD::vector<double> &w_) : w(w_)
139                { } 
140
141                // set f = x'(t)
142                void Ode(
143                        const double                &t, 
144                        const CppAD::vector<double> &x, 
145                        CppAD::vector<double>       &f)
146                {       f[0] = - w[0] * x[0];
147                        f[1] = + w[0] * x[0] - w[1] * x[1];     
148                }
149        };
150
151        // --------------------------------------------------------------
152        class Method_two {
153        private:
154                Fun_two F;
155        public:
156                // constructor
157                Method_two(const CppAD::vector<double> &w_) : F(w_)
158                { }
159                void step(
160                        double ta, 
161                        double tb, 
162                        CppAD::vector<double> &xa ,
163                        CppAD::vector<double> &xb ,
164                        CppAD::vector<double> &eb )
165                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
166                }
167                size_t order(void)
168                {       return 4; }
169        };
170}
171
172bool OdeErrControl_two(void)
173{       bool ok = true;     // initial return value
174
175        CppAD::vector<double> w(2);
176        w[0] = 10.;
177        w[1] = 1.;
178        Method_two method(w);
179
180        CppAD::vector<double> xi(2);
181        xi[0] = 1.;
182        xi[1] = 0.;
183
184        CppAD::vector<double> eabs(2);
185        eabs[0] = 1e-4;
186        eabs[1] = 1e-4;
187
188        // inputs
189        double ti   = 0.;
190        double tf   = 1.;
191        double smin = 1e-4;
192        double smax = 1.;
193        double scur = .5;
194        double erel = 0.;
195
196        // outputs
197        CppAD::vector<double> ef(2);
198        CppAD::vector<double> xf(2);
199        CppAD::vector<double> maxabs(2);
200        size_t nstep;
201
202       
203        xf = OdeErrControl(method,
204                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
205
206        double x0 = exp(-w[0]*tf);
207        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
208        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
209
210        double x1 = w[0] * (exp(-w[0]*tf) - exp(-w[1]*tf))/(w[1] - w[0]);
211        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
212        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
213
214        return ok;
215}
216/*
217Define
218$latex X : \R \rightarrow \R^2$$ by
219$latex \[
220\begin{array}{rcl}
221        X_0 (0)       & = & 1  \\
222        X_1 (0)       & = & 0  \\
223        X_0^{(1)} (t) & = & 2 \alpha t X_0 (t)  \\
224        X_1^{(1)} (t) & = &  1 / X_0 (t)
225\end{array}
226\] $$
227It follows that
228$latex \[
229\begin{array}{rcl}
230X_0 (t) & = &  \exp ( \alpha t^2 )  \\
231X_1 (t) & = &  \int_0^t \exp( - \alpha s^2 ) {\bf d} s \\
232& = & 
233\frac{1}{ \sqrt{\alpha} \int_0^{\sqrt{\alpha} t} \exp( - r^2 ) {\bf d} r
234\\
235& = & \frac{\sqrt{\pi}}{ 2 \sqrt{\alpha} {\rm erf} ( \sqrt{\alpha} t )
236\end{array}
237\] $$
238If $latex X_0 (t) < 0$$,
239we return $code nan$$ in order to inform
240$code OdeErrControl$$ that its is taking to large a step.
241
242*/
243
244# include <cppad/rosen_34.hpp>          // CppAD::Rosen34
245# include <cppad/cppad.hpp>
246
247namespace {
248        // --------------------------------------------------------------
249        class Fun_three {
250        private:
251                const double alpha_;
252                bool was_negative_;
253        public:
254                // constructor
255                Fun_three(double alpha) : alpha_(alpha), was_negative_(false)
256                { } 
257
258                // set f = x'(t)
259                void Ode(
260                        const double                &t, 
261                        const CppAD::vector<double> &x, 
262                        CppAD::vector<double>       &f)
263                {       f[0] = 2. * alpha_ * t * x[0];
264                        f[1] = 1. / x[0];       
265                        // case where ODE does not make sense
266                        if( x[0] < 0. || x[1] < 0. )
267                        {       was_negative_ = true;
268                                f[0] = CppAD::nan(0.);
269                        }
270                }
271                // set f_t = df / dt
272                void Ode_ind(
273                        const double                &t, 
274                        const CppAD::vector<double> &x, 
275                        CppAD::vector<double>       &f_t)
276                {
277                        f_t[0] =  2. * alpha_ * x[0];
278                        f_t[1] = 0.;   
279                        if( x[0] < 0. || x[1] < 0. )
280                        {       was_negative_ = true;
281                                f_t[0] = CppAD::nan(0.);
282                        }
283                }
284                // set f_x = df / dx
285                void Ode_dep(
286                        const double                &t, 
287                        const CppAD::vector<double> &x, 
288                        CppAD::vector<double>       &f_x)
289                {       double x0_sq = x[0] * x[0];
290                        f_x[0 * 2 + 0] = 2. * alpha_ * t;   // f0 w.r.t. x0
291                        f_x[0 * 2 + 1] = 0.;                // f0 w.r.t. x1
292                        f_x[1 * 2 + 0] = -1./x0_sq;         // f1 w.r.t. x0
293                        f_x[1 * 2 + 1] = 0.;                // f1 w.r.t. x1
294                        if( x[0] < 0. || x[1] < 0. )
295                        {       was_negative_ = true;
296                                f_x[0] = CppAD::nan(0.);
297                        }
298                }
299                bool was_negative(void)
300                {       return was_negative_; }
301
302        };
303
304        // --------------------------------------------------------------
305        class Method_three {
306        public:
307                Fun_three F;
308
309                // constructor
310                Method_three(double alpha) : F(alpha)
311                { }
312                void step(
313                        double ta, 
314                        double tb, 
315                        CppAD::vector<double> &xa ,
316                        CppAD::vector<double> &xb ,
317                        CppAD::vector<double> &eb )
318                {       xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb);
319                }
320                size_t order(void)
321                {       return 3; }
322        };
323}
324
325bool OdeErrControl_three(void)
326{       bool ok = true;     // initial return value
327
328        double alpha = 10.;
329        Method_three method(alpha);
330
331        CppAD::vector<double> xi(2);
332        xi[0] = 1.;
333        xi[1] = 0.;
334
335        CppAD::vector<double> eabs(2);
336        eabs[0] = 1e-4;
337        eabs[1] = 1e-4;
338
339        // inputs
340        double ti   = 0.;
341        double tf   = 1.;
342        double smin = 1e-4;
343        double smax = 1.;
344        double scur = 1.;
345        double erel = 0.;
346
347        // outputs
348        CppAD::vector<double> ef(2);
349        CppAD::vector<double> xf(2);
350        CppAD::vector<double> maxabs(2);
351        size_t nstep;
352
353       
354        xf = OdeErrControl(method,
355                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
356
357
358        double x0       = exp( alpha * tf * tf );
359        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
360        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
361
362        double root_pi    = sqrt( 4. * atan(1.));
363        double root_alpha = sqrt( alpha );
364        double x1 = CppAD::erf(alpha * tf) * root_pi / (2 * root_alpha);
365        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
366        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
367
368        ok &= method.F.was_negative();
369
370        return ok;
371}
372
373namespace {
374        // --------------------------------------------------------------
375        class Fun_four {
376        private:
377                 size_t n;   // dimension of the state space
378        public:
379                // constructor
380                Fun_four(size_t n_) : n(n_)
381                { } 
382
383                // given x(0) = 0
384                // solution is x_i (t) = t^(i+1)
385                void Ode(
386                        const double                &t, 
387                        const CppAD::vector<double> &x, 
388                        CppAD::vector<double>       &f)
389                {       size_t i;
390                        f[0] = CppAD::nan(0.);
391                        for(i = 1; i < n; i++)
392                                f[i] = (i+1) * x[i-1];
393                }
394        };
395
396        // --------------------------------------------------------------
397        class Method_four {
398        private:
399                Fun_four F;
400        public:
401                // constructor
402                Method_four(size_t n_) : F(n_)
403                { }
404                void step(
405                        double ta, 
406                        double tb, 
407                        CppAD::vector<double> &xa ,
408                        CppAD::vector<double> &xb ,
409                        CppAD::vector<double> &eb )
410                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
411                }
412                size_t order(void)
413                {       return 4; }
414        };
415}
416
417bool OdeErrControl_four(void)
418{       bool   ok = true;     // initial return value
419       
420        // construct method for n component solution
421        size_t  n = 6;       
422        Method_four method(n);
423
424        // inputs to OdeErrControl
425
426        // special case where scur is converted to ti - tf
427        // (so it is not equal to smin)
428        double ti   = 0.;
429        double tf   = .9;
430        double smin = .8;
431        double smax = 1.;
432        double scur = smin;
433        double erel = 1e-7;
434
435        CppAD::vector<double> xi(n);
436        CppAD::vector<double> eabs(n);
437        size_t i;
438        for(i = 0; i < n; i++)
439        {       xi[i]   = 0.;
440                eabs[i] = 0.;
441        }
442
443        // outputs from OdeErrControl
444        CppAD::vector<double> ef(n);
445        CppAD::vector<double> xf(n);
446       
447        xf = OdeErrControl(method,
448                ti, tf, xi, smin, smax, scur, eabs, erel, ef);
449
450        // check that Fun_four always returning nan results in nan
451        for(i = 0; i < n; i++)
452        {       ok &= CppAD::isnan(xf[i]);
453                ok &= CppAD::isnan(ef[i]);
454        }
455
456        return ok;
457}
458
459// ==========================================================================
460bool ode_err_control(void)
461{       bool ok = true;
462        ok     &= OdeErrControl_one();
463        ok     &= OdeErrControl_two();
464        ok     &= OdeErrControl_three();
465        ok     &= OdeErrControl_four();
466        return ok;
467}
468
Note: See TracBrowser for help on using the repository browser.