6
$\begingroup$

enter image description hereI have a system of non-linear differential equations: $\dot{x} = f(x,y,z),$ $\dot{y} = g(x,y,z)$ and $\dot{z} = h(x,y,z)$ with $0 \leq x,y,z \leq 1.$ The restpoints of the system are $(1,0,0),$ $(0,1,0),$ $(0,0,1)$ and $(0,0,0).$ I am trying to plot the basin of attractions of these restpoints and I was wondering if there is any plotting function in Mathematica that does this. Thanks a lot for the help in advance!

The code I am using is:

Clear[f, g, h, p, r, l, g, jac, u1, u2, u3, u4]

r = 1;  

G = {{3, 3, 3, 1}, {2.5, 2.5, 
   2.5, 0.5}, {2.5, 2.5, 2.5, 0.5}, {3.5, 3.5, 3.5, 1.5}};      
u1[x_, y_, z_] = 
  G[[1, 1]]*x +  G[[1, 2]]*y + G[[1, 3]]*z +  
   G[[1, 4]]*(1 - x - y - z) ;

u2[x_, y_, z_] = 
  G[[2, 1]]*x +  G[[2, 2]]*y + G[[2, 3]]*z +  
   G[[2, 4]]*(1 - x - y - z);

u3[x_, y_, z_] = 
  G[[3, 1]]*x +  G[[3, 2]]*y + G[[3, 3]]*z +  
   G[[3, 4]]*(1 - x - y - z) ;

u4[x_, y_, z_] = 
  G[[4, 1]]*x +  G[[4, 2]]*y + G[[4, 3]]*z +  
   G[[4, 4]]*(1 - x - y - z) ;

ualpha[x_, y_, z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]) + (z*
     u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

us[x_, y_, z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]);

ua[x_, y_, z_] = (z*u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

uc[x_, y_, z_] = (x*u1[x, y, z]) + (z*u3[x, y, z]);

ud[x_, y_, z_] = (y*u2[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

f[x_, y_, 
   z_] = ((1 - r)*x*u1[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - x;

g[x_, y_, 
   z_] = ((1 - r)*y*u2[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     ud[x, y, z]/((ualpha[x, y, z])^2)) - y;
h[x_, y_, 
   z_] = ((1 - r)*z*u3[x, y, z]/ualpha[x, y, z]) + (r*ua[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - z;

VectorPlot3D[{f[x, y, z], g[x, y, z], h[x, y, z]}, {x, 0, 1}, {y,0,1}, {z, 0, 1}]

The following code to produce trajectories from random initial points does not work. Any help fixing the code will be appreciated. Thanks!

nmax = 1000;
tmax = 100;
For[kn = 1, kn <= nmax, kn++, xinit = RandomReal[]; 
 yinit = RandomReal[{0, 1 - xinit}]; 
 zinit = RandomReal[{0, 1 - xinit - yinit}]; 
 solution = 
  NDSolve[{x'[t] == F1[x[t], y[t], z[t]], 
    y'[t] == F2[x[t], y[t], z[t]], z'[t] == F3[x[t], y[t], z[t]], 
    x[0] == xinit, y[0] == yinit, z[0] == zinit}, {x, y, z}, {t, 0, 
    tmax}]; Plot3D[{x[t], y[t], z[t]} /. solution, {t, 0, tmax}, {x, 
    0, 1}, PlotRange -> All, 
   BaseStyle -> Arrowheads[{0, .025, .025, 0}], 
   AxesLabel -> {"sc", "sd", "ac"}] /. Line -> Arrow]
$\endgroup$
5
  • $\begingroup$ Look for VectorPlot3D $\endgroup$ – Ulrich Neumann Jun 20 at 15:41
  • $\begingroup$ If you give your actual differential equations, or at least simplified versions thereof, people on this forum can experiment more broadly and come up with better proposals. $\endgroup$ – Roman Jun 20 at 16:01
  • $\begingroup$ I edited my post which has the code. The vector3dplot is producing a phase plot which is very hard to interpret. $\endgroup$ – egt123 Jun 20 at 16:21
  • $\begingroup$ Are you sure that all of the fixed points are stable? (I think you are missing one at {x -> 4/9, y -> 2/9, z -> 2/9}. I looks to me that only {0,0,0} and {1,0,0} are stable. $\endgroup$ – Craig Carter Jun 20 at 18:49
  • $\begingroup$ Yes, there is an interior rest point at (4/9,2/9,2/9). I checked the stability of the system using Jacobian matrix. The only stable rest points are (1,0,0) and (0,0,0). $\endgroup$ – egt123 Jun 20 at 19:05
5
$\begingroup$

We can have an attraction basin coarse image as follows

r = 1;
G = {{3, 3, 3, 1}, {2.5, 2.5, 2.5, 0.5}, {2.5, 2.5, 2.5, 0.5}, {3.5, 
3.5, 3.5, 1.5}};
u1[x_, y_, z_] := G[[1, 1]]*x + G[[1, 2]]*y + G[[1, 3]]*z + G[[1, 4]]*(1 - x - y - z);
u2[x_, y_, z_] := G[[2, 1]]*x + G[[2, 2]]*y + G[[2, 3]]*z + G[[2, 4]]*(1 - x - y - z);
u3[x_, y_, z_] := G[[3, 1]]*x + G[[3, 2]]*y + G[[3, 3]]*z + G[[3, 4]]*(1 - x - y - z);
u4[x_, y_, z_] := G[[4, 1]]*x + G[[4, 2]]*y + G[[4, 3]]*z + G[[4, 4]]*(1 - x - y - z);
ualpha[x_, y_, z_] := (x*u1[x, y, z]) + (y*u2[x, y, z]) + (z*
 u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);
us[x_, y_, z_] := (x*u1[x, y, z]) + (y*u2[x, y, z]);
ua[x_, y_, z_] := (z*u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);
uc[x_, y_, z_] := (x*u1[x, y, z]) + (z*u3[x, y, z]);
ud[x_, y_, z_] := (y*u2[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

f[x_, y_, z_] := ((1 - r)*x*u1[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*uc[x, y, z]/((ualpha[x, y, z])^2)) - x;
g[x_, y_, z_] := ((1 - r)*y*u2[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*ud[x, y, z]/((ualpha[x, y, z])^2)) - y;
h[x_, y_, z_] := ((1 - r)*z*u3[x, y, z]/ualpha[x, y, z]) + (r*ua[x, y, z]*uc[x, y, z]/((ualpha[x, y, z])^2)) - z;

ODEs = {x'[t] == f[x[t], y[t], z[t]], y'[t] == g[x[t], y[t], z[t]], z'[t] == h[x[t], y[t], z[t]]};

tmax = 30;
tol = 0.01;
p1 = {1, 0, 0};
P0 = {};
P1 = {};
nmax = 20000;
For[k = 1, k <= nmax, k++,
 rand = RandomReal[{0, 1}, 3];
 cinits = Thread[{x[0], y[0], z[0]} == rand];
 solODE = NDSolve[Join[ODEs, cinits], {x, y, z}, {t, 0, tmax}][[1]];
 prox = Evaluate[{x[t], y[t], z[t]} /. solODE] /. {t -> tmax};
 d0 = Norm[prox];
 d1 = Norm[p1 - prox];
 If[d0 > tol && d1 > tol, Continue[], If[d0 < d1, AppendTo[P0, rand], AppendTo[P1, rand]]]
 ]

Show[ConvexHullMesh[P0, MeshCellStyle -> {{2, All} -> Opacity[0.5, Orange]}], 
     ConvexHullMesh[P1, MeshCellStyle -> {{2, All} -> Opacity[0.5, Blue]}]]

enter image description here

$\endgroup$
5
$\begingroup$

I can recommend ParametricNDSolve[] to compute trajectories and ParametricPlot3D[] for visualization, for example,

Clear[f, g, h, p, r, l, g, jac, u1, u2, u3, u4]

r = 1;

G = {{3, 3, 3, 1}, {2.5, 2.5, 2.5, 0.5}, {2.5, 2.5, 2.5, 0.5}, {3.5, 
    3.5, 3.5, 1.5}};
u1[x_, y_, z_] = 
  G[[1, 1]]*x + G[[1, 2]]*y + G[[1, 3]]*z + G[[1, 4]]*(1 - x - y - z);

u2[x_, y_, z_] = 
  G[[2, 1]]*x + G[[2, 2]]*y + G[[2, 3]]*z + G[[2, 4]]*(1 - x - y - z);

u3[x_, y_, z_] = 
  G[[3, 1]]*x + G[[3, 2]]*y + G[[3, 3]]*z + G[[3, 4]]*(1 - x - y - z);

u4[x_, y_, z_] = 
  G[[4, 1]]*x + G[[4, 2]]*y + G[[4, 3]]*z + G[[4, 4]]*(1 - x - y - z);

ualpha[x_, y_, 
   z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]) + (z*
     u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

us[x_, y_, z_] = (x*u1[x, y, z]) + (y*u2[x, y, z]);

ua[x_, y_, z_] = (z*u3[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

uc[x_, y_, z_] = (x*u1[x, y, z]) + (z*u3[x, y, z]);

ud[x_, y_, z_] = (y*u2[x, y, z]) + ((1 - x - y - z)*u4[x, y, z]);

f[x_, y_, 
   z_] = ((1 - r)*x*u1[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - x;

g[x_, y_, 
   z_] = ((1 - r)*y*u2[x, y, z]/ualpha[x, y, z]) + (r*us[x, y, z]*
     ud[x, y, z]/((ualpha[x, y, z])^2)) - y;
h[x_, y_, 
   z_] = ((1 - r)*z*u3[x, y, z]/ualpha[x, y, z]) + (r*ua[x, y, z]*
     uc[x, y, z]/((ualpha[x, y, z])^2)) - z;

Code to compute and visualize trajectories

s = ParametricNDSolve[{x'[t] == f[x[t], y[t], z[t]], 
   y'[t] == g[x[t], y[t], z[t]], z'[t] == h[x[t], y[t], z[t]], 
   x[0] == x0, y[0] == 0., z[0] == 0.}, {x, y, z}, {t, 0, 100}, {x0}];

ParametricPlot3D[
 Table[{x[x0][t], y[x0][t], z[x0][t]} /. s, {x0, 0.05, 1, .05}], {t, 
  0, 100}, PlotRange -> All, ColorFunction -> Hue, 
 PlotTheme -> "Marketing", AxesLabel -> {"x", "y", "z"}]

Figure 1

In the picture above there are two points of attraction - $(0,0,0)$ and $(1,0,0)$, With initial condition y[0]==0.5, z[0]==0 we also have these points

s1 = ParametricNDSolve[{x'[t] == f[x[t], y[t], z[t]], 
    y'[t] == g[x[t], y[t], z[t]], z'[t] == h[x[t], y[t], z[t]], 
    x[0] == x0, y[0] == 0.5, z[0] == 0.}, {x, y, z}, {t, 0, 
    100}, {x0}];

ParametricPlot3D[
 Table[{x[x0][t], y[x0][t], z[x0][t]} /. s1, {x0, 0.05, 1, .05}], {t, 
  0, 100}, PlotRange -> All, ColorFunction -> Hue, 
 PlotTheme -> "Marketing", AxesLabel -> {"x", "y", "z"}]

Figure 2

In a case of 100 trajectories with initial condition $x(0)+y(0)+z(0)=1$ we can use this code

s2 = ParametricNDSolve[{x'[t] == f[x[t], y[t], z[t]], 
   y'[t] == g[x[t], y[t], z[t]], z'[t] == h[x[t], y[t], z[t]], 
   x[0] == (Cos[p1] Cos[p2])^2, y[0] == (Cos[p1] Sin[p2])^2, 
   z[0] == Sin[p1]^2}, {x, y, z}, {t, 0, 100}, {p1, p2}]


P1 = RandomReal[{0, Pi}, 10]; P2 = 
 RandomReal[{0, Pi}, 10]; ParametricPlot3D[
 Table[{x[p1, p2][t], y[p1, p2][t], z[p1, p2][t]} /. s2, {p1, 
    P1}, {p2, P2}] // Flatten, {t, 0, 100}, PlotRange -> All, 
 ColorFunction -> Hue, PlotTheme -> "Marketing", 
 AxesLabel -> {"x", "y", "z"}]

Figure 3

$\endgroup$
7
  • $\begingroup$ That was very helpful. Thanks a lot ! Is it also possible to divide the xyz space into two regions, where one region is the basin of attraction of the stable restpoint (1,0,0) and the second region is the basin of attraction of the stable restpoint (0,0,0)? $\endgroup$ – egt123 Jun 20 at 18:25
  • $\begingroup$ @egt123 Practically yes, we can divide region as Cesario did. $\endgroup$ – Alex Trounev Jun 20 at 20:31
  • $\begingroup$ Is it possible to draw the trajectories starting from (let's say) 1000 random points from the space (x,y,z) with x+y+z<=1 ? I tried something but the code does not work. I edited my original post and put the code which I have been trying to produce trajectories from random initial points. $\endgroup$ – egt123 Jun 22 at 15:39
  • $\begingroup$ @egt123 Do you need $0\le x,y,z \le 1$ or $x+y+z=1$? $\endgroup$ – Alex Trounev Jun 22 at 17:42
  • $\begingroup$ I need $x,y,z \geq 0$ and $x+y+z \leq 1.$ $\endgroup$ – egt123 Jun 22 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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