我正在尝试使我的 runge-kutta 四阶代码模块化。我不想每次使用它时都必须编写和声明代码,但是在.hpp和.cpp文件中声明它以分别使用它。但我遇到了一些问题。一般来说,我想求解 n 维方程组。为此,我使用两个函数:一个用于方程组,另一个用于龙格-库塔方法,如下所示:
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
and phi
是我想从 .txt 文件读取的全局变量。我已经设法做到这一点,但只是在包含所有代码的单个 .cpp 文件中。
Also, sistvar
是系统维度,也是全局变量。我试图将其作为参数输入F
。但它的编写方式似乎给出了错误,因为 sistvar 是一个常量,不能作为变量进行更改,而且我不能将变量放入数组的大小中。
此外,这两个函数具有相互依赖性,就像调用时一样F
inside 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