我想唯一的形式的弹簧质量方程： ，其中质量，阻尼，刚度，外部力和横向位移；使用龙格-库塔四阶法。我将 RK4 C++ 代码修改为：$M\ddot{Y}+B\dot{Y}+KY = F_y(t)$$M :$$B :$$K :$$F_y :$$Y :$

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

我面临的问题是：

• RK4 在函数 egFunction(void) 内部被调用，而这个 Function() 在 main 内部被调用，其中时间循环以迭代的形式运行，我如何为我的问题定义 ti、tf 到 RK4 函数。

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