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.
Best Answer
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
+xinf1and the+yinf2look a bit weird.