Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
ode.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_ODE_H
7#define THEORETICA_ODE_H
8
9#include "../algebra/vec.h"
10
11
12namespace theoretica {
13
15 namespace ode {
16
17
22 template<typename Vector = vec<real>>
24
27
30
31
34
35
39 ode_solution_t(size_t steps, const Vector& x0, real t0) {
40
41 t.resize(steps);
42 x.resize(steps);
43
44 t[0] = t0;
45 x[0] = x0;
46 }
47
48
49#ifndef THEORETICA_NO_PRINT
50
52 inline std::string to_string(const std::string& separator = " ") const {
53
54 if (t.size() != x.size()) {
55 TH_MATH_ERROR("ode_solution_t::to_string", t.size(), INVALID_ARGUMENT);
56 return "";
57 }
58
59 std::stringstream res;
60
61 for (unsigned int i = 0; i < t.size(); ++i) {
62
63 res << t[i] << separator;
64
65 for (unsigned int j = 0; j < x[i].size(); ++j) {
66
67 res << x[i][j];
68
69 if(j != x[i].size() - 1)
70 res << separator;
71 else
72 res << "\n";
73 }
74 }
75
76 return res.str();
77 }
78
79
81 inline operator std::string() {
82 return to_string();
83 }
84
85
88 inline friend std::ostream& operator<<(
89 std::ostream& out, const ode_solution_t<Vector>& obj) {
90 return out << obj.to_string();
91 }
92#endif
93
94 };
95
96
99
100
103
104
107
108
111
112
115
116
122 template<typename Vector>
123 using ode_function = std::function<Vector(real, const Vector&)>;
124
125
126 // Steppers
127 // (functions which compute one iteration of a method)
128
129
139 template<typename Vector, typename OdeFunction = ode_function<Vector>>
140 inline Vector step_euler(OdeFunction f, const Vector& x, real t, real h = 0.0001) {
141
142 return x + h * f(t, x);
143 }
144
145
155 template<typename Vector, typename OdeFunction = ode_function<Vector>>
156 inline Vector step_midpoint(OdeFunction f, const Vector& x, real t, real h = 0.0001) {
157
158 return x + h * f(t + h / 2.0, x + f(t, x) * h / 2.0);
159 }
160
161
171 template<typename Vector, typename OdeFunction = ode_function<Vector>>
172 inline Vector step_heun(OdeFunction f, const Vector& x, real t, real h = 0.001) {
173
174 const Vector x_p = x + h * f(t, x);
175 const real t_new = t + h;
176
177 return x + (f(t, x) + f(t_new, x_p)) * (h / 2.0);
178 }
179
180
191 template<typename Vector, typename OdeFunction = ode_function<Vector>>
192 inline Vector step_rk2(OdeFunction f, const Vector& x, real t, real h = 0.001) {
193
194 const Vector k1 = f(t, x);
195 const Vector k2 = f(t + (h / 2.0), x + k1 * (h / 2.0));
196
197 return x + k2 * h;
198 }
199
200
211 template<typename Vector, typename OdeFunction = ode_function<Vector>>
212 inline Vector step_rk4(OdeFunction f, const Vector& x, real t, real h = 0.01) {
213
214 const real half = h / 2.0;
215
216 const Vector k1 = f(t, x);
217 const Vector k2 = f(t + half, x + k1 * half);
218 const Vector k3 = f(t + half, x + k2 * half);
219 const Vector k4 = f(t + h, x + k3 * h);
220
221 return x + (k1 + k2 * 2.0 + k3 * 2.0 + k4) * (h / 6.0);
222 }
223
224
235 template<typename Vector, typename OdeFunction = ode_function<Vector>>
236 inline Vector step_k38(OdeFunction f, const Vector& x, real t, real h = 0.001) {
237
238 const Vector k1 = f(t, x);
239 const Vector k2 = f(t + (h / 3.0), x + k1 * (h / 3.0));
240 const Vector k3 = f(t + (h * 2.0 / 3.0), x + h * (-k1 / 3.0 + k2));
241 const Vector k4 = f(t + h, x + h * (k1 - k2 + k1));
242
243 return x + (k1 + 3.0 * k2 + 3.0 * k3 + k4) * (h / 8.0);
244 }
245
246
257 template<typename Vector, typename OdeFunction = ode_function<Vector>>
259 OdeFunction f, const Vector& x0, real t0,
260 const Vector& x1, real t1, real h = 0.001) {
261
262 return x1 + h * ((3. / 2.) * f(t1, x1) - f(t0, x0) / 2.);
263 }
264
265
276 template<typename Vector, typename OdeFunction = ode_function<Vector>>
278 OdeFunction f, const Vector& x0, real t0, const Vector& x1, real t1,
279 const Vector& x2, real t2, real h = 0.001) {
280
281 return x2 + h * (
282 (23. / 12.) * f(t2, x2) - (4. / 3.) * f(t1, x1) + (3. / 8.) * f(t0, x0)
283 );
284 }
285
286
287 // Solvers
288 // (functions which solve numerically an ODE over an interval)
289
290
307 template <
308 typename Vector, typename OdeFunction = ode_function<Vector>,
309 typename StepFunction = std::function<Vector(OdeFunction, real, const Vector&)>
310 >
313 OdeFunction f, const Vector& x0, real t0, real tf,
314 StepFunction step, real stepsize = 0.001) {
315
316 if (tf < t0) {
317 TH_MATH_ERROR("ode::solve_fixstep", tf, INVALID_ARGUMENT);
318 ode_solution_t<Vector> err; err.t = Vector(nan());
319 return err;
320 }
321
322 const unsigned int steps = floor((tf - t0) / stepsize);
323 unsigned int total_steps = steps;
324 unsigned int i;
325
326 if (abs(t0 + steps * stepsize - tf) > MACH_EPSILON)
327 total_steps++;
328
329 // Initialize solution structure
331
332
333 // Iterate over each step of the numerical method
334 for (i = 1; i <= steps; ++i) {
335 solution.x[i] = step(f, solution.x[i - 1], solution.t[i - 1], stepsize);
336 solution.t[i] = solution.t[i - 1] + stepsize;
337 }
338
339
340 // Additional shorter step if the stepsize
341 // does not cover exactly the time interval
342 if (steps != total_steps) {
343 solution.x[i] = step(
344 f, solution.x[i - 1], solution.t[i - 1], tf - solution.t[i - 1]);
345 solution.t[i] = tf;
346 }
347
348 return solution;
349 }
350
351
363 template <
364 typename Vector, typename OdeFunction = ode_function<Vector>
365 >
372
373
385 template <
386 typename Vector, typename OdeFunction = ode_function<Vector>
387 >
394
395
407 template <
408 typename Vector, typename OdeFunction = ode_function<Vector>
409 >
416
417
430 template <
431 typename Vector, typename OdeFunction = ode_function<Vector>
432 >
439
440
453 template <
454 typename Vector, typename OdeFunction = ode_function<Vector>
455 >
462
463
475 template <
476 typename Vector, typename OdeFunction = ode_function<Vector>
477 >
484 }
485}
486
487#endif
void resize(size_t n) const
Compatibility function to allow for allocation or resizing of dynamic vectors.
Definition vec.h:416
TH_CONSTEXPR unsigned int size() const
Returns the size of the vector (N)
Definition vec.h:406
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
Vector step_euler(OdeFunction f, const Vector &x, real t, real h=0.0001)
Compute one step of Euler's method for ordinary differential equations.
Definition ode.h:140
Vector step_k38(OdeFunction f, const Vector &x, real t, real h=0.001)
Compute one step of Kutta's 3/8 rule method for ordinary differential equations.
Definition ode.h:236
ode_solution_t< vec4 > ode_solution4d
The solution of an ODE in 4 variables.
Definition ode.h:114
Vector step_rk2(OdeFunction f, const Vector &x, real t, real h=0.001)
Compute one step of the Runge-Kutta method of 2nd order for ordinary differential equations.
Definition ode.h:192
std::function< Vector(real, const Vector &)> ode_function
A function representing a system of differential equations, taking as input the time (independent var...
Definition ode.h:123
Vector step_adams3(OdeFunction f, const Vector &x0, real t0, const Vector &x1, real t1, const Vector &x2, real t2, real h=0.001)
Compute one step of the Adams-Bashforth linear multistep method of 3rd order for ordinary differentia...
Definition ode.h:277
ode_solution_t< Vector > solve_midpoint(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.0001)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:388
Vector step_adams2(OdeFunction f, const Vector &x0, real t0, const Vector &x1, real t1, real h=0.001)
Compute one step of the Adams-Bashforth linear multistep method of 2nd order for ordinary differentia...
Definition ode.h:258
ode_solution_t< vec3 > ode_solution3d
The solution of an ODE in 3 variables.
Definition ode.h:110
ode_solution_t< vec2 > ode_solution2d
The solution of an ODE in 2 variables.
Definition ode.h:106
ode_solution_t< Vector > solve_k38(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.0001)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:478
ode_solution_t< Vector > solve_rk4(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.01)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:456
ode_solution_t< vec< real > > ode_solution
The solution of an ODE with any number of variables.
Definition ode.h:98
Vector step_midpoint(OdeFunction f, const Vector &x, real t, real h=0.0001)
Compute one step of the midpoint method for ordinary differential equations.
Definition ode.h:156
ode_solution_t< Vector > solve_fixstep(OdeFunction f, const Vector &x0, real t0, real tf, StepFunction step, real stepsize=0.001)
Integrate an ordinary differential equation using any numerical algorithm with a constant step size,...
Definition ode.h:312
Vector step_rk4(OdeFunction f, const Vector &x, real t, real h=0.01)
Compute one step of the Runge-Kutta method of 4th order for ordinary differential equations.
Definition ode.h:212
ode_solution_t< Vector > solve_heun(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.0001)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:410
ode_solution_t< real > ode_solution1d
The solution of an ODE in 1 variable.
Definition ode.h:102
Vector step_heun(OdeFunction f, const Vector &x, real t, real h=0.001)
Compute one step of Heun's method for ordinary differential equations.
Definition ode.h:172
ode_solution_t< Vector > solve_rk2(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.0001)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:433
ode_solution_t< Vector > solve_euler(OdeFunction f, const Vector &x0, real t0, real tf, real stepsize=0.0001)
Integrate an ordinary differential equation over a certain domain with the given initial conditions u...
Definition ode.h:366
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
double real
A real number, defined as a floating point type.
Definition constants.h:198
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
TH_CONSTEXPR int floor(real x)
Compute the floor of x Computes the maximum integer number that is smaller than x.
Definition real_analysis.h:271
The base type for the solution of an ODE, holding a vector of the values of the time (independent va...
Definition ode.h:23
friend std::ostream & operator<<(std::ostream &out, const ode_solution_t< Vector > &obj)
Stream the ODE solution in string representation to an output stream (std::ostream)
Definition ode.h:88
vec< real > t
A vector of the time values (independent variable).
Definition ode.h:26
ode_solution_t(size_t steps, const Vector &x0, real t0)
Prepare the structure for integration by specifying the number of total steps and the initial conditi...
Definition ode.h:39
ode_solution_t()
Default constructor.
Definition ode.h:33
std::string to_string(const std::string &separator=" ") const
Convert the ODE solution to string representation.
Definition ode.h:52
vec< Vector > x
A vector of the phase space values (dependent variable).
Definition ode.h:29