0
$\begingroup$

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;
    }
    
$\endgroup$
  • $\begingroup$ Can you just change the Function() definition to have the times as inputs? Something like: Function( double initialTime, double endTime ) $\endgroup$ – spektr May 7 '16 at 16:11
  • $\begingroup$ Yes, I can. But what will be the expression/value for the initial and final time? Is there anyway to use the function as defined? $\endgroup$ – TheCoder May 7 '16 at 18:13
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ @Chooward : Kindly see my updated question, where I define the variables under questions and how they are called/used inside the Function() $\endgroup$ – TheCoder May 9 '16 at 8:26
  • $\begingroup$ @TheCoder Made some changes based on your added info. $\endgroup$ – spektr May 9 '16 at 16:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.