# Equation of motion by RK4 method

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;
}

• Can you just change the Function() definition to have the times as inputs? Something like: Function( double initialTime, double endTime ) – spektr May 7 '16 at 16:11
• 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? – TheCoder May 7 '16 at 18:13

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.

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