Differential Equation

Mathematical and statistical software often relies on sequential computations. Examples are likelihood evaluations where it typically is necessary to loop over the rows of the data, or solving ordinary differential equations where numerical approximations are based on looping over the evolving time. When using high-level languages such as R or python such calculations can be very slow unless the algorithms can be vectorized. Fortunately, it is straight-forward to make the implementations in C/C++ and subsequently make an interface to R and python (Rcpp and pybind11). ...

A small illustration on using the armadillo C++ linear algebra library for solving an ordinary differential equation of the form \[ X’(t) = F(t,X(t),U(t)).\] The abstract super class Solver defines the methods solve (for approximating the solution in user-defined time-points) and solveint (for interpolating user-defined input functions on a finer grid). As an illustration a simple Runge-Kutta solver is derived in the class RK4. The first step is to define the ODE, here a simple one-dimensional ODE \(X’(t) = \theta\cdot\{U(t)-X(t)\}\) with a single input \(U(t)\): rowvec dX(const rowvec &input, // time (first element) and additional input variables const rowvec &x, // state variables const rowvec &theta) { // parameters rowvec res = { theta(0)*theta(1)*(input(1)-x(0)) }; return( res ); } The ODE may then be solved using the following syntax odesolver::RK4 MyODE(dX); arma::mat res = MyODE.solve(input, init, theta); with the step size defined implicitly by input (first column is the time variable and the following columns the optional different input variables) and boundary conditions defined by init. ...