pyRKIntegrator 2.2.0

Description:

pyRKIntegrator 2.2.0

Runge-Kutta and Runge-Kutta-Nystrom integrators
This project contains a series of adaptive generic Runge-Kutta and Runge-Kutta-Nystrom integrators developed in C++, with a Python interface.
The code has been optimized under AVX/AVX2 vectorization directives using the Intel(r) C++ compiler (icpc). Although this is the default compilation mode, it is not strictly required as the makefile will automatically default to GNU if the Intel compilers are not found.
The default optimization mode is fast although it can be changed using the variable OPTL, e.g.,
make OPTL=2

Deployment
A Makefile is provided within the tool to automate the installation for easiness of use for the user. To install the tool simply create a virtual environment as stated below or use the system Python. Once this is done simply type:
make

This will install all the requirements and install the package to your active python. To uninstall simply use
make uninstall

The previous operations can be done one step at a time using
make requirements

to install all the requirements;
make python

to compile and;
make install

to install the tool.
Virtual environment
The package can be installed in a Python virtual environement to avoid messing with the system Python installation.
Next, we will use Conda for this purpose.
Assuming that Conda is already installed, we can create a virtual environment with a specific python version and name (my_env) using
conda create -n my_env python=3.8

The environment is placed in ~/.conda/envs/my_env.
Next we activate it be able to install packages using conda itself or another Python package manager in the environment directory:
conda activate my_env

Then just follow the instructions as stated above.
Runge-Kutta methods
The Runge-Kutta methods are a family of numerical integrators for ODEs. They require a function so that

dydx = f(x,y)

The C++ API requires defining a the function to integrate as
void odefun(double x, double y[], int n, double dydx[])

where y and dydx have n components already allocated. The Runge-Kutta integrator is called as follows
RK_PARAM rkp = rkdefaults(xspan); // Standard Runge-Kutta parameters
RK_OUT rko = odeRK(testfun,xspan,y0,n,rkp); // Runge-Kutta integrator

where rkp contains the ODE integration parameters (see RK.h) and rko is the output structure containing:

A return value, rko.retval, that is > 0 in case of success or < 0 in case of failue. A value of 0 indicates that the odeRK routine has not been run.
The number of integration steps rko.n.
An error value rko.err.
The solution, rko.x and rko.y.
The available Runge-Kutta schemes are:
Euler-Heun 1(2) (eulerheun12)
Bogacki-Shampine 2(3) (bogackishampine23)
Fehlberg 4(5) (fehlberg45)
Cash-Karp 4(5) (cashkarp45)
Dormand-Prince 4-5 (dormandprince45)
Calvo 5(6) (calvo56)
Dormand-Prince 7(8) (dormandprince78)
Curtis 8(10) (curtis810)
Hiroshi 9(12) (hiroshi912)

Runge-Kutta-Nystrom methods
The Runge-Kutta-Nystrom methods are a family of numerical integrators for smooth second order ODEs. They require a function so that

dy2dx2 = f(x,y)

The C++ API requires defining a the function to integrate as
void odefun(double x, double y[], int n, double dy2dx2[])

where rkp contains the ODE integration parameters (see RK.h) and rko is the output structure containing:

A return value, rko.retval, that is > 0 in case of success or < 0 in case of failue. A value of 0 indicates that the odeRK routine has not been run.
The number of integration steps rko.n.
An error value rko.err.
The solution, rko.x, rko.y and rko.dy.
The available Runge-Kutta-Nystrom schemes are:
Runge-Kutta-Nystrom 3(4) (rkn34)
Runge-Kutta-Nystrom 4(6) (rkn46)
Runge-Kutta-Nystrom 6(8) (rkn68)
Runge-Kutta-Nystrom 10(12) (rkn1012)

C and C++ implementations
There is a C and a C++ implementation of the algorithms. The language can be chosen in the Makefile by setting the USE_CPP variable. The compilation of the examples and the python tools will follow.
make USE_CPP=ON/OFF

The Python interface
The Python interface allows using the Runge-Kutta and Runge-Kutta-Nystrom integrators using Python's cython. The wrapper can be compiled using
make python

Note that this will create a .so that must be in the same directory of RungeKutta.py. Otherwise, the Python interface is accessed without the speedup of compiled code.
In Python, the Runge-Kutta module should be first using
import pyRKIntegrator as rk

The function to integrate must be defined as
def odefun(x,y,n,dydx): # For Runge-Kutta
def odefun(x,y,n,dy2dx): # For Runge-Kutta-Nystrom

where y and dydx or dy2dx are vectors that have already been allocated (of size n). Access to the parameters structure is provided by the odeset class:
params = rk.odeset() # For Runge-Kutta and Runge-Kutta-Nystrom

This will initiate using default parameters, which can be changed by using key arguments or by modifying the class fields. The access to the integratos is:
x,y,err = rk.odeRK(odefun,xspan,y0,odeset) # For Runge-Kutta
x,y,dy,err = rk.odeRKN(odefun,xspan,y0,dy0,odeset) # For Runge-Kutta-Nystrom

Examples 1 to 4 include examples of advanced features to use with the integrator.
Event and output functions
Event detection and output functions are a characteristic of this set of integrators. They can be set using the RK_PARAM structure in C++
RK_PARAM rkp;
rkp.eventfcn = eventfun;
rkp.outputfcn = outputfun;

or the odeset class in Python
odeset = rk.odeset(eventfun=eventfun,outputfun=outputfun)

Output functions
Output functions allow retrieving the integrated values right after a successful integration step. A value of either 1 or 0 must be returned to continue or stop the integration. The implementation in C++ is
int outputfun(double x, double y[], int n) {
/* Some code */
return 1; // or 0 to stop integration
}

and in Python
def outputfun(x,y,n):
# Some code
return 1 # or 0 to stop integration

Event functions
Event functions allow for solving a problem so that

g(x) - val = 0

using a root solving algorithm. This is useful, for example, to detect when the integrator reaches a certain value (see example2.py for a more comprehensive usage). The implementation in C++ is
int eventfun(double x, double y[], int n, double value[1],int direction[1]) {
/* Some code */
return 1; // or 0 to stop integration
}

and in Python
def eventfun(x,y,n,value,direction):
# Some code
return 1 # or 0 to stop integration

Where

value is a mathematical expression describing the event. An event occurs when value(i) is equal to zero.
direction

0 if all zeros are to be located.
+1 locates only zeros where the event function is increasing.
-1 locates only zeros where the event function is decreasing.
Regarding direction, only 0 is implemented.