为 n 维系统实现模块化 Runge-kutta 4 阶方法
Posted
技术标签:
【中文标题】为 n 维系统实现模块化 Runge-kutta 4 阶方法【英文标题】:Implementing modular Runge-kutta 4th order method for a n-dimension system 【发布时间】:2019-12-07 23:29:42 【问题描述】:我正在尝试使我的 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 并且不能作为变量进行更改,而且我不能将变量放入数组的大小中。
此外,这两个函数具有相互依赖关系,因为需要在rk4
、eq
中调用F
号码。
您能给我一些建议吗?我已经搜索并阅读了有关此的书籍,但找不到答案。这可能是一项简单的任务,但我在 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
【问题讨论】:
我非常怀疑该代码是否可以用作 RK4 方法的代码。您只得到一个订单 1 方法,而不是预期的订单 4 方法。测试它,计算一个非标量测试问题的误差,用不同步长的已知解决方案。 您正在使用 C++。您是否考虑过将 ODE 系统实现为一个导出维度和 ODE 函数的类(通过某个接口或基类)? 好的,振荡器是一个循环耦合系统。您是否还以数字方式验证了该方法的顺序?正如我所说,该方法仍然具有 1 阶,并且在第一步之后数组将填充 1 阶正确值,因此如果步长足够小,您仍然可以获得视觉上正确的结果。 对于一些包含大量样板的 C++ 代码,请查看***.com/q/49738918/3088138 的答案,或者people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html 的更简单的实现 通过一些巧妙的指针处理,您可以在循环中压缩 RK4 的阶段,使用系数数组c= 0,0.5, 0.5, 1;
和 b=1,2,2,1;
同时评估导数并更新下一个状态向量。
【参考方案1】:
使用向量函数,其中向量算术取自 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();
将向量传递给要填充的函数显然需要在 Eigen 中进行一些更高的模板考虑。
【讨论】:
以上是关于为 n 维系统实现模块化 Runge-kutta 4 阶方法的主要内容,如果未能解决你的问题,请参考以下文章
matlab实现RK45(Runge-Kutta45ode45)求解器算法
四阶龙格库塔(Runge-Kutta)求解微分方程-多种编程语言