# Could anyone help me to compute the following code?

The following code works good when the value of " r < 21 ". But when r > 20, I get the following error message

General::ovfl: Overflow occurred in computation.
General::stop: Further output of General::ovfl will be suppressed during this calculation.



This is the code which I need to compute for r > 20,

is = .85 {450, 450};
ip = {{60, 10}, {20, 10}};
f[t_, x_, y_, z_] := -x + y;
g[t_, x_, y_, z_] := -x - 0.01 y - y z;
p[t_, x_, y_, z_] := -0.5 z + x^2 - r;
r = 21;
t[0] = 0;
x[0] = 0.1;
y[0] = 0.1;
z[0] = 0.1;
a[0] = 0;
b[0] = -0.1;
c[0] = -0.1;
d[0] = 0.1;
tmax = 20000;
amax = 20000;
h = 0.01;
Do[{t[n] = t[0] + h n, k1 = h f[t[n], x[n], y[n], z[n]];
l1 = h g[t[n], x[n], y[n], z[n]];
m1 = h p[t[n], x[n], y[n], z[n]];
k2 = h f[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
l2 = h g[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
m2 = h p[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
k3 = h f[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
l3 = h g[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
m3 = h p[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
k4 = h f[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
l4 = h g[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
m4 = h p[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
x[n + 1] = x[n] + 1/6 (k1 + 2 k2 + 2 k3 + k4);
y[n + 1] = y[n] + 1/6 (l1 + 2 l2 + 2 l3 + l4);
z[n + 1] = z[n] + 1/6 (m1 + 2 m2 + 2 m3 + m4);}, {n, 0, tmax}]
T2 = Table[{x[i], z[i]}, {i, 0, tmax}];
pp1 = ListLinePlot[T2, PlotStyle -> {Blur, Thin}]



• @Bill Thank you for your suggestion. I gave all the input values as rationals instead of decimals but the code did not compute even after two hours. Still it is displaying as running. Nov 8, 2020 at 9:25
• Ok, i will wait to get the output. No can't reduce the number of iterations. Nov 8, 2020 at 10:19
• Did you ever obtain a solution by for r = 21 by waiting to get the output? Nov 11, 2020 at 3:38

The r = 21 solution in the question becomes exceeding large at abut n = 895, and this has nothing to do with lack of Precision. I determined this as follows.

r = 21; h = 1/100;
f[t_, x_, y_, z_] := -x + y;
g[t_, x_, y_, z_] := -x - y/100 - y z;
p[t_, x_, y_, z_] := -z/2 + x^2 - r;


Some savings in computation time can be obtained by

t[n] = h n;
k1 = h f[t[n], x[n], y[n], z[n]];
l1 = h g[t[n], x[n], y[n], z[n]];
m1 = h p[t[n], x[n], y[n], z[n]];
k2 = h f[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
l2 = h g[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
m2 = h p[t[n] + h/2, x[n] + k1/2, y[n] + l1/2, z[n] + m1/2];
k3 = h f[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
l3 = h g[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
m3 = h p[t[n] + h/2, x[n] + k2/2, y[n] + l2/2, z[n] + m2/2];
k4 = h f[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
l4 = h g[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
m4 = h p[t[n] + h, x[n] + k3, y[n] + l3, z[n] + m3];
xp = Simplify[x[n] + 1/6 (k1 + 2 k2 + 2 k3 + k4)]
yp = Simplify[y[n] + 1/6 (l1 + 2 l2 + 2 l3 + l4)]
zp = Simplify[z[n] + 1/6 (m1 + 2 m2 + 2 m3 + m4)]


The expressions xp, yp, and zp are not reproduced here because they are fairly large. The structure of the question is that of a recurrence relation, which can be evaluated by

tab = RecurrenceTable[{x[n + 1] == xp, y[n + 1] == yp, z[n + 1] == zp,
x[0] == .1100000, y[0] == .1100000, z[0] == .1100000}, {x, y, z}, {n, 1, 895}];
N@tab[[-1]]
(* {-1.509523880701748*10^10286856148010, 1.684787554556877*10^17442606938024,
4.184274805718905*10^12133439661425} *)


which requires just less than two minutes to complete. A plot of the solution is

ListLinePlot[Transpose@tab, PlotRange -> {-400, 400}, ImageSize -> Large,
LabelStyle -> {15, Bold, Black}, PlotLegends -> Placed[{x, y, z}, {.3, .2}]]


These curves are insensitive to the Precision of {x[0], y[0], z[0]} so long as it is at least 1000. Setting these initial values to 1/10 (infinite precision) is impractical, because the computation is so slow. Below is a plot of the Precision of the solution.

ListLinePlot[Transpose@Map[Precision, tab, {2}], ImageSize -> Large,
LabelStyle -> {15, Bold, Black}, PlotLegends -> Placed[{x, y, z}, {.3, .2}]]


Clearly, the Precision is more than adequate when the solution begins growing rapidly. This appears to be the true nature of the solution for r = 21. Perhaps, the oscillatory equilibrium of the recurrence relation becomes unstable with increasing r.

Work-Around

The code in the question is a fourth-order explicit Runge-Kutta solution of the the third order system of ODEs,

x'[t] == -x[t] + y[t];
y'[t] == -x[t] - y[t]/100 - y[t] z[t];
z'[t] == -z[t]/2 + x[t]^2 - r;


The explicit Runge-Kutta algorithm is unstable for some numerical parameters, and evidently that is the case here. If it is necessary for some reason to solve these ODEs using this particular algorithm, then decreasing h improves stability. For h = 1/200, the r = 21` solution is

Of course, there are many faster and more stable algorithms.