-
Notifications
You must be signed in to change notification settings - Fork 1.2k
/
scalar_initial_value_problem_test.cc
231 lines (199 loc) · 8.89 KB
/
scalar_initial_value_problem_test.cc
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
#include "drake/systems/analysis/scalar_initial_value_problem.h"
#include <gtest/gtest.h>
#include "drake/common/test_utilities/expect_throws_message.h"
#include "drake/common/unused.h"
#include "drake/systems/analysis/integrator_base.h"
#include "drake/systems/analysis/runge_kutta2_integrator.h"
namespace drake {
namespace systems {
namespace analysis {
namespace {
// Checks scalar IVP solver usage with multiple integrators.
GTEST_TEST(ScalarInitialValueProblemTest, UsingMultipleIntegrators) {
// Accuracy upper bound, as not all the integrators used below support
// error control.
const double kAccuracy = 1e-2;
// The initial time t₀, for IVP definition.
const double kInitialTime = 0.0;
// The initial state x₀, for IVP definition.
const double kInitialState = 0.0;
// The default parameters 𝐤₀, for IVP definition.
const VectorX<double> kParameters = VectorX<double>::Constant(1, 1.0);
// Instantiates a generic IVP for test purposes only,
// using a generic ODE dx/dt = -x + k₁, that does not
// model (nor attempts to model) any physical process.
ScalarInitialValueProblem<double> ivp(
[](const double& t, const double& x,
const VectorX<double>& k) -> double {
unused(t);
return -x + k[0];
}, kInitialState, kParameters);
// Testing against closed form solution of above's scalar IVP, which
// can be written as x(t; [k₁]) = k₁ + (x₀ - k₁) * e^(-(t - t₀)).
const double t0 = kInitialTime;
const double x0 = kInitialState;
const VectorX<double>& k1 = kParameters;
const double t1 = kInitialTime + 1.0;
EXPECT_NEAR(ivp.Solve(t0, t1), k1[0] + (x0 - k1[0]) * std::exp(-(t1 - t0)),
kAccuracy);
// Replaces default integrator.
const double kMaximumStep = 0.1;
const IntegratorBase<double>& default_integrator = ivp.get_integrator();
IntegratorBase<double>* configured_integrator =
ivp.reset_integrator<RungeKutta2Integrator<double>>(kMaximumStep);
EXPECT_NE(configured_integrator, &default_integrator);
EXPECT_EQ(configured_integrator, &ivp.get_integrator());
// Specifies a different parameter vector, but leaves both
// initial time and state as defaults.
const Eigen::Vector2d k2{1, 5.0};
ScalarInitialValueProblem<double> ivp2(
[](const double& t, const double& x,
const VectorX<double>& k) -> double {
unused(t);
return -x + k[0];
}, kInitialState, k2);
const double t2 = kInitialTime + 0.3;
// Testing against closed form solution of above's scalar IVP,
// which can be written as x(t; [k₁]) = k₁ + (x₀ - k₁) * e^(-(t - t₀)).
EXPECT_NEAR(ivp2.Solve(t0, t2), k2[0] + (x0 - k2[0]) * std::exp(-(t2 - t0)),
kAccuracy);
}
// Parameterized fixture for testing accuracy of scalar IVP solutions.
class ScalarInitialValueProblemAccuracyTest
: public ::testing::TestWithParam<double> {
protected:
void SetUp() {
integration_accuracy_ = GetParam();
}
// Expected accuracy for numerical integral
// evaluation in the relative tolerance sense.
double integration_accuracy_{0.};
};
// Accuracy test of the solution for the stored charge Q in an RC
// series circuit excited by a sinusoidal voltage source E(t),
// where dQ/dt = (E(t) - Q / Cs) / Rs and Q(t₀; [Rs, Cs]) = Q₀.
TEST_P(ScalarInitialValueProblemAccuracyTest, StoredCharge) {
// The initial time t₀.
const double kInitialTime = 0.0;
// The initial stored charge Q₀ at time t₀.
const double kInitialStoredCharge = 0.0;
// The resistance Rs of the resistor and the capacitance Cs
// of the capacitor.
const double kLowestResistance = 1.0;
const double kHighestResistance = 10.0;
const double kResistanceStep = 1.0;
const double kLowestCapacitance = 1.0;
const double kHighestCapacitance = 10.0;
const double kCapacitanceStep = 1.0;
const double kTotalTime = 1.0;
const double kTimeStep = 0.1;
const double Q0 = kInitialStoredCharge;
const double t0 = kInitialTime;
const double tf = kTotalTime;
for (double Rs = kLowestResistance; Rs <= kHighestResistance ;
Rs += kResistanceStep) {
for (double Cs = kLowestCapacitance; Cs <= kHighestCapacitance ;
Cs += kCapacitanceStep) {
const Eigen::Vector2d parameters{Rs, Cs};
// Instantiates the stored charge scalar IVP.
ScalarInitialValueProblem<double> stored_charge_ivp(
[](const double& t, const double& q,
const VectorX<double>& k) -> double {
const double Rs_ = k[0];
const double Cs_ = k[1];
return (std::sin(t) - q / Cs_) / Rs_;
}, kInitialStoredCharge, parameters);
IntegratorBase<double>& inner_integrator =
stored_charge_ivp.get_mutable_integrator();
inner_integrator.set_target_accuracy(integration_accuracy_);
const std::unique_ptr<ScalarDenseOutput<double>> stored_charge_approx =
stored_charge_ivp.DenseSolve(t0, tf);
const double tau = Rs * Cs;
const double tau_sq = tau * tau;
for (double t = kInitialTime; t <= kTotalTime; t += kTimeStep) {
// Tests are performed against the closed form
// solution for the scalar IVP described above, which is
// Q(t; [Rs, Cs]) = 1/Rs * (τ²/ (1 + τ²) * e^(-t / τ) +
// τ / √(1 + τ²) * sin(t - arctan(τ)))
// where τ = Rs * Cs for Q(t₀ = 0; [Rs, Cs]) = Q₀ = 0.
const double solution = (
tau_sq / (1. + tau_sq) * std::exp(-t / tau)
+ tau / std::sqrt(1. + tau_sq)
* std::sin(t - std::atan(tau))) / Rs;
EXPECT_NEAR(stored_charge_ivp.Solve(t0, t),
solution, integration_accuracy_)
<< "Failure solving dQ/dt = (sin(t) - Q / Cs) / Rs using Q(t₀ = "
<< t0 << "; [Rs, Cs]) = " << Q0 << " for t = " << t << ", Rs = "
<< Rs << " and Cs = " << Cs << " to an accuracy of "
<< integration_accuracy_;
EXPECT_NEAR(stored_charge_approx->EvaluateScalar(t),
solution, integration_accuracy_)
<< "Failure approximating the solution for"
<< " dQ/dt = (sin(t) - Q / Cs) / Rs using Q(t₀ = "
<< t0 << "; [Rs, Cs]) = " << Q0 << " for t = " << t
<< ", Rs = " << Rs << " and Cs = " << Cs
<< " to an accuracy of " << integration_accuracy_
<< " with solver's continuous extension.";
}
}
}
}
// Accuracy test of the solution for population growth N, described
// by dN/dt = r * N and N(t₀; r) = N₀.
TEST_P(ScalarInitialValueProblemAccuracyTest, PopulationGrowth) {
// The initial time t₀.
const double kInitialTime = 0.0;
// The initial population N₀ at time t₀.
const double kInitialPopulation = 10.0;
// The Malthusian parameter r that shapes the growth of the
// population.
const double kLowestMalthusParam = 0.1;
const double kHighestMalthusParam = 1.0;
const double kMalthusParamStep = 0.1;
const double kTotalTime = 1.0;
const double kTimeStep = 0.1;
const double N0 = kInitialPopulation;
const double t0 = kInitialTime;
const double tf = kTotalTime;
for (double r = kLowestMalthusParam; r <= kHighestMalthusParam;
r += kMalthusParamStep) {
ScalarInitialValueProblem<double> population_growth_ivp(
[](const double& t, const double& n,
const VectorX<double>& k) -> double {
unused(t);
return k[0] * n;
},
kInitialPopulation, Vector1<double>{r});
IntegratorBase<double>& inner_integrator =
population_growth_ivp.get_mutable_integrator();
inner_integrator.set_target_accuracy(integration_accuracy_);
const std::unique_ptr<ScalarDenseOutput<double>> population_growth_approx =
population_growth_ivp.DenseSolve(t0, tf);
for (double t = kInitialTime; t <= kTotalTime; t += kTimeStep) {
// Tests are performed against the closed form
// solution for the IVP described above, which is
// N(t; r) = N₀ * e^(r * t).
const double solution = N0 * std::exp(r * t);
EXPECT_NEAR(population_growth_ivp.Solve(t0, t),
solution, integration_accuracy_)
<< "Failure solving dN/dt = r * N using N(t₀ = "
<< t0 << "; r) = " << N0 << " for t = " << t
<< " and r = " << r << " to an accuracy of "
<< integration_accuracy_;
EXPECT_NEAR(population_growth_approx->EvaluateScalar(t),
solution, integration_accuracy_)
<< "Failure approximating the solution for dN/dt = r * N"
<< " using N(t₀ = " << t0 << "; r) = " << N0 << " for t = "
<< t << " and r = " << r << " to an accuracy of "
<< integration_accuracy_ << " with solver's continuous extension.";
}
}
}
INSTANTIATE_TEST_SUITE_P(
IncreasingAccuracyScalarInitialValueProblemTests,
ScalarInitialValueProblemAccuracyTest,
::testing::Values(1e-1, 1e-2, 1e-3, 1e-4, 1e-5));
} // namespace
} // namespace analysis
} // namespace systems
} // namespace drake