# 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, 2014 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. Jun 2, 2014 at 22:00
• oops, shame over me :)
– eldo
Jun 2, 2014 at 22:03
• Have a look at reference.wolfram.com/mathematica/tutorial/… Jun 3, 2014 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. Jun 4, 2014 at 9:18

## 1 Answer

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. Jun 4, 2014 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;. Jun 4, 2014 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. Jun 4, 2014 at 15:27
• Apologies! Having had a closer look, I see that you are in fact correct. Thank you Holger, I appreciate your help. Jun 5, 2014 at 8:57