Runge-Kutta 4th order method cumulative errors when iterating

Question

I was trying to get a numerical solution on a simple chase problem

(moving target+rocket with constant speed module)

Every iteration my speed module decreases a bit, adding up the error; And after a couple of hundred iteration error blows up and speed drops dramatically.

However, this is not the case with Euler method(code below big block) and it is only popping up when using RK4 method.

I am not sure where the error is and why it's happening, so any input is 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 is in third column, as you can see everything is 200, not the case with RK4, where even first iteration it drops to ~199.99985

Floating point values on computers tend to have rounding errors. These rounding errors compound. See [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken). — Some programmer dude, Nov 17 '16 at 14:34
@Someprogrammerdude then why is Euler method working even with small step size? I still think that there's a problem with my RK4 method — Philipp Balakin, Nov 17 '16 at 14:37
Please use an inline function for vx and vy. All those brackets make it completely unreadable. — Martin Bonner supports Monica, Nov 17 '16 at 14:47
Also, you say "speed module is in third column" - what third column? — Martin Bonner supports Monica, Nov 17 '16 at 14:48
@MartinBonner in the output file; Regarding vx&vy: https://puu.sh/slliu/808a54c15f.png — Philipp Balakin, Nov 17 '16 at 15:12
Why do you define operators for the vector class if you then do not use them to pass the state to the derivative function? Instead of two `vx,vy` functions you could have one that takes a `Vector` argument and returns a `Vector`. — Lutz Lehmann, Nov 17 '16 at 16:09

Hans Olsson · Accepted Answer · 2016-11-17T15:13:52.523

1

You use

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

But it should be:

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

and similarly later on. Alternatively you could multiple all k-values by dt:

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

Those two variants are even more important for implicit Runge-Kutta methods

edited Nov 17 '16 at 15:13

answered Nov 17 '16 at 15:10

Hans Olsson

11,123
15
38

Does it really matter if I use it there instead of in `posP.x = posP.x + dt*((kx1 + 2.0*(kx2+kx3) + kx4)/6.0);`? – Philipp Balakin Nov 17 '16 at 15:13
1

Yes. You always have to multiply a velocity by dt when you add it to x – Hans Olsson Nov 17 '16 at 15:15
Well it is multiplied later on, in the piece of code I just sent to you – Philipp Balakin Nov 17 '16 at 15:16
1

But you have to multiply by dt consistently. See a standard reference such as Hairer, Noersett, Wanner "Solving Ordinary Differential Equations I", or https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods – Hans Olsson Nov 17 '16 at 15:31

Runge-Kutta 4th order method cumulative errors when iterating

1 Answers1