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$
    – eldo
    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$
    – InquisitiveInquirer
    Jun 2, 2014 at 22:00
  • $\begingroup$ oops, shame over me :) $\endgroup$
    – eldo
    Jun 2, 2014 at 22:03
  • $\begingroup$ Have a look at reference.wolfram.com/mathematica/tutorial/… $\endgroup$
    – Sascha
    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$
    – InquisitiveInquirer
    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$
    – InquisitiveInquirer
    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$
    – Holger Schmitz
    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$
    – InquisitiveInquirer
    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$
    – InquisitiveInquirer
    Jun 5, 2014 at 8:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.