2
$\begingroup$

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}]

Could anyone please help me to solve this issue.

$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$
    – vicky
    Nov 8, 2020 at 9:25
  • $\begingroup$ Ok, i will wait to get the output. No can't reduce the number of iterations. $\endgroup$
    – vicky
    Nov 8, 2020 at 10:19
  • $\begingroup$ Did you ever obtain a solution by for r = 21 by waiting to get the output? $\endgroup$
    – bbgodfrey
    Nov 11, 2020 at 3:38

1 Answer 1

1
$\begingroup$

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] == .1`100000, y[0] == .1`100000, z[0] == .1`100000}, {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}]]

enter image description here

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}]]

enter image description here

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

enter image description here

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

$\endgroup$

Your Answer

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

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