source: trunk/test_more/ode_err_control.cpp @ 990

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

trunk: Fix infinite loop in OdeErrControl?.

whats_new_07.omh: user's view of the changes.
ode_err_control.cpp: add test case for the infinite loop.
config.h: update version for subversion distribution option.
ode_err_control.hpp: fix infinite loop here.

File size: 10.4 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/*
216Define
217$latex X : \R \rightarrow \R^2$$ by
218$latex \[
219\begin{array}{rcl}
220        X_0 (0)       & = & 1  \\
221        X_1 (0)       & = & 0  \\
222        X_0^{(1)} (t) & = & 2 \alpha t X_0 (t)  \\
223        X_1^{(1)} (t) & = &  1 / X_0 (t)
224\end{array}
225\] $$
226It follows that
227$latex \[
228\begin{array}{rcl}
229X_0 (t) & = &  \exp ( \alpha t^2 )  \\
230X_1 (t) & = &  \int_0^t \exp( - \alpha s^2 ) {\bf d} s \\
231& = & 
232\frac{1}{ \sqrt{\alpha} \int_0^{\sqrt{\alpha} t} \exp( - r^2 ) {\bf d} r
233\\
234& = & \frac{\sqrt{\pi}}{ 2 \sqrt{\alpha} {\rm erf} ( \sqrt{\alpha} t )
235\end{array}
236\] $$
237If $latex X_0 (t) < 0$$,
238we return $code nan$$ in order to inform
239$code OdeErrControl$$ that its is taking to large a step.
240
241*/
242
243# include <cppad/rosen_34.hpp>          // CppAD::Rosen34
244# include <cppad/cppad.hpp>
245
246namespace {
247        // --------------------------------------------------------------
248        class Fun_three {
249        private:
250                const double alpha_;
251                bool was_negative_;
252        public:
253                // constructor
254                Fun_three(double alpha) : alpha_(alpha), was_negative_(false)
255                { } 
256
257                // set f = x'(t)
258                void Ode(
259                        const double                &t, 
260                        const CppAD::vector<double> &x, 
261                        CppAD::vector<double>       &f)
262                {       f[0] = 2. * alpha_ * t * x[0];
263                        f[1] = 1. / x[0];       
264                        // case where ODE does not make sense
265                        if( x[0] < 0. || x[1] < 0. )
266                        {       was_negative_ = true;
267                                f[0] = CppAD::nan(0.);
268                        }
269                }
270                // set f_t = df / dt
271                void Ode_ind(
272                        const double                &t, 
273                        const CppAD::vector<double> &x, 
274                        CppAD::vector<double>       &f_t)
275                {
276                        f_t[0] =  2. * alpha_ * x[0];
277                        f_t[1] = 0.;   
278                        if( x[0] < 0. || x[1] < 0. )
279                        {       was_negative_ = true;
280                                f_t[0] = CppAD::nan(0.);
281                        }
282                }
283                // set f_x = df / dx
284                void Ode_dep(
285                        const double                &t, 
286                        const CppAD::vector<double> &x, 
287                        CppAD::vector<double>       &f_x)
288                {       double x0_sq = x[0] * x[0];
289                        f_x[0 * 2 + 0] = 2. * alpha_ * t;   // f0 w.r.t. x0
290                        f_x[0 * 2 + 1] = 0.;                // f0 w.r.t. x1
291                        f_x[1 * 2 + 0] = -1./x0_sq;         // f1 w.r.t. x0
292                        f_x[1 * 2 + 1] = 0.;                // f1 w.r.t. x1
293                        if( x[0] < 0. || x[1] < 0. )
294                        {       was_negative_ = true;
295                                f_x[0] = CppAD::nan(0.);
296                        }
297                }
298                bool was_negative(void)
299                {       return was_negative_; }
300
301        };
302
303        // --------------------------------------------------------------
304        class Method_three {
305        public:
306                Fun_three F;
307
308                // constructor
309                Method_three(double alpha) : F(alpha)
310                { }
311                void step(
312                        double ta, 
313                        double tb, 
314                        CppAD::vector<double> &xa ,
315                        CppAD::vector<double> &xb ,
316                        CppAD::vector<double> &eb )
317                {       xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb);
318                }
319                size_t order(void)
320                {       return 3; }
321        };
322}
323
324bool OdeErrControl_three(void)
325{       bool ok = true;     // initial return value
326
327        double alpha = 10.;
328        Method_three method(alpha);
329
330        CppAD::vector<double> xi(2);
331        xi[0] = 1.;
332        xi[1] = 0.;
333
334        CppAD::vector<double> eabs(2);
335        eabs[0] = 1e-4;
336        eabs[1] = 1e-4;
337
338        // inputs
339        double ti   = 0.;
340        double tf   = 1.;
341        double smin = 1e-4;
342        double smax = 1.;
343        double scur = 1.;
344        double erel = 0.;
345
346        // outputs
347        CppAD::vector<double> ef(2);
348        CppAD::vector<double> xf(2);
349        CppAD::vector<double> maxabs(2);
350        size_t nstep;
351
352       
353        xf = OdeErrControl(method,
354                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
355
356
357        double x0       = exp( alpha * tf * tf );
358        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
359        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
360
361        double root_pi    = sqrt( 4. * atan(1.));
362        double root_alpha = sqrt( alpha );
363        double x1 = CppAD::erf(alpha * tf) * root_pi / (2 * root_alpha);
364        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
365        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
366
367        ok &= method.F.was_negative();
368
369        return ok;
370}
371
372namespace {
373        // --------------------------------------------------------------
374        class Fun_four {
375        private:
376                 size_t n;   // dimension of the state space
377        public:
378                // constructor
379                Fun_four(size_t n_) : n(n_)
380                { } 
381
382                // given x(0) = 0
383                // solution is x_i (t) = t^(i+1)
384                void Ode(
385                        const double                &t, 
386                        const CppAD::vector<double> &x, 
387                        CppAD::vector<double>       &f)
388                {       size_t i;
389                        f[0] = CppAD::nan(0.);
390                        for(i = 1; i < n; i++)
391                                f[i] = (i+1) * x[i-1];
392                }
393        };
394
395        // --------------------------------------------------------------
396        class Method_four {
397        private:
398                Fun_four F;
399        public:
400                // constructor
401                Method_four(size_t n_) : F(n_)
402                { }
403                void step(
404                        double ta, 
405                        double tb, 
406                        CppAD::vector<double> &xa ,
407                        CppAD::vector<double> &xb ,
408                        CppAD::vector<double> &eb )
409                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
410                }
411                size_t order(void)
412                {       return 4; }
413        };
414}
415
416bool OdeErrControl_four(void)
417{       bool   ok = true;     // initial return value
418       
419        // construct method for n component solution
420        size_t  n = 6;       
421        Method_four method(n);
422
423        // inputs to OdeErrControl
424
425        // special case where scur is converted to ti - tf
426        // (so it is not equal to smin)
427        double ti   = 0.;
428        double tf   = .9;
429        double smin = .8;
430        double smax = 1.;
431        double scur = smin;
432        double erel = 1e-7;
433
434        CppAD::vector<double> xi(n);
435        CppAD::vector<double> eabs(n);
436        size_t i;
437        for(i = 0; i < n; i++)
438        {       xi[i]   = 0.;
439                eabs[i] = 0.;
440        }
441
442        // outputs from OdeErrControl
443        CppAD::vector<double> ef(n);
444        CppAD::vector<double> xf(n);
445       
446        xf = OdeErrControl(method,
447                ti, tf, xi, smin, smax, scur, eabs, erel, ef);
448
449        // check that Fun_four always returning nan results in nan
450        for(i = 0; i < n; i++)
451        {       ok &= CppAD::isnan(xf[i]);
452                ok &= CppAD::isnan(ef[i]);
453        }
454
455        return ok;
456}
457
458// ==========================================================================
459bool OdeErrControl(void)
460{       bool ok = true;
461        ok     &= OdeErrControl_one();
462        ok     &= OdeErrControl_two();
463        ok     &= OdeErrControl_three();
464        ok     &= OdeErrControl_four();
465        return ok;
466}
467
Note: See TracBrowser for help on using the repository browser.