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.

$\endgroup$
9
  • $\begingroup$ Maybe related: link $\endgroup$
    – eldo
    Jun 2 '14 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$ Jun 2 '14 at 22:00
  • $\begingroup$ oops, shame over me :) $\endgroup$
    – eldo
    Jun 2 '14 at 22:03
  • $\begingroup$ Have a look at reference.wolfram.com/mathematica/tutorial/… $\endgroup$
    – Sascha
    Jun 3 '14 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$ Jun 4 '14 at 9:18
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.

$\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$ Jun 4 '14 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$ Jun 4 '14 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$ Jun 4 '14 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$ Jun 5 '14 at 8:57

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .