GSL ODE求解器报错：积分限/步长方向不一致问题求助

使用GSL库求解洛伦兹力ODE的C++代码，单个粒子模拟正常，但粒子数大于1时返回错误：ERROR: integration limits and/or step direction not consistent。代码通过整数参数控制粒子数量，每步输出粒子位置(x,y,z)和速度(vx,vy,vz)，初始值由GSL随机函数生成。

核心问题出在积分器驱动的生命周期管理与时间变量重置：

• 第一次粒子循环结束时调用了gsl_odeiv2_driver_free(drv);，直接释放了驱动对象；第二次循环时drv已变为野指针，继续使用会导致时间步状态混乱，触发积分区间错误
• 第一个粒子模拟结束后，时间变量t停在接近tf的位置，第二个粒子模拟时未重置t，导致t_next从t0+dt开始，与t的当前值冲突，加剧了步长方向不一致的问题

1. 调整驱动对象的生命周期：将驱动的创建放在粒子循环外，释放放在循环结束后，避免重复创建/释放导致的野指针问题
2. 重置时间变量：每个粒子初始化完成后，将t重置为t0，确保积分从初始时间开始
3. 状态数组隔离：将system数组的声明移到粒子循环内，避免残留上一个粒子的状态

#include<iostream>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

using namespace std;

int lorentzODE (double t, const double s[], double f[], void *params);

int main(){

   int noPrimaries = 2;
   double q = 1.0;                           // 粒子电荷
   double m = 1.0;                           // 粒子质量
   double B[3] = {0.0, 0.0, 1.0};            // 磁场
   double E[3] = {0.1, 0.0, 0.5};            // 电场

   double x0[3];                 // 初始位置
   double v0[3];                 // 初始速度

   const gsl_rng_type * T = gsl_rng_ranlxs0;
   gsl_rng * r = gsl_rng_alloc(T);
   gsl_rng_set(r, (unsigned long) time(NULL));

   double t0 = 0.0;
   double tf = 100.0;
   double dt = 0.05;

   int status;                       // 驱动函数返回状态
   double paramsB[8];
   paramsB[0] = q;
   paramsB[1] = m;
   paramsB[2] = B[0];
   paramsB[3] = B[1];
   paramsB[4] = B[2];
   paramsB[5] = E[0];
   paramsB[6] = E[1];
   paramsB[7] = E[2];

   const string &solverMethod = "RK8PD";
   const double h = 1.0e-06;
   const double epsAbs = 1.0e-08;
   const double epsRel = 1.0e-10;

   // 积分器配置
   double t, t_next;
   gsl_odeiv2_system odeSystem;
   odeSystem.function = lorentzODE;
   odeSystem.dimension = 6;
   odeSystem.params = paramsB;

   // 创建积分器驱动（放在循环外）
   gsl_odeiv2_driver *drv;
   drv = gsl_odeiv2_driver_alloc_y_new(&odeSystem, gsl_odeiv2_step_rk8pd, h, epsAbs, epsRel);

   for(int i=0;i<noPrimaries;i++){
      // 每个粒子单独声明状态数组，避免状态残留
      double system[6];

      x0[0] = gsl_rng_uniform (r);
      x0[1] = gsl_rng_uniform (r);
      x0[2] = gsl_rng_uniform (r);
      v0[0] = gsl_rng_uniform (r);
      v0[1] = gsl_rng_uniform (r);
      v0[2] = gsl_rng_uniform (r);
      system[0] = x0[0]; system[1] = x0[1]; system[2] = x0[2];
      system[3] = v0[0]; system[4] = v0[1]; system[5] = v0[2];

      // 重置时间变量为初始值
      t = t0;

      for(t_next = t0+dt; t_next<tf; t_next += dt){
         
         status = gsl_odeiv2_driver_apply(drv, &t, t_next, system);
         if(status != GSL_SUCCESS){
            printf("Error: status = %d\n", status);
            break;
         }
         printf("%.5e %.5e %.5e %.5e %.5e %.5e %.5e\n", t, system[0], system[1], system[2], system[3], system[4], system[5]);
      }
   }
   // 释放驱动（放在循环结束后）
   gsl_odeiv2_driver_free(drv);
   gsl_rng_free(r); // 释放随机数生成器
}

int lorentzODE (double t, const double s[], double f[], void *params){
  (void)(t); /* 避免未使用参数警告 */
  double *lparams = (double *)params;

  double q  = lparams[0];
  double m  = lparams[1];
  double mu = q/m;

  double Bx = lparams[2];
  double By = lparams[3];
  double Bz = lparams[4];

  double Ex = lparams[5];
  double Ey = lparams[6];
  double Ez = lparams[7];

  f[0] = s[3];
  f[1] = s[4];
  f[2] = s[5];
  f[3] = mu*(Bz*s[4] - By*s[5] + Ex);
  f[4] = mu*(Bx*s[5] - Bz*s[3] + Ey);
  f[5] = mu*(By*s[3] - Bx*s[4] + Ez);
  return GSL_SUCCESS;
}

内容的提问来源于stack exchange，提问作者lghizoni