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
145
146
147
148
149
150
151
152
153
154
155
156
|
// 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_gear_control.cpp}
OdeGearControl: Example and Test
################################
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}
The example tests OdeGearControl using the relations above:
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end ode_gear_control.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
# include <cppad/utility/ode_gear_control.hpp> // CppAD::OdeGearControl
namespace {
// --------------------------------------------------------------
class Fun {
private:
CPPAD_TESTVECTOR(double) w;
public:
// constructor
Fun(const CPPAD_TESTVECTOR(double) &w_) : w(w_)
{ }
// set f = x'(t)
template <class Scalar>
void Ode(
const Scalar &t,
const CPPAD_TESTVECTOR(Scalar) &x,
CPPAD_TESTVECTOR(Scalar) &f)
{ f[0] = - w[0] * x[0];
f[1] = + w[0] * x[0] - w[1] * x[1];
}
void Ode_dep(
const double &t,
const CPPAD_TESTVECTOR(double) &x,
CPPAD_TESTVECTOR(double) &f_x)
{ using namespace CppAD;
size_t n = x.size();
CPPAD_TESTVECTOR(AD<double>) T(1);
CPPAD_TESTVECTOR(AD<double>) X(n);
CPPAD_TESTVECTOR(AD<double>) F(n);
// set argument values
T[0] = t;
size_t i, j;
for(i = 0; i < n; i++)
X[i] = x[i];
// declare independent variables
Independent(X);
// compute f(t, x)
this->Ode(T[0], X, F);
// define AD function object
ADFun<double> fun(X, F);
// compute partial of f w.r.t x
CPPAD_TESTVECTOR(double) dx(n);
CPPAD_TESTVECTOR(double) df(n);
for(j = 0; j < n; j++)
dx[j] = 0.;
for(j = 0; j < n; j++)
{ dx[j] = 1.;
df = fun.Forward(1, dx);
for(i = 0; i < n; i++)
f_x [i * n + j] = df[i];
dx[j] = 0.;
}
}
};
}
bool OdeGearControl(void)
{ bool ok = true; // initial return value
using CppAD::NearEqual;
double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
CPPAD_TESTVECTOR(double) w(2);
w[0] = 10.;
w[1] = 1.;
Fun F(w);
CPPAD_TESTVECTOR(double) xi(2);
xi[0] = 1.;
xi[1] = 0.;
CPPAD_TESTVECTOR(double) eabs(2);
eabs[0] = 1e-4;
eabs[1] = 1e-4;
// return values
CPPAD_TESTVECTOR(double) ef(2);
CPPAD_TESTVECTOR(double) maxabs(2);
CPPAD_TESTVECTOR(double) xf(2);
size_t nstep;
// input values
size_t M = 5;
double ti = 0.;
double tf = 1.;
double smin = 1e-8;
double smax = 1.;
double sini = eps99;
double erel = 0.;
xf = CppAD::OdeGearControl(F, M,
ti, tf, xi, smin, smax, sini, eabs, erel, ef, maxabs, nstep);
double x0 = exp(-w[0]*tf);
ok &= NearEqual(x0, xf[0], 1e-4, 1e-4);
ok &= NearEqual(0., ef[0], 1e-4, 1e-4);
double x1 = w[0] * (exp(-w[0]*tf) - exp(-w[1]*tf))/(w[1] - w[0]);
ok &= NearEqual(x1, xf[1], 1e-4, 1e-4);
ok &= NearEqual(0., ef[1], 1e-4, 1e-4);
return ok;
}
// END C++
|
