c++ - 为 n 维系统实现模块化 Runge-kutta 4 阶方法

Question

我正在尝试使我的 runge-kutta 四阶代码模块化。我不想每次使用时都必须编写和声明代码，但是 在 .hpp 和 .cpp 文件中声明它以单独使用 .但我遇到了一些问题。一般来说，我想求解一个 n 维方程组。为此，我使用了两个函数:一个用于方程组，另一个用于 runge-kutta 方法，如下所示:

double F(double t, double x[], int eq)
{
    // System equations
    if      (eq == 0) { return (x[1]); }
    else if (eq == 1) { return (gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]); }
    else if (eq == 2) { return (-kappa * x[1] - phi * x[2]); }
    else { return 0; }
}

void rk4(double &t, double x[], double step)
{
    double x_temp1[sistvar], x_temp2[sistvar], x_temp3[sistvar];        
    double k1[sistvar], k2[sistvar], k3[sistvar], k4[sistvar];      

    int j;                                                      

    for (j = 0; j < sistvar; j++) 
    {
        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
    }
    for (j = 0; j < sistvar; j++)
    {
        k4[j] = step * F(t + step, x_temp3, j);
    }
    for (j = 0; j < sistvar; j++) 
    {
        x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
    }

    t += step;
}

上面的代码有效并且已经过验证。但是它有一些依赖关系，因为它使用一些全局变量来工作:
gama , OMEGA , zeta , alpha , beta , chi , kappa和 phi是我想从 .txt 文件中读取的全局变量。我已经设法做到了，但是只在一个包含所有代码的 .cpp 文件中。

另外， sistvar是系统维度，也是一个全局变量。我正在尝试将其作为参数输入 F .但是它的编写方式似乎给出了错误，因为 sistvar 是一个 const 并且不能作为变量进行更改，而且我不能将变量放入数组的大小中。

此外，这两个函数具有相互依赖关系，就像调用 F 时一样。里面 rk4 , eq需要号码。

你能给我一些建议吗？我已经搜索并阅读了有关此的书籍，但找不到答案。这可能是一项简单的任务，但我在 c/c++ 编程语言方面相对较新。

提前致谢!

* 已编辑(尝试使用 std::vector 实现)*

double F(double t, std::vector<double> x, int eq)
{
    // System Equations
    if (eq == 0) { return (x[1]); }
    else if (eq == 1) { return (gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]); }
    else if (eq == 2) { return (-kappa * x[1] - phi * x[2]); }
    else { return 0; }
}

double rk4(double &t, std::vector<double> &x, double step, const int dim)
{
    std::vector<double> x_temp1(dim), x_temp2(dim), x_temp3(dim);
    std::vector<double> k1(dim), k2(dim), k3(dim), k4(dim);

    int j;

    for (j = 0; j < dim; j++) {
        x_temp1[j] = x[j] + 0.5*(k1[j] = step * F(t, x, j));
    }
    for (j = 0; j < dim; j++) {
        x_temp2[j] = x[j] + 0.5*(k2[j] = step * F(t + 0.5 * step, x_temp1, j));
    }
    for (j = 0; j < dim; j++) {
        x_temp3[j] = x[j] + (k3[j] = step * F(t + 0.5 * step, x_temp2, j));
    }
    for (j = 0; j < dim; j++) {
        k4[j] = step * F(t + step, x_temp3, j);
    }

    for (j = 0; j < dim; j++) {
        x[j] += (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
    }

    t += step;

    for (j = 0; j < dim; j++) {
        return x[j];
    }
}

vector                 array

2.434 s   |        |   0.859 s
2.443 s   |        |   0.845 s
2.314 s   |        |   0.883 s
2.418 s   |        |   0.884 s
2.505 s   |        |   0.852 s
2.428 s   |        |   0.923 s
2.097 s   |        |   0.814 s
2.266 s   |        |   0.922 s
2.133 s   |        |   0.954 s
2.266 s   |        |   0.868 s
_______                _______
average = 2.330 s      average = 0.880 s

Answer 1

使用 vector 函数，其中 vector 算术取自 Eigen3

#include <eigen3/Eigen/Dense>
using namespace Eigen;

与问题中讨论的部分相同的部分可能看起来像(受 function pointer with Eigen 启发)

VectorXd Func(const double t, const VectorXd& x)
{ // equations for solving simple harmonic oscillator
    Vector3d dxdt;
    dxdt[0] = x[1];
    dxdt[1] = gama * sin(OMEGA*t) - zeta * x[1] - alpha * x[0] - beta * pow(x[0], 3) - chi * x[2]; 
    dxdt[2] = -kappa * x[1] - phi * x[2]; 
    return dxdt;
}

MatrixXd RK4(VectorXd Func(double t, const VectorXd& y), const Ref<const VectorXd>& y0, double t, double h, int step_num)
{
    MatrixXd y(y0.rows(), step_num );
    VectorXd k1, k2, k3, k4;
    y.col(0) = y0;

    for (int i=1; i<step_num; i++){
        k1 = Func(t, y.col(i-1));
        k2 = Func(t+0.5*h, y.col(i-1)+0.5*h*k1);
        k3 = Func(t+0.5*h, y.col(i-1)+0.5*h*k2);
        k4 = Func(t+h, y.col(i-1)+h*k3);
        y.col(i) = y.col(i-1) + (k1 + 2*k2 + 2*k3 + k4)*h/6;
        t = t+h;
    }
    return y.transpose();
}

将 vector 传递给要填充的函数显然需要在 Eigen 中进行一些更高的模板考虑。