Overflow problem implementing Euler's method

Questions
, ,
1

I tried to solve the following problem by the recursive method using the following code of the improved Euler’s method, but it returns me a message ‘‘Overflow occurred in computation’’

X'=Z+(Y-\alpha)X
Y'=1-\beta Y-X^2
Z'=-X-\gamma Z

with initial conditions (X(0),Y(0),Z(0)=(1,2,3).

Where

X: interest rate

Υ:investment demand

Z: price index

\alpha: savings,
\beta: cost per investment,
\gamma: the absolute value of the elasticity of demand

I tried to solve by the recursive method using the following code of the improved Euler’s method

Q[a_, b_, h_, N_] := (u[0] = a; v[0] = b; 
  Do[{u[n + 1] = 
     u[n] + h*
       F[u[n] + (h/2)*F[u[n], v[n]], 
        v[n] + (h/2)*
          G[u[n], v[
            n]]],                                                     \
          
    v[n + 1] = 
     v[n] + h*
       G[u[n] + (h/2)*F[u[n], v[n]], v[n] + (h/2)*G[u[n], v[n]]]}, {n,
     0, N}])
Open in Wolfram Cloud ↗
2

If you lower the number of iterations from 1,000 to just 10 you can see that u[i] gets very big very quickly:

In[179]:= Q[1, 3, 2, 1, 10]

In[180]:= Table[u[i], {i, 10}]

Out[180]= {6.79, 485.3302009, 8.035888983*10^9, 6.013776873*10^38, 1.886256893*10^154,1.825639218351085*10^616, 1.602032875976091*10^2464, 9.49942447295617*10^9855, 1.174358712296279*10^39423, 2.74292980559411*10^157691}

Not sure why this is, but maybe ‘h’ needs to be chosen smaller (usually h is a small number in these sorts of algorithms)

3

Athanasios,

as Arnoud Buzing pointed out your system is strongly diverging - but you could clearly see from this discussion that (at least with this initial conditions) this is not supposed to happen! Obviously the used step size h=1 is much too big, try e.g. h=0.1 .

As a side note: It is convention in WL that own variables, etc. should not start with a capital letter - for good reason! In your code you have:

Do[(*...some code ...*), {n, 0, N}]

But N is the name of a WL-function (giving numerical value)! I am surprised that this is working at all!

4

I have tried to do what you have advised me. Could you take a look please and tell me your opinion?
Thank you in advance.

Open in Wolfram Cloud ↗
5

Not sure what opinion you’re looking for, but that looks beautiful to me?

Problem solved?

6

Yes, that problem has been solved. I am so grateful for your help. Both of you

7

You can probably simplify things by vectorizing your code.

Instead of doing everything in triplicate and having lots of code repetition, you can perhaps just do this:

\[Alpha] = 0.9;
\[Beta] = 0.2;
\[Gamma] = 0.4;

f[{x_, y_, z_}] := {z + (y - \[Alpha])*x, 1 - \[Beta]*y - x^2, -x - \[Gamma]*z}

Q[a_, b_, c_, h_, NN_] := (
 u[0] = {a, b, c};
 Do[
   u[n + 1] = u[n] + h*f[u[n] + h/2*f[u[n]]],
   {n, 0, NN}]);

Q[1, 3, 2, 0.1, 200];
X = Interpolation[Table[{n, u[n]}, {n, 0, 2000}]];
ParametricPlot3D[X[t], {t, 0, 2000}]
8

Ok, thanks. But is not wrong the code above right?

9

Your code is fine, just trying to show you how to make it potentially simpler.

10

Thank you very much! I will keep it in mind!

11

Is it possible to contrast on the same axis system the two ParametricPlots from the first two cases on the [0,200]?

12

Yes, for example with this code:

MakePlot[{\[Alpha]_, \[Beta]_, \[Gamma]_}, {a_, b_, c_}, h_, max_Integer] := Module[{u, if},
  f[{x_, y_, z_}] := {z + (y - \[Alpha])*x, 1 - \[Beta]*y - x^2, -x - \[Gamma]*z};
  u[0] = {a, b, c};
  Do[u[n + 1] = u[n] + h*f[u[n] + h/2*f[u[n]]], {n, 0, max}];
  if = Interpolation[Table[{n, u[n]}, {n, 0, max/h}]];
  ParametricPlot3D[if[t], {t, 0, 2000}]
]

You can do this:

plot1 = MakePlot[{0.9, 0.2, 1.2}, {1, 3, 2}, 0.1, 200]
plot2 = MakePlot[{0.9, 0.2, 0.4}, {1, 3, 2}, 0.1, 200]

And then combine them:

Show[{plot1, plot2}]
13

That is very cool

14

Here is another idea, using AnimationVideo. Same MakePlot function, but with an explicit PlotRange to keep that constant for any parameter choice:

MakePlot[{\[Alpha]_, \[Beta]_, \[Gamma]_}, {a_, b_, c_}, h_, 
  max_Integer] := Module[{u, if},
  f[{x_, y_, z_}] := {z + (y - \[Alpha])*x, 
    1 - \[Beta]*y - x^2, -x - \[Gamma]*z};
  u[0] = {a, b, c};
  Do[u[n + 1] = u[n] + h*f[u[n] + h/2*f[u[n]]], {n, 0, max}];
  if = Interpolation[Table[{n, u[n]}, {n, 0, max/h}]];
  ParametricPlot3D[if[t], {t, 0, 2000}, 
   PlotRange -> {{-3, 3}, {-3, 3}, {-3, 3}}]
  ]

Then:

AnimationVideo[
 MakePlot[{0.9, 0.2, t}, {1, 1, 1}, 0.02, 1000],
 {t, 0.0, 2.0}]

Generates this video:

https://www.wolframcloud.com/obj/2d831a2f-4c28-4881-82bc-3bc95c88bbad

15

Amazing!!! Could I do the same for the

Plot[{X[t],Y[t],Z[t]},{t,0,200},PlotLegends->"Expressions"]
16

I can not do the comparison of two graphs like that

Open in Wolfram Cloud ↗