C++中数组实现四阶龙格-库塔的变量声明及异常排查

问题概述 #

模拟耦合振子链时，使用数组结合四阶龙格-库塔（RK4）算法处理多振子，当前因变量声明和函数调用逻辑错误导致计算结果发散：

• 当f2函数依赖索引i时结果发散，不依赖i的测试函数（如返回1）运行正常
• 推测问题与指针使用不当有关，无指针时程序无法运行

现有代码 #

头文件 func_solition.h #

#include <iostream>
#include <cmath>
using namespace std;

double f1(double theta1) {
    return theta1;
}

double f2(double t, double theta[], double alpha, double beta, int i) {
    //return sin(theta[i])-beta*(theta[i+1]+theta[i-1]-2*theta[i]);
    return sin(theta[i]);
}

主文件 #

#include "func_solition.h"
#include <iostream>
using namespace std;

int main() {
// Paramètres intégrale 
    double dt = 0.01;
    double t = 0;
    double TF = 20;

    // CI pendules
    int n = 10;                      // ATTENTION ON A DEUX PENDULES POUR LES BC, PREVOIR +2 PENDULES
    double theta[n];
    double theta1[n];
    for (int i = 0; i < n; i++) {
        theta[i] = 0;
    }
    theta[1] = 3.141592/2;
    for (int i = 0; i < n; i++) {
        theta1[i] = 0;
    }
    theta1[1] = 0;

    // Constantes équation du mouvement
    double alpha = 0.1;
    double beta = 0.01;

    // Variables pour rk4
    double tmp1=0.0, tmp2=0.0;
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;

    cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;

    // rk4
    while (t <= TF) {
        for (int i = 1; i < n-2; ++i) {
            k1 = dt * f1(theta1[i]);
            l1 = dt * f2(t, theta, alpha, beta, i);

            k2 = dt * f1(theta1[i] + 0.5 * l1);
            tmp2 = theta[i] + 0.5 * k1;
            tmp1 = t + 0.5 * dt;
            l2 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k3 = dt * f1(theta1[i] + 0.5 * l2);
            tmp2 = theta[i] + 0.5 * k2;
            tmp1 = t + 0.5 * dt;
            l3 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k4 = dt * f1(theta1[i] + l3);
            tmp2 = theta[i] + k3;
            tmp1 = t + dt;
            l4 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            theta[i] = theta[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            theta1[i] = theta1[i] + (l1 + 2 * l2 + 2 * l3 + l4) / 6;

        }

// Boundary conditions
        theta[n-1] = 0;
        theta1[n-1] = 0;
        theta[0] = 0;
        theta1[0] = 0;

        t = t + dt;
        cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
        
    }
}

核心问题分析 #

1. RK4中间步骤的数组传递错误：
   计算l2/l3/l4时，传递&tmp2（单个double变量的地址）给f2的theta[]参数，但f2中尝试访问theta[i]——此时指针指向单个值而非完整数组，访问theta[i]会导致越界访问未定义内存，这是结果发散的根本原因。
2. 循环范围错误：
   主循环for (int i = 1; i < n-2; ++i)，当n=10时仅循环到i=7，无法处理倒数第二个振子（最后一个为边界条件），正确范围应为i < n-1。
3. 变长数组的非标准性：
   double theta[n];是C变长数组（VLA），属于非标准扩展，建议改用std::vector<double>保证可移植性。

修复方案 #

1. 修正RK4中间步骤的状态传递逻辑 #

RK4每个阶段需要构造完整的中间状态数组，复制当前theta数组并修改对应索引值后，再传递给f2。

2. 调整循环范围 #

将循环条件改为for (int i = 1; i < n-1; ++i)，确保所有非边界振子都被处理。

3. 替换变长数组为vector（可选但推荐） #

修复后的主文件代码 #

#include "func_solition.h"
#include <iostream>
#include <vector>
using namespace std;

int main() {
    // Paramètres intégrale 
    double dt = 0.01;
    double t = 0;
    double TF = 20;

    // CI pendules
    int n = 10;
    vector<double> theta(n, 0.0);
    vector<double> theta1(n, 0.0);
    theta[1] = 3.141592/2;
    theta1[1] = 0;

    // Constantes équation du mouvement
    double alpha = 0.1;
    double beta = 0.01;

    // Variables pour rk4
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;

    cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;

    // rk4
    while (t <= TF) {
        // 保存当前状态，避免RK4中间步骤修改原数组影响其他振子计算
        vector<double> theta_current = theta;
        vector<double> theta1_current = theta1;

        for (int i = 1; i < n-1; ++i) {
            k1 = dt * f1(theta1_current[i]);
            l1 = dt * f2(t, theta_current.data(), alpha, beta, i);

            // 构造k2对应的中间theta状态
            vector<double> theta_k2 = theta_current;
            theta_k2[i] += 0.5 * k1;
            k2 = dt * f1(theta1_current[i] + 0.5 * l1);
            l2 = dt * f2(t + 0.5*dt, theta_k2.data(), alpha, beta, i);

            // 构造k3对应的中间theta状态
            vector<double> theta_k3 = theta_current;
            theta_k3[i] += 0.5 * k2;
            k3 = dt * f1(theta1_current[i] + 0.5 * l2);
            l3 = dt * f2(t + 0.5*dt, theta_k3.data(), alpha, beta, i);

            // 构造k4对应的中间theta状态
            vector<double> theta_k4 = theta_current;
            theta_k4[i] += k3;
            k4 = dt * f1(theta1_current[i] + l3);
            l4 = dt * f2(t + dt, theta_k4.data(), alpha, beta, i);

            // 更新原数组
            theta[i] += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            theta1[i] += (l1 + 2 * l2 + 2 * l3 + l4) / 6;
        }

        // Boundary conditions
        theta[n-1] = 0;
        theta1[n-1] = 0;
        theta[0] = 0;
        theta1[0] = 0;

        t += dt;
        cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
    }
}

额外说明 #

• 修复后，RK4每个阶段的中间状态都是完整数组，f2可正确访问theta[i]（以及后续恢复耦合项时的theta[i+1]/theta[i-1]），避免内存越界。
• 使用vector替代变长数组，既符合C++标准，又能自动管理内存，避免潜在栈溢出问题。
• 保存当前状态到theta_current和theta1_current，确保每个振子的RK4计算都基于同一时刻的初始状态，不会被其他振子的中间更新干扰。

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