c++ - Runge-Kutta 四阶法迭代时的累积误差

Question

我试图获得一个简单追逐问题的数值解

(动靶+定速模块火箭)

每次迭代我的速度模块都会减少一点，将误差加起来；在几百次迭代之后，错误爆发，速度急剧下降。

但是，使用 Euler 方法(大块下方的代码)时，情况并非如此，它仅在使用 RK4 方法时弹出。

我不确定错误在哪里以及为什么会发生，所以欢迎任何意见

#include <fstream>
#include <iostream>
#include <vector>
#include <cmath>
#include <cstring>
#define vx(t,x,y) n*V*((t)*(V)-(x))/pow(((t)*(V)-(x))*((t)*(V)-(x))+((h)-(y))*((h)-(y)),0.5)
#define vy(t,y,x) n*V*((h)-(y))/pow(((t)*(V)-(x))*((t)*(V)-(x))+((h)-(y))*((h)-(y)),0.5)
using namespace std;
class Vector {
public:
    double x,y;
    Vector(double xx, double yy):x(xx), y(yy){};
    virtual ~Vector(){}
    Vector operator-() {return Vector(-x,-y);};
    friend Vector operator-(const Vector &, const Vector &); 
    friend Vector operator+(const Vector &, const Vector &); 
    Vector operator*(double l){return Vector(x*l,y*l);};
    friend Vector operator*(double, const Vector &);
    Vector operator/(double l){return Vector(x/l,y/l);};
    void operator+=(const Vector & v ){ x+=v.x; y+=v.y;};
    void operator-=(const Vector & v ){ x-=v.x; y-=v.y;};
    void operator/=(const Vector & v ){ x/=v.x; y/=v.y;};
    friend ostream & operator<<(ostream & os,const Vector & v){os<<"("<<v.x<<", "<<v.y<<")";return os;};
    double norm() {return sqrt(x*x+y*y);};
};

Vector operator-(const Vector & v1, const Vector & v2){
    return Vector(v1.x-v2.x,v1.y-v2.y);
}
Vector operator+(const Vector & v1, const Vector & v2){
    return Vector(v1.x+v2.x,v1.y+v2.y);
}
Vector operator*(double l, const Vector & v){
    return Vector(v.x*l,v.y*l);
}

int main() {
    Vector posP(0,0);
    double V=100.,t = 0,dt = pow(10.,-2),vx,vy,h=1000.,x,y,n=2.,v;
    double kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4;
    Vector posE(0,h);
    FILE *fp;
    fp = fopen("/Users/Philipp/Desktop/a.dat","w");
    while(posP.y<(h)){
        posE.x=posE.x+V*dt;
        x=posP.x;y=posP.y;
        kx1 = vx(t,x,y);
        ky1 = vy(t,y,x);
        kx2 = vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);
        ky2 = vy(t+dt/2.0,y+ky1/2.0,x+kx1/2.0);
        kx3 = vx(t+dt/2.0,x+kx2/2.0,y+ky2/2.0);
        ky3 = vy(t+dt/2.0,y+ky2/2.0,x+kx2/2.0);
        kx4 = vx(t+dt,x+kx3,y+ky3);
        ky4 = vy(t+dt,y+ky3,x+kx3);
        posP.x = posP.x + dt*((kx1 + 2.0*(kx2+kx3) + kx4)/6.0);
        posP.y = posP.y + dt*((ky1 + 2.0*(ky2+ky3) + ky4)/6.0);
        v=sqrt(((kx1 + 2.0*(kx2+kx3) + kx4)/(6.0))*((kx1 + 2.0*(kx2+kx3) + kx4)/(6.0))+((ky1 + 2.0*(ky2+ky3) + ky4)/(6.0))*((ky1 + 2.0*(ky2+ky3) + ky4)/(6.0)));
        t+=dt;
        if ((posE-posP).norm()<1) break;
        fprintf(fp,"%lf %lf %lf %lf \n",posP.x, posP.y, v, t);
    }
    fclose(fp);
    return 0;
}

欧拉方法

//Euler cycle
    while(posP.y<(h)) {
        posE.x=posE.x+V*dt;
        x=posP.x;y=posP.y;
        vx=vx(t,x,y);
        vy=vy(t,y,x);
        posP.x=posP.x+vx*dt;
        posP.y=posP.y+vy*dt;
        t+=dt;
        if ((posE-posP).norm()<0.1) break;
        fprintf(fp,"%lf %lf %lf \n",posP.x, posP.y,vx*vx+vy*vy, t);

//速度模块在第三列，你可以看到一切都是 200，而不是 RK4 的情况，即使是第一次迭代它也会下降到 ~199.99985

    }

Answer 1

你用

kx2 = vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);

但应该是:

kx2 = vx(t+dt/2.0,x+kx1/2.0*dt,y+ky1/2.0*dt);

以后也类似。或者，您可以将所有 k 值乘以 dt:

kx2 = dt*vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);

这两个变体对于隐式 Runge-Kutta 方法更为重要