File integrators.hpp

Header includes various numerical methods for estimating solutions to stochastic and deterministic ODEs. Methods included in this are: -ODEs: RK4 -SDEs: Euler, Heun

namespace driver

High level interface to differential equation integrators.

Drivers wrap the single step integrators and provide an interface for simulating differential equations for multiple steps from the initial condition. Drivers also handle memory management of the work arrays needed by the integrators.

Author
Oliver Laslett
Date
2017

Functions

void rk4(double *states, const double *initial_state, const std::function<void(double *, const double *, const double)> derivs, const size_t n_steps, const size_t n_dims, const double step_size, )
void eulerm(double *states, const double *initial_state, const double *wiener_process, const std::function<void(double *, const double *, const double)> drift, const std::function<void(double *, const double *, const double)> diffusion, const size_t n_steps, const size_t n_dims, const size_t n_wiener, const double step_size, )
void heun(double *states, const double *initial_state, const double *wiener_process, const std::function<void(double *, double *, const double *, const double)> sde, const size_t n_steps, const size_t n_dims, const size_t n_wiener, const double step_size, )
void implicit_midpoint(double *x, const double *x0, const double *dw, const std::function<void(double *, double *, double *, double *, const double *, const double, const double)> sde, const size_t n_dim, const size_t w_dim, const size_t n_steps, const double t0, const double dt, const double eps, const size_t max_iter, )
namespace integrator

Numerical methods for differential equations.

Numerical methods for simulating the time evolution of deterministic and stochastic ordinary nonlinear differential equations. All integrators compute a single step of the solution.

Author
Oliver Laslett
Date
2017

Functions

void rk4(double *next_state, double *k1, double *k2, double *k3, double *k4, const double *current_state, const std::function<void(double *, const double *, const double)> derivs, const size_t n_dims, const double t, const double h, )
void heun(double *next_state, double *drift_arr, double *trial_drift_arr, double *diffusion_matrix, double *trial_diffusion_matrix, const double *current_state, const double *wiener_steps, const std::function<void(double *, double *, const double *, const double)> sde, const size_t n_dims, const size_t wiener_dims, const double t, const double step_size, )
void eulerm(double *states, double *diffusion_matrix, const double *initial_state, const double *wiener_process, const std::function<void(double *, const double *, const double)> drift, const std::function<void(double *, const double *, const double)> diffusion, const size_t n_dims, const size_t n_wiener, const double t, const double step_size, )
template <class CSTATES, class CDIFF>
void milstein(CSTATES &next_state, const CSTATES &current_state, const CSTATES &drift, const CDIFF &diffusion, const CSTATES &wiener_increments, const double step_size)
int implicit_midpoint(double *x, double *dwm, double *a_work, double *b_work, double *adash_work, double *bdash_work, double *x_guess, double *x_opt_tmp, double *x_opt_jac, lapack_int *x_opt_ipiv, const double *x0, const double *dw, const std::function<void(double *, double *, double *, double *, const double *, const double, const double)> sde, const size_t n_dim, const size_t w_dim, const double t, const double dt, const double eps, const size_t max_iter, )
void rk45(double *next_state, double *temp_state, double *k1, double *k2, double *k3, double *k4, double *k5, double *k6, double *h_ptr, double *t_ptr, const double *current_state, const std::function<void(double *, const double *, const double)> ode, const size_t n_dims, const double eps, )

RK45 Cash-Karp adaptive step deterministic ODE solver.

Solves ODEs :) returns the next state and the time step to be used in the next instance what about the actual time step that it used? That’s also important!

namespace ck_butcher_table

Cash-Karp parameter table for RK45.

Variables

constexpr double c11 = 0.2
constexpr double c21 = 3.0/40.0
constexpr double c22 = 9.0/40.0
constexpr double c31 = 3.0/10.0
constexpr double c32 = -9.0/10.0
constexpr double c33 = 6.0/5.0
constexpr double c41 = -11.0/54.0
constexpr double c42 = 2.5
constexpr double c43 = -70.0/27.0
constexpr double c44 = 35.0/27.0
constexpr double c51 = 1631.0/55296.0
constexpr double c52 = 175.0/512.0
constexpr double c53 = 575.0/13824.0
constexpr double c54 = 44275.0/110592.0
constexpr double c55 = 253.0/4096.0
constexpr double hc1 = 0.2
constexpr double hc2 = 0.3
constexpr double hc3 = 0.6
constexpr double hc4 = 1.0
constexpr double hc5 = 7.0/8.0
constexpr double x11 = 37.0/378.0
constexpr double x13 = 250.0/621.0
constexpr double x14 = 125.0/594.0
constexpr double x16 = 512.0/1771.0
constexpr double x21 = 2825.0/27648.0
constexpr double x23 = 18575.0/48384.0
constexpr double x24 = 13525.0/55296.0
constexpr double x25 = 277.0/14336.0
constexpr double x26 = 0.25