Exploding Runge Kutta Method

Question

I've been attempting to build a Runge Kutta fourth order integrator to model simple projectile motion. My code is as follows

double rc4(double initState, double (*eqn)(double,double),double now,double dt)
{
        double k1 = eqn(initState,now);
        double k2 = eqn(initState + k1*dt/2.0,now + dt/2.0);
        double k3 = eqn(initState + k2*dt/2.0,now + dt/2.0);
        double k4 = eqn(initState + k3*dt, now + dt);

        return initState + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4);
}

This is called within a while loop

while (time <= duration && yPos >=0)
                {

                        xPos = updatePosX(xPos,vx,timeStep);
                        yPos = updatePosY(yPos,vy,timeStep);


                        vx = rc4(vx,updateVelX,time,timeStep);
                        vy = rc4(vy,updateVelY,time,timeStep);

                        cout << "x Pos: " << xPos <<"\t y Pos: " << yPos << endl;

                        time+=timeStep;

                        myFile << xPos << "  " << yPos << "  " << vx << "  " << vy << endl;

                }

However, contrary to what should happen my results simply blow up. What's going on here?

Did you mean to call `rc4(updateVelX, vx, time, timeStep)` instead of `rc4(vx, updateVelX, time, timeStep) (notice the inversion of the first 2 arguments)? — SleuthEye, Jan 19 '15 at 03:06
I have attempted to fix the code as shown above, and am still running into the same issue. — user3277807, Jan 19 '15 at 03:35
What do you mean by it "blows up"? Infinite loop? There isn't enough code present to know what's going on. Where is duration initialized? The updatePos functions? — TriHard8, Jan 19 '15 at 03:45
Duration is initialized with input. I mean the update positron functions blow up. For full code please see.http://stackoverflow.com/questions/28013383/implementation-of-runge-kutta-fourth-order-in-c/28013421#28013421 — user3277807, Jan 19 '15 at 04:05

Lutz Lehmann · Answer 1 · 2015-01-29T14:56:40.617

Your rk4 code looks right. But only for scalar differential equations.

What you most certainly have is a system of coupled differential equations in a dimension greater than 1. Here you have to apply the integration method in its vector form. That is, x,y,vx,vy are combined into a 4 dimensional (phase) state vector and the system function is vector valued, k1,...k4 are vectors etc.

As an advanced note, time <= duration is sensible to rounding errors accumulated in the repetitions of time+=timeStep;. Better use time <= duration-timeStep/2 to have time at the end of the loop close to duration.

Reading the code on the closed previous question I see that you have problems with the idea of a differential equation. You should not use the result of the Euler step as acceleration in the RK4 implementation. The system for ballistic motion without air friction is

dotx = vx
doty = vy
dotvx = 0
dotvy = -g

which you would have to implement in vector form as something like

eqn(t, [x,y,vx,vy]) // where X = array of double of dimension 4
   { return [vx,vy,0,-g]; }

Exploding Runge Kutta Method

1 Answers1