龙格-库塔四阶

Runge-Kutta 4th-Order

现在我正在努力让我的 Runge-Kutta 发挥作用。实际上我必须解决一个比这个 ODE 系统复杂得多的 ODE 系统,但是 runge-kutta-function 是一样的。我正在解决这个简单的系统,以检查模拟是否收敛到解析解。

首先,如果您发现我的代码中有任何错误,我想问问您。问题是,我得到了一个解决方案,但它不同于解析解(通过 symbolab 解决)。这是我的代码:

#include <iostream>
#include <vector>
#include <fstream>

using namespace std;

double fx(double x, double y){
    return y;
}

double fy(double x, double y){
    return x + y;
}

void runge_kutta(double dt, double &x, double &y){
    vector<double> k(4), l(4);

    #define dt_half dt/2.0
    #define dt_six dt/6.0

    k.at(0)=fx(x, y);
    l.at(0)=fy(x,y);

    k.at(1)=fx(x+dt_half*k.at(0), y+dt_half*l.at(0));
    l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));

    k.at(2)=fx(x+dt_half*k.at(1), y+dt_half*l.at(1));
    l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));

    k.at(3)=fx(x+dt*k.at(2), y+dt*l.at(2));
    l.at(3)=fy(x+dt*k.at(2), y+dt*l.at(2));

    x+=dt_six*(k.at(0)+2.0*k.at(1)+2.0*k.at(2)+k.at(3));
    y+=dt_six*(l.at(0)+2.0*l.at(1)+2.0*l.at(2)+l.at(3));


}

double t=10.0;
double dt=1e-3;

double x=0.0;
double y=1.0;

ofstream out;

int main(){
    out.open("out.txt");
    out << "#t, x, y" << endl;
    for(double time=0.0;time<t;time+=dt){
        cout<<time<<endl;
        out << time << ' '
            << x << ' '
            << y << ' '
            << endl;
        runge_kutta(dt, x, y);
    }
    out.close();  
}

你能看到任何错误吗?

非常感谢您的帮助!

仔细观察这些线条

l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));

我认为他们应该使用 fy 而不是 fx。

我已经强调了 RK4 代码中的几个不一致之处

任何处理 y 的东西都应该调用 fy 并且有导数 l 任何处理 x[=22= 的东西] 有导数 k.