Implementing modular Runge-kutta 4th order method for a n-dimension system

Question

i'm trying to make my runge-kutta 4th order code modular. I don't want to have to write and declare the code everytime I use it, but declare it in a .hpp and a .cpp file to use it separetely. But i'm having some problems. Generally I want to solve a n-dimension system of equations. For that I use two functions: one for the system of equations and another for the runge-kutta method as follows:

double F(double t, double x[], int eq)
{
    // System equations
    if      (eq == 0) { return (x[1]); }
    else if (eq == 1) { return (gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]); }
    else if (eq == 2) { return (-kappa * x[1] - phi * x[2]); }
    else { return 0; }
}

void rk4(double &t, double x[], double step)
{
    double x_temp1[sistvar], x_temp2[sistvar], x_temp3[sistvar];        
    double k1[sistvar], k2[sistvar], k3[sistvar], k4[sistvar];      

    int j;                                                      

    for (j = 0; j < sistvar; j++) 
    {
        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        k4[j] = step * F(t + step, x_temp3, j);
    }
    for (j = 0; j < sistvar; j++) 
    {
        x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
    }

    t += step;
}

The above code works and it is validated. However it has some dependencies as it uses some global variables to work:

gama, OMEGA, zeta, alpha, beta, chi, kappa and phi are global variables that I want to read from a .txt file. I already manage to do that, however only in a single .cpp file with all code included.

Also, sistvar is the system dimension and also a global variable. I'm trying to enter it as an argument in F. But the way it is written seems to give errors as sistvar is a const and can't be changed as a variable and I can't put variables inside an array's size.

In addition, the two functions has an interdependency as when a call F inside rk4, eq number is needeed.

Could you give me tips in how to do that? I already searched and read books about this and could not find an answer for it. It is probably an easy task but i'm relatively new in c/c++ programming languages.

Thanks in advance!

* EDITED (Tried to implement using std::vector)*

double F(double t, std::vector<double> x, int eq)
{
    // System Equations
    if (eq == 0) { return (x[1]); }
    else if (eq == 1) { return (gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]); }
    else if (eq == 2) { return (-kappa * x[1] - phi * x[2]); }
    else { return 0; }
}

double rk4(double &t, std::vector<double> &x, double step, const int dim)
{
    std::vector<double> x_temp1(dim), x_temp2(dim), x_temp3(dim);
    std::vector<double> k1(dim), k2(dim), k3(dim), k4(dim);

    int j;

    for (j = 0; j < dim; j++) {
        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
    }
    for (j = 0; j < dim; j++) {
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
    }
    for (j = 0; j < dim; j++) {
        x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
    }
    for (j = 0; j < dim; j++) {
        k4[j] = step * F(t + step, x_temp3, j);
    }

    for (j = 0; j < dim; j++) {
        x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
    }

    t += step;

    for (j = 0; j < dim; j++) {
        return x[j];
    }
}

vector                 array

2.434 s   |        |   0.859 s
2.443 s   |        |   0.845 s
2.314 s   |        |   0.883 s
2.418 s   |        |   0.884 s
2.505 s   |        |   0.852 s
2.428 s   |        |   0.923 s
2.097 s   |        |   0.814 s
2.266 s   |        |   0.922 s
2.133 s   |        |   0.954 s
2.266 s   |        |   0.868 s
_______                _______
average = 2.330 s      average = 0.880 s

I highly doubt that the code works as a code for the RK4 method. You only get an order 1 method, not the expected order 4 method. Test it, compute the error for a non-scalar test problem with known solution for different step sizes. — Lutz Lehmann
You are using C++. Have you thought about implementing the ODE system as a class that exports the dimension and ODE function (via some interface or base class)? — Lutz Lehmann
Ok, the oscillator is a circularly coupled system. Did you also validate numerically the order of the method? As I said, the method has still order 1, and the arrays will be filled with order 1 correct values after the first step, so you still get visually correct results if the step size is small enough. — Lutz Lehmann
For some C++ code with lots of boilerplate look at the answer to stackoverflow.com/q/49738918/3088138, or a more bare-bones implementation at people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html — Lutz Lehmann
With some clever pointer handling you could compactify the stages of RK4 in a loop, using coefficient arrays c={ 0,0.5, 0.5, 1}; and b={1,2,2,1}; to simultaneously evaluate derivatives and update the next state vector. — Lutz Lehmann

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2019-12-11T12:57:32

Using vector function where the vector arithmetic is taken from Eigen3

#include <eigen3/Eigen/Dense>
using namespace Eigen;

of the same parts as discussed in the question could look like (inspired by function pointer with Eigen)

VectorXd Func(const double t, const VectorXd& x)
{ // equations for solving simple harmonic oscillator
    Vector3d dxdt;
    dxdt[0] = x[1];
    dxdt[1] = gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]; 
    dxdt[2] = -kappa * x[1] - phi * x[2]; 
    return dxdt;
}

MatrixXd RK4(VectorXd Func(double t, const VectorXd& y), const Ref<const VectorXd>& y0, double t, double h, int step_num)
{
    MatrixXd y(y0.rows(), step_num );
    VectorXd k1, k2, k3, k4;
    y.col(0) = y0;

    for (int i=1; i<step_num; i++){
        k1 = Func(t, y.col(i-1));
        k2 = Func(t+0.5*h, y.col(i-1)+0.5*h*k1);
        k3 = Func(t+0.5*h, y.col(i-1)+0.5*h*k2);
        k4 = Func(t+h, y.col(i-1)+h*k3);
        y.col(i) = y.col(i-1) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
        t = t+h;
    }
    return y.transpose();
}

Passing a vector to a function to be filled apparently requires some higher template contemplations in Eigen.

Implementing modular Runge-kutta 4th order method for a n-dimension system

1 Answers