Runge-Kutta to find maximum of system of ODEs

Question

I'm trying to understand a code for finding the minimum and maximum for a system of ODEs. Here is the code:

RK[f_, t_, x_, h_, n_] := 
 Module[{i, s = t, k1, k2, k3, k4, y = x, h2 = h/2},
  For[i = 1, i <= n, i++, k1 = h f[s, y];
   s += h2;
   k2 = h f[s, y + k1/2];
   k3 = h f[s, y + k2/2];
   s += h2;
   k4 = h f[s, y + k3];
   y += (k1 + 2 (k2 + k3) + k4)/6];
  y]

g[t_, x_] := {x[[2]], -x[[1]]/(x[[1]]^2 + x[[3]]^2)^(3/2), 
   x[[4]], -x[[3]]/(x[[1]]^2 + x[[3]]^2)^(3/2)};

lst = {};
dt = 0.1;
n = 100;

x = {0, -.5, 1., -.5};

Do[lst = Append[lst, {t, (x[[1]]^2 + x[[3]]^2)^(1/2)}];
x = RK[g, t, x, dt/n, n]
, {t, 0, 10, dt}]

x0 = {0, -.5, 1., -.5};
t = .9; (*Initial guess*)
Do[x = RK[g, 0, x0, t/n, n];
 t = t - x[[1]]/x[[2]];
 Print[t, " ", x], {i, 0, 5}]

ListPlot[lst, Joined -> True]

The Range-Kutta is fourth order and I understand how it works. The system of equations is modeling a 2-body system under a gravitational force. I have the vector x = {x, x', y, y'}. The system has the initial values of x(0) = 0, x'(0) = -0.5, y(0) = 1, y'(0) = -.5. I have plotted the values of Sqrt[x^2 + y^2] to plot the distance of the object over time and have the following:

I'm trying to find the minimum and maximum distances using Newton's method, as seen in the last few lines of the code. The code currently works for finding the minimum distance, which it calculates to be at t = 0.90715, which is a distance of 0.14286, which agrees with the plot. My problem is finding the maximum distance now using this code. I thought I would just need to change the initial guess to 2.6, which is where the maximum occurs on the plot. This does not work in the code and the newton method resolves at t = 3.4115 which is neither a maximum nor minimum. I believe what the code is doing is trying different values of t because the Range-Kutta will give the answer to the system at that time t. I'm not sure why it is not calculating the maximum distance.

Answer 1

The question uses Newton's method based on t = t - x[[1]]/x[[2]]. Conceptually, the correction term is f/f'. Hence, the question actually is solving for x[[1]] == 0. To obtain the extremes of (x[[1]]^2 + x[[3]]^2)^(1/2), we must solve for the temporal derivative of x[[1]]^2 + x[[3]]^2 set equal to zero. The corresponding Newton's method then becomes

t = t -  (x[[1]] x[[2]] + x[[3]] x[[4]])/(x[[2]]^2 + x[[4]]^2 
    - x[[1]]^2 /(x[[1]]^2 + x[[3]]^2)^(3/2) - x[[3]]^2 /(x[[1]]^2 + x[[3]]^2)^(3/2))

Using it with no other changes in the question's code and an initial guess of 0.9 yields t == 0.894415, and with an initial guess of 2.6 yields t == 2.57204. The values of t agree well with the plot in the question.

These results can, of course, be obtained in a straightforward manner using higher level Mathematica functions:

Clear[t];
s = Flatten@NDSolve[{x1''[t] == -x1[t]/(x1[t]^2 + x3[t]^2)^(3/2), 
    x3''[t] == -x3[t]/(x1[t]^2 + x3[t]^2)^(3/2), x1[0] == 0, 
    x1'[0] == -1/2, x3[0] == 1, x3'[0] == -1/2}, {x1[t], x3[t]}, {t, 0, 10}];

Plot[Sqrt[x1[t]^2 + x3[t]^2] /. s, {t, 0, 10}, AxesOrigin -> {0, 0}]

Last@FindMinimum[Evaluate[(x1[t]^2 + x3[t]^2) /. s], {t, 0.9}]
(* {t -> 0.894414} *)

Last@FindMaximum[Evaluate[(x1[t]^2 + x3[t]^2) /. s], {t, 2.6}]
(* {t -> 2.60448} *)

The results obtained from the question's modified code agree reasonably well with the more accurate results just obtained, and better agreement can be obtain by using larger n, say n == 1000.