# 3d system: stream plot and function

I have a 3D ODE-system $$x'=f_1(x,y,z),\quad y'=f_2(x,y,z),\quad z'=f_3(x,y,z)$$

living on the sphere and an equilibrium $$(x_0,y_0,z_0)$$. I would like to plot a function $$V(x,y,z)$$ to see whether t could be a Lyapunov function for this equilibrium. To this end, it would be great to see the stream plot and the function V.

To be more concrete: I have $$x'=x(-x+f(x,y,z)),\quad y'=y(x-y+f(x,y,z)),\quad z'=z(y-z+f(x,y,z))$$

with $$f(x,y,z)=x^3-xy^2+y^3-yz^2+z^3$$ and the equilibrium $$(a,2a,3a)$$ with $$a=\sqrt{1/14}$$.

Here $$(x,y,z)\in S^2$$ (one-sphere).

Would like to see the stream plot and to plot $$V(x,y,z)=(x-a)^2+(y-2a)^2+(z-3a)^2$$.

In particular, I would like to see graphically, if $$V$$ can be a global Lyapunov-function for $$(a,2a,3a)$$ on $$M:=\{(x,y,z)\in S^2: x,y,z>0\}$$. This means that $$V>0$$ on $$M\setminus\{(a,2a,3a)\}$$ and $$V'<0$$ on $$M\setminus\{(a,2a,3a)\}$$. Maybe one can see on a plot that this is not the case? Then any further analysis would be superfluous...

I tried this for the stream plot (only plotting x and y).

z = (1-x^2-y^2)^(1/2)
s= x^3 - x*y^2+y^3-y*z^2+z^3
xdot = x * (-x+s)
ydot = y* (x-y+s)

StreamPlot[{xdot,ydot},{x,0,1},{y,0,1}]


• Hello, welcome to Mathematica.SE. Then, please show us the code text rather than the screenshot of it, so we can easily test. BTW, sadly there's no StreamPlot3D in Mathematica at the moment. Related: 1. mathematica.stackexchange.com/q/137809/1871 2. mathematica.stackexchange.com/q/123137/1871 There may be more… – xzczd Mar 29 '20 at 13:30
• I added the code and removed the screenshot. - But maybe it is possible to plot everything on the sphere? How can I plot the function $V(x,y,z)$ to see where on the sphere it is positive and decreasing? – Scuderi Mar 29 '20 at 13:36
• You mean this?: mathematica.stackexchange.com/q/65401/1871 – xzczd Mar 29 '20 at 13:38
• Yes, it would be great to have something like that. I also saw this post bu it was too difficult to me to understand how to use this for my system. – Scuderi Mar 29 '20 at 13:40
• However, this then does not capture the $z$-dynamics but is only a projection of the x-y-dynamics onto the sphere, isn't it? – Scuderi Mar 29 '20 at 13:47

OK, after reading Alex's answer and your post in math.SE, I think I somewhat (at least partly) understand the question. First, as to your doubt in the comment:

However, this then does not capture the $$z$$ dynamics but is only a projection of the $$x$$-$$y$$-dynamics onto the sphere, isn't it?

The answer is no, because the independent variables in that post are $$\theta$$ and $$\phi$$, the angle components in spherical coordinate. And we can draw the same plot for your vector with proper coordinate transform:

f = x^3 - x y^2 + y^3 - y z^2 + z^3;
xdot = x (-x + f);
ydot = y (x - y + f);
zdot = z (y - z + f);
vector = {xdot, ydot, zdot};

a = 1/Sqrt[14];

point = {a, 2 a, 3 a};

transformedvector =
TransformedField["Cartesian" -> "Spherical",
vector, {x, y, z} -> {r, theta, phi}] /. r -> 1 // Simplify
(*
{0,
1/4 Cos[theta] Sin[theta] (4 Cos[theta] + (-2 Cos[phi] - 2 Cos[3 phi] - 7 Sin[phi] +
Sin[3 phi]) Sin[theta]),
Cos[phi] (2 Cos[phi] - Sin[phi]) Sin[phi] Sin[theta]^2}
*)


As one can see, after substituting $$r=1$$ into the transformed vector, $$r$$ component of the vector becomes $$0$$, which indicates all the vectors at $$r=1$$ are on the unit ball.

Then we plot it in 2D:

plot = StreamPlot[
transformedvector // Rest // Evaluate, {theta, 0, Pi}, {phi, 0, 2 Pi}]~Show~
Graphics@{Orange, PointSize@Large,
Point@Rest@CoordinateTransform["Cartesian" -> "Spherical", point]}


Well, in my opinion this illustration is already good enough, but if you insist on visualising on the ball:

func = {theta, phi} \[Function]
Evaluate@CoordinateTransform[
"Spherical" -> "Cartesian", {1, theta, phi}]

plot3D = Graphics3D[(plot[[1]] /. (head : Arrow | Point)[z_] :>
head[z /. {x_?NumericQ, y_} :> func @@ {x, y}])]


plot3D~Show~Graphics3D@Ball[{0, 0, 0}, 0.98]


Visualization of the solution in the form of trajectories on the sphere

eq = {x[t] (-x[t] + x[t]^3 - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 +
z[t]^3),
y[t] (x[t] + x[t]^3 - y[t] - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 +
z[t]^3),
z[t] (x[t]^3 + y[t] - x[t] y[t]^2 + y[t]^3 - z[t] - y[t] z[t]^2 +
z[t]^3)};

sol = ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]},
x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c],
z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, 20}, {c, b}]

a = 1/Sqrt[14.]; Show[
Graphics3D[{{Green, Ball[]}, {Orange, PointSize[.05],
Point[{a, 2 a, 3 a}]}}, Boxed -> False],
ParametricPlot3D[
Evaluate[Table[sol[Pi/12, b], {b, 0, 2 Pi, .1}]], {t, 0, 10},
PlotRange -> All],
ParametricPlot3D[
Evaluate[Table[sol[Pi/3, b], {b, 0, 2 Pi, .1}]], {t, 0, 10},
PlotRange -> All]]


Let's check if the point {a,2 a,3 a} is Lyapunov stable. We linearize the equation in a neighborhood of this point

eq1 = eq /. {x[t] -> a + e x1[t], y[t] -> 2 a + e y1[t],
z[t] -> 3 a + e z1[t]};
s1 = Series[eq1, {x1[t], 0, 1}, {y1[t], 0, 1}, {z1[t], 0, 1}] ;
eqL = s1 // Normal;
eql = Series[eqL, {e, 0, 1}] // Normal; eqlp = Chop[eql /. e -> 1]

(*Out[]= {-0.286351 x1[t] - 0.0190901 y1[t] + 0.286351 z1[t],
0.496342 x1[t] - 0.572703 y1[t] + 0.572703 z1[t], -0.0572703 x1[t] +
0.744513 y1[t] + 0.0572703 z1[t]}*)


Matrix of linear system X'[t] =A.X

A = CoefficientArrays[eqlp, {x1[t], y1[t], z1[t]}] // Normal // Last

(*Out[]= {{-0.286351, -0.0190901, 0.286351}, {0.496342, -0.572703,
0.572703}, {-0.0572703, 0.744513, 0.0572703}}*)


Finally check

LyapunovSolve[
Transpose[
A], -{{1, 0, 0}, {0, 2, 0}, {0, 0,
3}}] // PositiveDefiniteMatrixQ

(*Out[]= False*)


Therefore, the system is unstable. Eigenvalues

Eigenvalues[A]

(*Out[]= {-0.801784, 0.534522, -0.534522}*)


Solution close to the point $$(a,2a,3a)$$. Code for v.12.1:

sol = ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]},
x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c],
z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, 100}, {c, b}];
Show[Graphics3D[{Green, Ball[]}, Boxed -> False],
ParametricPlot3D[
Evaluate[Table[sol[Pi/12, b], {b, .01, Pi/2, .02}]], {t, 0, 30},
PlotRange -> All],
ParametricPlot3D[
Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .02}]], {t, 0, 30},
PlotRange -> All]]
{ParametricPlot3D[
Evaluate[Table[sol[Pi/12, b], {b, .01, Pi/2, .01}]], {t, 0, 30},
PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}],
ParametricPlot3D[
Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .01}]], {t, 0, 30},
PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}]}


We see that the trajectories leave the sphere near the point $$(a, 2a, 3a)$$.Code for v.12.0:

eq = {x[t] (-x[t] + x[t]^3 - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 +
z[t]^3),
y[t] (x[t] + x[t]^3 - y[t] - x[t] y[t]^2 + y[t]^3 - y[t] z[t]^2 +
z[t]^3),
z[t] (x[t]^3 + y[t] - x[t] y[t]^2 + y[t]^3 - z[t] - y[t] z[t]^2 +
z[t]^3)}; tm = 23; sol =
ParametricNDSolveValue[{eq == {x'[t], y'[t], z'[t]},
x[0] == Cos[b] Sin[c], y[0] == Sin[b] Sin[c],
z[0] == Cos[c]}, {x[t], y[t], z[t]}, {t, 0, tm}, {c, b}];
a = 1/Sqrt[14];

Show[Graphics3D[{Green, Opacity[.4], Sphere[]},
PlotRange -> {{0, 1}, {0, 1}, {1/4, 1}}, Boxed -> False],
ParametricPlot3D[
Evaluate[Table[sol[Pi/6, b], {b, .01, Pi/2, .02}]], {t, 0, tm},
PlotRange -> All],
ParametricPlot3D[
Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .02}]], {t, 0, tm},
PlotRange -> All]]

Show[{ParametricPlot3D[
Evaluate[Table[sol[Pi/6, b], {b, .01, Pi/2, .01}]], {t, 0, tm},
PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}],
ParametricPlot3D[
Evaluate[Table[sol[Pi/3, b], {b, 0.01, Pi/2, .01}]], {t, 0, tm},
PlotRange -> All, Boxed -> False, AxesLabel -> {x, y, z}]}]


• Wow! The plot suggests that (a,2a,3a) is globally stable on M, doesn‘t it? So there Must be some Lyapunov-function on M. Is there a way to plot functions V(x,y,z) in this plot or another to see if they are possible candidates? – Scuderi Mar 29 '20 at 21:07
• @Scuderi See update to my answer. – Alex Trounev Mar 29 '20 at 21:49
• I agree that in $3d$ the point $(a,2a,3a)$ is a saddle and thus unstable. But isn't this one dimension too much? Because we are only living on the sphere. In some sense the point must be stable as the plot suggests. – Scuderi Mar 30 '20 at 9:06
• In a numerical solution, we cannot get too close to this point. Therefore, there is an instability that is visible with good resolution - see update. – Alex Trounev Mar 30 '20 at 15:17