Equation of motion by RK4 method

Question

I would like to sole the spring mass equation of the form: $M\ddot{Y}+B\dot{Y}+KY = F_y(t)$, where $M :$ mass, $B :$ Damping, $K : $ stiffness, $F_y :$ external force and $Y :$ displacement in transverse direction; using Runge-Kutta 4th-order method. I modify RK4 C++ code as:

/*
Definition of the y'(t) = f1(t,y,y') = y' by the definition
*/
double f1(double t, double y, double v)
{
  double d1y;
  d1y = v;
  return d1y;
}

/*
Definition of the y"(t) = f2(t,y,y')
*/
double f2(double t, double y, double v, double Fy)
{
 double d2y;
 double M = 0.03715; // 10.0;
 double B = 0.0;
 double K = 0.149; 
 d2y = (-Fy-B*v-K*x)/M; 
 return d2y;
}

/*-----------------------------------------------------------------------
Program to solve 1st order Ordinary Differential Equations
y'(t) = d1y(t,y)
y"(t) = d2y(t,y)
method: 4th-order Runge-Kutta method */

double rk4_2nd(double(*d1x)(double, double, double),
           double(*d2x)(double, double, double, double),
           double ti, double xi, double vi, double tf, double Fy,
           double& xf, double& vf)
{
  double h,t,k1x,k2x,k3x,k4x,k1v,k2v,k3v,k4v;

  h = tf-ti;
  t = ti;

  k1x = h*d1x(t,xi,vi);
  k1v = h*d2x(t,xi,vi,Fy);

  k2x = h*d1x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0);
  k2v = h*d2x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0,Fy);

  k3x = h*d1x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0);
  k3v = h*d2x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0,Fy);

  k4x = h*d1x(t+h,xi+k3x,vi+k3v);
  k4v = h*d2x(t+h,xi+k3x,vi+k3v,Fy);

  xf = xi + (k1x + 2.0*(k2x+k3x) + k4x)/6.0;
  vf = vi + (k1v + 2.0*(k2v+k3v) + k4v)/6.0;

  return 0.0;
}

The problem I am facing is that:

• RK4 is called inside a function e.g.Function(void) and this Function() is called inside main where a time loop runs as form of iterations, How can I define ti, tf, to RK4 function for my problem.

   void Function() {
  
   Y0[n] = Y[n] - YC; // n =  number of objects i.e. n = 1 at the moment
   // Y : current displacement; YC: Initial/reference position = constant
   V0[n] = V[n]; // Velocity at (t-1) = velocity at (t)
   FY[n] += -k * (Y0[n]); // external force; k : spring constant
  
   rk4_2nd(f1, f2, t0, Y0[n], V0[n], t1, FY[n], Y[n], V[n]); 
  } 
  
  int main() {
    N = 0;
  do {
    N = N + 1;
    Function();
   } while (N < NMAX);
   return 0;
  }
  

Answer 1

This update to your code should be what you want if I understand how you want to implement it.

void Function(double currentTime, double nextTime) {
 for( int i = 0; i < n; i++ ){
    Y0[i] = Y[i] - YC; // n =  number of objects i.e. n = 1 at the moment
    // Y : current displacement; YC: Initial/reference position = constant
    V0[i] = V[i]; // Velocity at (t-1) = velocity at (t)
    FY[i] += -k * (Y0[i]); // external force; k : spring constant

    rk4_2nd_me1(f1, f2, currentTime, Y0[i], V0[i], nextTime, FY[i], Y[i], V[i]); 
 }
} 

int main() {
  double currentTime = 0.0, finalTime = 1.0;
  double dt = (finalTime-currentTime)/static_cast<double>(NMAX);
  int N = 0;
  do {
    currentTime = N*dt;
    N = N + 1;
    Function(currentTime,currentTime+dt);
  } while (N < NMAX);
 return 0;
}

Now based on this, I assume that Y0, V0, Y, V, FY, YC are all global variables.