# Runge-Kutta 2nd Order ODE Solver

Suppose I have a 2nd order ODE of the form y''(t) = 1/y with y(0) = 0 and y'(0) = 10, and want to solve it using a Runge-Kutta solver. I've read that we need to convert the 2nd order ODE into two 1st order ODEs, but I'm having trouble doing that at the moment and am hoping someone here might be able to help. This is my code thus far:

Remove["Global*"]
(*dy/dt=*)f[t_, y_] := 1/y;
(*d^2y/dt^2=*)g[t_, y_, yd_] :=???;
t[0] = 0;
y[0] = 0;
yd[0] = 10;
tmax = 1000;
h = 0.01;

Do[
{t[n] = t[0] + h n,

k1 = h f[t[n], y[n], yd[n]];
l1 = h g[t[n], y[n], yd[n]];

k2 = h f[t[n] + h/2, y[n] +  k1/2, yd[n] + l1/2];
l2 = h g[t[n] + h/2, y[n] +  k1/2, yd[n] + l1/2];

k3 = h f[t[n] + h/2, y[n] + k2/2, yd[n] + l2/2];
l3 = h g[t[n] + h/2, y[n] + k2/2, yd[n] + l2/2];

k4 = h f[t[n] + h, y[n] + k3, yd[n] + l3];
l4 = h g[t[n] + h, y[n] + k3, yd[n] + l3];

y[n + 1] = y[n] + 1/6 (k1 + 2 k2 + 2 k3 + k4);
yd[n + 1] = yd[n] + 1/6 (l1 + 2 l2 + 2 l3 + l4);

}, {n, 0, tmax}]


As you can see by the question marks for the function g[t_,y_,yd_], I don't know how I can set it in such a way that y''(t) = 1/y. Do I feed the results of y[n+1] into g when running the algorithm? Any help would be appreciated.

-
Maybe related: link – eldo Jun 2 '14 at 21:46
That's the question I asked yesterday :) I want to move on from a plain system of 1st order ODEs to a 2nd order ODE using a system of equations. – user7388 Jun 2 '14 at 22:00
oops, shame over me :) – eldo Jun 2 '14 at 22:03
Have a look at reference.wolfram.com/mathematica/tutorial/… – Sascha Jun 3 '14 at 19:33
Thanks Sascha, but unfortunately that doesn't show explicitly how to convert a 2nd order ODE to two 1st order ODEs using RK4. It seems to give more of an overview of possible solvers used by NDSolve. – user7388 Jun 4 '14 at 9:18

You can specify the functions like this

(*dy/dt=*)f[t_, y_, yd_] := yd;
(*d^2y/dt^2=*)g[t_, y_, yd_] := 1/y;
t[0] = 0;
y[0] = 1;
yd[0] = 10;
tmax = 1000;
h = 0.01;
Do[{t[n] = t[0] + h n, k1 = h f[t[n], y[n], yd[n]];
l1 = h g[t[n], y[n], yd[n]];
k2 = h f[t[n] + h/2, y[n] + k1/2, yd[n] + l1/2];
l2 = h g[t[n] + h/2, y[n] + k1/2, yd[n] + l1/2];
k3 = h f[t[n] + h/2, y[n] + k2/2, yd[n] + l2/2];
l3 = h g[t[n] + h/2, y[n] + k2/2, yd[n] + l2/2];
k4 = h f[t[n] + h, y[n] + k3, yd[n] + l3];
l4 = h g[t[n] + h, y[n] + k3, yd[n] + l3];
y[n + 1] = y[n] + 1/6 (k1 + 2 k2 + 2 k3 + k4);
yd[n + 1] = yd[n] + 1/6 (l1 + 2 l2 + 2 l3 + l4);}, {n, 0, tmax}]


The function f is the derivative of y and therefore equal to yd. The function g is the derivative of yd which means it is the second derivative of the function you are looking for. Here you can specify the right hand side of the equation y''=1/y.

Also, you shouldn't specify y(0) = 0 in your initial conditions, because then 1/y is not defined.

-
Unfortunately that doesn't seem to do the trick. I've also tried setting f[t_, y_, yd_] := y; and that didn't work either. – user7388 Jun 4 '14 at 14:26
I have edited the answer to include the complete code. This works fine when I run it. Note that you have to set y[0]=1;. – Holger Schmitz Jun 4 '14 at 14:56
If you try g[t_, y_, yd_] := -y (Hooke's Law) and then plot the output using Soln1 = Table[{t[i], y[i]}, {i, 0, tmax}];, Soln2 = Table[{t[i], yd[i]}, {i, 0, tmax}]; ListLinePlot[Soln1], and ListLinePlot[Soln2] and compare this to the results from NDSolve using y'[t] == -y[t] and y''[t] == -y[t]` you'll see that something is slightly awry. – user7388 Jun 4 '14 at 15:27
Apologies! Having had a closer look, I see that you are in fact correct. Thank you Holger, I appreciate your help. – user7388 Jun 5 '14 at 8:57