使用C与GSL库求解阻尼摆ODE系统初期结果异常求助

我用C语言结合GSL库实现了简单阻尼摆的常微分方程（ODE）求解代码，但运行时前十几轮计算结果异常，出现73.609这类不符合初始设置的位置值；当位置数值达到100后，计算结果恢复正常且可正常绘图。相关代码及输出结果如下：

struct params{
    double g;
    double length;
    double damping;
    double startPosition;
    int startVelocity;
};

int func(double t, const double y[], double f[], void *params){
    (void) (t);
    struct params *init  = (struct params*)params;
    // 速度导数（位置的一阶导）
    f[0] = y[1];
    // 加速度（位置的二阶导）
    f[1] = -(init->damping)*y[1] - (init->g / init->length)*sin(y[0]);
    return GSL_SUCCESS;
}

int main(){
    // 参数初始化
    struct params *inputs = malloc(sizeof(struct params));
    
    if (inputs != NULL){
        inputs->g=9.8;
        inputs->length=2;
        inputs->damping=0.1;
        inputs->startPosition = 0.2;
        inputs->startVelocity = 0;
    }

    int i;
    double t = 0.0, t1 = 100;
    double dt = 1;
    const int step_total = 101;

    double* timeX = malloc(step_total*sizeof(double));
    double* yOne = malloc(step_total*sizeof(double));
    double* yTwo = malloc(step_total*sizeof(double));

    clock_t start_time = clock(); // 记录开始时间
    double timeout = 5.0; // 超时时间（秒）

    int SCREEN_WIDTH = 1280;
    int SCREEN_HEIGHT = 720;


    // 系统初始化
    gsl_odeiv2_system system = {func, NULL, 2, inputs};

    // ODE求解器初始化
    const gsl_odeiv2_step_type *type = gsl_odeiv2_step_rk8pd;
    gsl_odeiv2_step *step = gsl_odeiv2_step_alloc(type, 2);
    gsl_odeiv2_control *control = gsl_odeiv2_control_y_new(1e-6,0.0);
    gsl_odeiv2_evolve *evolve = gsl_odeiv2_evolve_alloc(2);

    double y[2] = {inputs->startPosition, inputs->startVelocity};

    //printf ( " Step method is '%s'\n", gsl_odeiv2_step_name(step));
    i = 0; 
    while(t < t1){

        int status = gsl_odeiv2_evolve_apply(evolve, control, step, &system, &t, t1, &dt, y);
        
        if (status != GSL_SUCCESS){
            printf("error, return value=%d\n", status);
            break;
        }
        //if (y[0] <= 0) break;
        //printf("At time t=%.3f s, position=%.3f m, velocity=%.3f m/s\n", t, y[0], y[1]);
        timeX[i] = t;
        yOne[i] = y[0];
        yTwo[i] = y[1];
        i++;


        clock_t current_time = clock();
        double elapsed_time = (double)(current_time - start_time) / CLOCKS_PER_SEC;
        if (elapsed_time > timeout){
            printf(" Infinite loop timeout reached. \n");
            break;
        } 

    }
    int total = i-1;
    i = 0;
    gsl_odeiv2_evolve_reset(evolve);
    gsl_odeiv2_control_free(control);
    gsl_odeiv2_step_free(step);
    gsl_odeiv2_evolve_free(evolve);
}

输出结果示例：

At time t=0.575 s, position=73.609 m, velocity=0.004 m/s
At time t=1.185 s, position=74.530 m, velocity=0.000 m/s
……
At time t=16.675 s, position=100.000 m, velocity=0.000 m/s
At time t=17.320 s, position=0.071 m, velocity=-0.103 m/s
……

请问为何系统会出现这类异常初始值？

原因分析与解决办法 #

核心原因：数组越界导致内存损坏 #

代码中step_total被设置为101，意味着最多只能存储101个数据点，但你使用的rk8pd是自适应步长求解器，它会根据计算精度自动缩小步长（比如初始步长设为1，但实际可能调整到0.5甚至更小），导致实际计算步数远超过101。此时timeX、yOne、yTwo数组会发生越界写入，破坏内存中其他区域的数据——比如覆盖了inputs结构体里的阻尼系数、重力加速度等参数。

当阻尼系数被意外修改为负数时，原阻尼振荡系统会变成增幅振荡系统，位置数值会迅速发散增大；后续内存写入操作可能又覆盖回正确的参数值，计算就恢复了正常的阻尼振荡行为，这就是你看到数值先异常增大后恢复的原因。

解决措施 #

1. 扩大数组分配大小：根据预估的最大步数调整step_total，比如改为1000或更大，确保能容纳所有求解产生的数据点；或者改用动态数组（如每次达到容量时重新扩容）。
2. 添加越界检查：在循环中加入if (i >= step_total) break;，避免数组越界写入。
3. 统一参数类型：将struct params中的startVelocity从int改为double，保持和y[1]的类型一致，避免潜在的类型转换问题。

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