# Plotting Basins of Attraction of 3-dimensional system

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

• Look for VectorPlot3D Jun 20, 2021 at 15:41
• 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. Jun 20, 2021 at 16:01
• I edited my post which has the code. The vector3dplot is producing a phase plot which is very hard to interpret. Jun 20, 2021 at 16:21
• 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. Jun 20, 2021 at 18:49
• 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). Jun 20, 2021 at 19:05

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


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


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


• 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)? Jun 20, 2021 at 18:25
• @egt123 Practically yes, we can divide region as Cesario did. Jun 20, 2021 at 20:31
• 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. Jun 22, 2021 at 15:39
• @egt123 Do you need $0\le x,y,z \le 1$ or $x+y+z=1$? Jun 22, 2021 at 17:42
• I need $x,y,z \geq 0$ and $x+y+z \leq 1.$ Jun 22, 2021 at 17:50

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