Solving coupled 2nd order ODEs with Runge-Kutta 4

I'm having a hard time figuring out how coupled 2nd order ODEs should be solved with the RK4 method. This is the system I'm given:

$x'' = f(t, x, y, x', y')$
$y'' = g(t, x, y, x', y')$

I'll use the notation $u = x',\ w = y'$, thus $u' = x'',\ w' = y''$.
I am also given the initial conditions $x_0 = x(0),\ y_0 = y(0)$ and $u_0 = u(x_0),\ w_0 = w(y_0)$.

For the specific initial values that I have (which are not that relevant here), when I plot $x$ vs. $y$ I should get a closed loop. However, what I get is this spiral below:

I am certain that I should get a closed loop, so I must be applying RK4 incorrectly. Here's how I did it (you can view the MATLAB function here if you want).

Firstly, here's how I do the final approximations (substitute $a$ for $x$, $y$, $u$, or $w$):

$a_{i+1} = a_i + \frac{1}{6}(k_{a1} + 2k_{a2} + 2k_{a3} + k_{a4})$

To calculate the intermediate approximations (the $k$s), I do this:

$k_{x1} = u_i \cdot \Delta t$
$k_{y1} = w_i \cdot \Delta t$
$k_{u1} = f(t, x_i, y_i, u_i, w_i) \cdot \Delta t$
$k_{w1} = g(t, x_i, y_i, u_i, w_i) \cdot \Delta t$

$k_{x2} = (u_i + \frac{k_{u1}}{2}) \cdot \Delta t$
$k_{y2} = (w_i + \frac{k_{w1}}{2}) \cdot \Delta t$
$k_{u2} = f(t+\frac{\Delta t}{2}, x_i+\frac{k_{x1}}{2}, y_i+\frac{k_{y1}}{2}, u_i+\frac{k_{u1}}{2}, w_i+\frac{k_{w1}}{2}) \cdot \Delta t$
$k_{w2} = g(t+\frac{\Delta t}{2}, x_i+\frac{k_{x1}}{2}, y_i+\frac{k_{y1}}{2}, u_i+\frac{k_{u1}}{2}, w_i+\frac{k_{w1}}{2}) \cdot \Delta t$

$k_{x3} = (u_i + \frac{k_{u2}}{2}) \cdot \Delta t$
$k_{y3} = (w_i + \frac{k_{w2}}{2}) \cdot \Delta t$
$k_{u3} = f(t+\frac{\Delta t}{2}, x_i+\frac{k_{x2}}{2}, y_i+\frac{k_{y2}}{2}, u_i+\frac{k_{u2}}{2}, w_i+\frac{k_{w2}}{2}) \cdot \Delta t$
$k_{w3} = g(t+\frac{\Delta t}{2}, x_i+\frac{k_{x2}}{2}, y_i+\frac{k_{y2}}{2}, u_i+\frac{k_{u2}}{2}, w_i+\frac{k_{w2}}{2}) \cdot \Delta t$

And finally...

$k_{x4} = (u_i + k_{u3}) \cdot \Delta t$
$k_{y4} = (w_i + k_{w3}) \cdot \Delta t$
$k_{u4} = f(t+\Delta t, x_i+k_{x3}, y_i+k_{y3}, u_i+k_{u3}, w_i+k_{w3}) \cdot \Delta t$
$k_{w4} = g(t+\Delta t, x_i+k_(x3), y_i+K_{y3}, u_i+k_{u3}, w_i+k_{w3}) \cdot \Delta t$

I can't figure out where the mistake is. Can anyone please help me? Btw, you'll notice that in the MATLAB code, the functions $f$ and $g$ don't have a $t$ parameter — that's because it does not contribute to the value of the functions so I left it out.

Solution 1:

When I try to solve the ODE in your Matlab file with the built-in solver ode45, I get a very similar picture. So I think your implementation of RK4 is fine. I don't know what makes you that certain that you should get closed loops, but I'd suggest you take a good look at the ODEs and make sure that these are the correct equations. To me, the +x in f1 and the +y in f2 look a bit weird.

Solution 2:

I have done this before and with some simplifications I arrived at the following 2nd order scheme for RK4:

Function Derivatives (where ${\rm h} = \Delta t$)

K0x = f(t,x,y,xp,yp)            
K0y = g(t,x,y,xp,yp)
Q1x = xp + (h/2)*K0x            
Q1y = yp + (h/2)*K0y
K1x = f(t+h/2,x+(h/2)*xp,y+(h/2)*yp,Q1x,Q1y)    
K1y = g(t+h/2,x+(h/2)*xp,y+(h/2)*yp,Q1x,Q1y)
Q2x = xp + (h/2)*K1x           
Q2y = yp + (h/2)*K1y
K2x = f(t+h/2,x+(h/2)*Q1x,y+(h/2)*Q1y,Q2x,Q2y)   
K2y = g(t+h/2,x+(h/2)*Q1x,y+(h/2)*Q1y,Q2x,Q2y)
Q3x = xp + h*K2x         
Q3y = yp + h*K2y
K3x = f(t+h, x+h*Q2x, y+h*Q2y, Q3x,Q3y)      
K3y = g(t+h, x+h*Q2x, y+h*Q2y, Q3x,Q3y)

Simulation Step

x  -> x  + h*(xp+(h/6)*(K0x + K1x + K2x))
y  -> y  + h*(yp+(h/6)*(K0y + K1y + K2y))
xp -> xp + (h/6)*(K0x + 2*K1x + 2*K2x + K3x)
yp -> yp + (h/6)*(K0y + 2*K1y + 2*K2y + K3y)

Here is a link to a similar answer in SO.