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

[luguecos]

Hello, 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 NULL; }
}

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

    for (int j = 0; j < sistvar; j++) 
    {
        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
        x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
        k4[j] = step * F(t + step, x_temp3, 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.

This is a way I tried but it didn't work:

double system_eq(double t, double x[], int eq)
{
    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 NULL; }
}

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

    for (int j = 0; j < sistvar; j++) {

        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(time, x, j));
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(time + 0.5 * step, x_temp1, j));
        x_temp3[j] = x[j] + (k3[j] = step * F(time + 0.5 * step, x_temp2, j));
        k4[j] = step * F(time + step, x_temp3, j);


        x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6;

    }

    t += step;
}

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!

[salem c]

Wrap things up in a class.

#include <iostream>
#include <cmath>
using namespace std;

class RungeKutta {
  // All your globals are now class members
  double gama, OMEGA, zeta, alpha, beta, chi, kappa, phi;
  int sistvar;

  // These could be here, or back in your rk4 function
  double *x_temp1;
  double *x_temp2;
  double *x_temp3;
  double *k1;
  double *k2;
  double *k3;
  double *k4;
public:
  RungeKutta(double g, double O, double z,
             double a, double b, double c,
             double k, double p, int s)
    : gama(g),
      OMEGA(O),
      zeta(z),
      alpha(a),
      beta(b),
      chi(c),
      kappa(k),
      phi(p),
      sistvar(s) {
      // create all those arrays based on sistvar size
      x_temp1 = new double[sistvar];
      x_temp2 = new double[sistvar];
      x_temp3 = new double[sistvar];
      k1 = new double[sistvar];
      k2 = new double[sistvar];
      k3 = new double[sistvar];
      k4 = new double[sistvar];
    }
  ~RungeKutta() {
    // and clean them up
    delete [] x_temp1;
    delete [] x_temp2;
    delete [] x_temp3;
    delete [] k1;
    delete [] k2;
    delete [] k3;
    delete [] k4;
  }
  double F(double t, double x[], int eq);
  void rk4(double &t, double x[], double step);
};

double RungeKutta::F(double t, double x[], int eq) {
  // gama etc are now your class member variables
  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; }  //!! NULL is a pointer, not a double
}

void RungeKutta::rk4(double &t, double x[], double step) {
  for (int j = 0; j < sistvar; j++)
  {
    x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
    x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
    x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
    k4[j] = step * F(t + step, x_temp3, j);

    x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
  }

  t += step;
}

int main()
{
  // create an instance, passing in your 9 globals.
  RungeKutta myrk(1,2,3,4,5,6,7,8,9);
  double t = 0;
  double x[10];
  myrk.rk4(t,x,1.0);  // call the rk4 method
}

The next move might be to replace all that manual memory management by using say
http://www.cplusplus.com/reference/vector/vector/vector/
std::vector<double> x_temp1(sistvar);