c++ - Runge-Kutta 4th order method cumulative errors when iterating -

i trying numerical solution on simple chase problem

(moving target+rocket constant speed module)

every iteration speed module decreases bit, adding error; , after couple of hundred iteration error blows , speed drops dramatically.

however, not case euler method(code below big block) , popping when using rk4 method.

i not sure error , why it's happening, input appreciated

#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 method

//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); 

//speed module in third column, can see 200, not case rk4, first iteration drops ~199.99985

    } 

you use

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

but should be:

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

and later on. alternatively multiple k-values dt:

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

those 2 variants more important implicit runge-kutta methods