Рунге-Кутта метод 4-го порядка кумулятивные ошибки при переборе

Question

Я пытался получить численное решение на простой задаче преследования

(движущейся мишени + ракеты с постоянным модулем скорости)

Каждая итерация мой модуль скорости уменьшается бит, добавив ошибку; И после пары сот итерационной ошибки взрывается, и скорость резко падает.

Однако это не относится к методу Эйлера (код ниже большого блока), и он появляется только при использовании метода 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); 

и так же позже. В качестве альтернативы вы можете несколько все к-значения по DT:

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

Эти два варианта еще более важным для неявных методов Рунге-Кутта