# Runge-Kutta to find maximum of system of ODEs

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.

• The code as written determines the value of t for which x[[1]] vanishes, and it does so successfully. By the way, you may wish to delete your earlier comments, which no longer are relevant. Also, you may wish to use NDSolve and Maximize instead of your own routines. – bbgodfrey Nov 21 '16 at 3:09
• I'm trying to use my own routines to better understand numerical analysis. Do you see a way to use this code to maximize x[[1]]? – glassseahorse Nov 21 '16 at 3:18

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.

• Looks good! thank you. How did you calculate the temporal derivatives? – glassseahorse Nov 21 '16 at 13:05
• First, I differentiated (x1[t]^2 + x3[t]^2) once and twice using D. Then, I replaced x1[t] by x[[1]], x1'[t] by x[[2]], and x1''[t] by the second term in the definition of g, and did the same for the x3 terms in the derivatives. Hope this helps. – bbgodfrey Nov 21 '16 at 13:09
• Why did you differentiate (x1[t]^2 + x3[t]^2)  and not Sqrt[(x1[t]^2 + x3[t]^2)]? – glassseahorse Nov 21 '16 at 14:02
• They both have the same values of t for extremes, and not taking Sqrt results in a simpler expression. Since you apparently are doing all this as a learning exercise, try using the Sqrt` expression and see what happens. – bbgodfrey Nov 21 '16 at 14:11