1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin ode_err_maxabs.cpp}
OdeErrControl: Example and Test Using Maxabs Argument
#####################################################
Define
:math:`X : \B{R} \rightarrow \B{R}^2` by
.. math::
:nowrap:
\begin{eqnarray}
X_0 (t) & = & - \exp ( - w_0 t ) \\
X_1 (t) & = & \frac{w_0}{w_1 - w_0} [ \exp ( - w_0 t ) - \exp( - w_1 t )]
\end{eqnarray}
It follows that :math:`X_0 (0) = 1`, :math:`X_1 (0) = 0` and
.. math::
:nowrap:
\begin{eqnarray}
X_0^{(1)} (t) & = & - w_0 X_0 (t) \\
X_1^{(1)} (t) & = & + w_0 X_0 (t) - w_1 X_1 (t)
\end{eqnarray}
Note that :math:`X_1 (0)` is zero an if :math:`w_0 t` is large,
:math:`X_0 (t)` is near zero.
This example tests OdeErrControl using the *maxabs* argument.
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end ode_err_maxabs.cpp}
*/
// BEGIN C++
# include <cstddef> // for size_t
# include <cmath> // for exp
# include <cppad/utility/ode_err_control.hpp> // CppAD::OdeErrControl
# include <cppad/utility/near_equal.hpp> // CppAD::NearEqual
# include <cppad/utility/vector.hpp> // CppAD::vector
# include <cppad/utility/runge_45.hpp> // CppAD::Runge45
namespace {
// --------------------------------------------------------------
class Fun {
private:
CppAD::vector<double> w;
public:
// constructor
Fun(const CppAD::vector<double> &w_) : w(w_)
{ }
// set f = x'(t)
void Ode(
const double &t,
const CppAD::vector<double> &x,
CppAD::vector<double> &f)
{ f[0] = - w[0] * x[0];
f[1] = + w[0] * x[0] - w[1] * x[1];
}
};
// --------------------------------------------------------------
class Method {
private:
Fun F;
public:
// constructor
Method(const CppAD::vector<double> &w_) : F(w_)
{ }
void step(
double ta,
double tb,
CppAD::vector<double> &xa ,
CppAD::vector<double> &xb ,
CppAD::vector<double> &eb )
{ xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
}
size_t order(void)
{ return 4; }
};
}
bool OdeErrMaxabs(void)
{ bool ok = true; // initial return value
CppAD::vector<double> w(2);
w[0] = 10.;
w[1] = 1.;
Method method(w);
CppAD::vector<double> xi(2);
xi[0] = 1.;
xi[1] = 0.;
CppAD::vector<double> eabs(2);
eabs[0] = 0.;
eabs[1] = 0.;
CppAD::vector<double> ef(2);
CppAD::vector<double> xf(2);
CppAD::vector<double> maxabs(2);
double ti = 0.;
double tf = 1.;
double smin = .5;
double smax = 1.;
double scur = .5;
double erel = 1e-4;
bool accurate = false;
while( ! accurate )
{ xf = OdeErrControl(method,
ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs);
accurate = true;
size_t i;
for(i = 0; i < 2; i++)
accurate &= ef[i] <= erel * maxabs[i];
if( ! accurate )
smin = smin / 2;
}
double x0 = exp(-w[0]*tf);
ok &= CppAD::NearEqual(x0, xf[0], erel, 0.);
ok &= CppAD::NearEqual(0., ef[0], erel, erel);
double x1 = w[0] * (exp(-w[0]*tf) - exp(-w[1]*tf))/(w[1] - w[0]);
ok &= CppAD::NearEqual(x1, xf[1], erel, 0.);
ok &= CppAD::NearEqual(0., ef[1], erel, erel);
return ok;
}
// END C++
|
