Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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:

Spiral

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.

share|improve this question
add comment

2 Answers

up vote 1 down vote accepted

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.

share|improve this answer
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.