2

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
Luguecos
  • 131
  • 1
  • 1
  • 9
  • 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 Dec 08 '19 at 11:07
  • 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 Dec 08 '19 at 11:10
  • You are using C++. Have you contemplated using the std::vector class to implement the state vectors? – Lutz Lehmann Dec 08 '19 at 11:11
  • Hello Dr. Lehmann, thank you for your response. As I mentioned, the first code works and it is already validated. What I coudn't manage to do is to make it modular (as function in a separate .cpp file or a class) and split it from the main code. I thought about implementing with std::vector but I would need to change too much of the code to do it. I'm looking to do minor changes in the code to accomplish this goal as the code already works and it is validated. – Luguecos Dec 08 '19 at 14:41
  • No, the code might compile and run, but it does not implement RK4 for systems. If you validate with scalar equations or uncoupled systems you might miss that. Just think about what the value of `x_temp1[2]` is at the time it is used in the round `j=1` of the loop. – Lutz Lehmann Dec 08 '19 at 15:37
  • Dr. Lehmann, **if** I **understand** you **correctly**, what you are asking me to think about is a impossible situation as `x_temp1[2]` is only possible when `j = 2`. `j` defines which equation of the system the `F` function returns. – Luguecos Dec 08 '19 at 15:56
  • I validated this code with linear oscillator system, duffing oscillator system and some other specific coupled systems. – Luguecos Dec 08 '19 at 16:12
  • 1
    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 Dec 08 '19 at 21:04
  • 1
    For some C++ code with lots of boilerplate look at the answer to https://stackoverflow.com/q/49738918/3088138, or a more bare-bones implementation at https://people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html – Lutz Lehmann Dec 08 '19 at 21:40
  • Hello Dr. Lehmann. I did take a look in the material you presented and I really can't see the difference between my code and the other 2 listed. Mine seems to be only more compact, but does the same thing. It calculates the first inclination, with the value of the first it calculates the second and so on until k4. For all equations (j= 0, 1 and 2 - in this case). The only difference is that I calculate in a single `for` loop while theirs use one loop for each inclination. – Luguecos Dec 09 '19 at 19:46
  • No, think about it again. In the first traversal of the loop you compute everything for `j=0`, that is, evaluating the first equation and updating the first components of the vectors. All the other vector components have "stale" values from the previous step (or random values in the first step). This is really not what you want for a vector version. You can not avoid the separate loops, you can only separate them out as `axpy` (the vector equivalent of FMA) operations. – Lutz Lehmann Dec 09 '19 at 19:55
  • I used to write this code with 4 `for` loops, but recently I changed to a more compact version with 1 loop. Is that the problem? (I edited the post and put also the older version in the end of the question). – Luguecos Dec 09 '19 at 19:56
  • 1
    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 Dec 09 '19 at 19:56
  • Yes, that is the correct way that gives you order 4. – Lutz Lehmann Dec 09 '19 at 19:57
  • Now I understand. I did some calculations and it is clarified. Thank you for your help! I'm going to write the corrected version in the thread. – Luguecos Dec 09 '19 at 20:07
  • I found a topic that implements the C++ hints I gave above, https://stackoverflow.com/q/51840543/3088138. The only issue here is that the k vectors are constantly allocated and deleted, it would be better to re-use the once-allocated vectors by passing them as argument to `Func`. (Or not, if you are building a "dense output" solution where the k vectors constitute the interpolation coefficients.) – Lutz Lehmann Dec 09 '19 at 23:38
  • I took a look in the link you provided and after reading it I tried to implement my code using std::vector as it can be sized in running time. I almost tried to implement a class with it, however I decided to do a performence test before and turns out that the vector version takes 2.64x to run in comparison to the array version. I'm doing something wrong? I did run the code using Intel c++ compiler with -O2 optimization in both cases. The vector version is in the EDITED part of the post as the data of the performance test. – Luguecos Dec 11 '19 at 06:01
  • Yes, that is normal. std::vector is more an array list than a mathematical vector. Array access there is a method of the class, with bounds check and automatic extension of the array. That is why I looked for the last link for a code that uses a vector in the linear algebra sense. The boost::eigen library does that. – Lutz Lehmann Dec 11 '19 at 07:51

1 Answers1

0

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51