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; }
Improve this question
$\endgroup$
2
  • $\begingroup$ Can you just change the Function() definition to have the times as inputs? Something like: Function( double initialTime, double endTime ) $\endgroup$ Commented May 7, 2016 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$ Commented May 7, 2016 at 18:13

1 Answer 1

Reset to default
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.

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.