2
$\begingroup$

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.

Improve this question
$\endgroup$
9
  • $\begingroup$ Maybe related: link $\endgroup$ Commented Jun 2, 2014 at 21:46
  • $\begingroup$ 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. $\endgroup$ Commented Jun 2, 2014 at 22:00
  • $\begingroup$ oops, shame over me :) $\endgroup$ Commented Jun 2, 2014 at 22:03
  • $\begingroup$ Have a look at reference.wolfram.com/mathematica/tutorial/… $\endgroup$ Commented Jun 3, 2014 at 19:33
  • $\begingroup$ 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. $\endgroup$ Commented Jun 4, 2014 at 9:18

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$ Commented Jun 4, 2014 at 14:26
  • $\begingroup$ 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;. $\endgroup$ Commented Jun 4, 2014 at 14:56
  • $\begingroup$ 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. $\endgroup$ Commented Jun 4, 2014 at 15:27
  • $\begingroup$ Apologies! Having had a closer look, I see that you are in fact correct. Thank you Holger, I appreciate your help. $\endgroup$ Commented Jun 5, 2014 at 8:57

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.