# Plotting a 3d surface along the vector field flow

Consider the following example of the Lorentz system

f[x_, y_, z_] := σ (y - x);
g[x_, y_, z_] := x (ρ - z) - y;
h[x_, y_, z_] := x y - β z;
set = {σ -> 10, ρ -> 28, β -> 8/3};
tf = 50;
sol = NDSolve[{
x'[t] == f[x[t], y[t], z[t]] /. set,
y'[t] == g[x[t], y[t], z[t]] /. set,
z'[t] == h[x[t], y[t], z[t]] /. set,
x == 1,
y == 1,
z == 1}
, {x, y, z}
, {t, 0, tf}
, MaxSteps -> ∞
];
ParametricPlot3D[
Evaluate[ {x[t], y[t], z[t]} /. sol ]
, {t, 0, tf}
, BoxRatios -> {1, 1, 1}
, PlotRange -> All
] How can the surface along the direction of the vector field flow (around the two attractors) be plotted? (as in the very inaccurate example below) • There are things to do after your question is answered. It's a good idea to stay vigilant for some time, because better approaches may come later improving over previous replies. Experienced users may point alternatives, caveats or limitations. New users should test answers before voting and wait 24 hours before accepting the best one. But at some point more answers become unlikely, and it's time to accept and vote. Participation is essential for the site, please do your part. – rhermans Jul 6 '18 at 23:02
• Probably what you really want is this: find the complex-conjugate eigenvectors of the Jacobian at the two fixed points, form their sum and difference to get two real vectors that span the plane through each fixed point. Then draw the planes. This is the well-defined way of doing things without guessing based on the numerical points you're plotting. – Jens Jul 7 '18 at 0:13

Here's a somewhat crude idea, using your own definitions:

res = Table[Evaluate[{x[t], y[t], z[t]} /. First@sol], {t, 0, 50, .01}];
split = GatherBy[res, #[] >= 0 &];

NonlinearModelFit[#, a x + b y + c, {a, b, c}, {x, y}] & /@ split;

Show[
ListPointPlot3D[split, BoxRatios -> {1, 1, 1}],
Plot3D[
Evaluate@Through[%[x, y]],
{x, -20, 20}, {y, -20, 20},
PlotStyle -> Opacity[0.5]
]
] In essence, I split the points in the attractor in two groups, the ones with positive vs. negative $x$-coordinate, then find a best-fit plane through each group, and display the results.

Similar approach to @MarcoB but using FindClusters and Fit.

Using sol from OP.

res = Table[Evaluate[{x[t], y[t], z[t]} /. First@sol], {t, 0, 50, .01}];
cluster = FindClusters[res, 2, Method -> "KMeans"];


Then

Show[
ListPointPlot3D[cluster],
Plot3D[
Evaluate[Fit[#, {1, x, y}, {x, y}]],
Evaluate[
Sequence @@
MapThread[Prepend, {Most[MinMax /@ Transpose@#], {x, y}}]],
PlotStyle -> None,
MeshStyle -> Opacity[.25]
] & /@ cluster,
Boxed -> False,
AxesEdge -> {{-1, -1}, {1, -1}, {1, -1}}
] Hope this helps.

OK, here is the approach I mentioned in the comments:

Clear[x, y, z];
f[x_, y_, z_] := σ (y - x);
g[x_, y_, z_] := x (ρ - z) - y;
h[x_, y_, z_] := x y - β z;
set = {σ -> 10, ρ -> 28, β -> 8/3};
tf = 50;
sol = NDSolve[{x'[t] == f[x[t], y[t], z[t]] /. set,
y'[t] == g[x[t], y[t], z[t]] /. set,
z'[t] == h[x[t], y[t], z[t]] /. set, x == 1, y == 1,
z == 1}, {x, y, z}, {t, 0, tf}, MaxSteps -> Infinity];

fixedPoints = {x, y, z} /.
Rest[Solve[{f[x, y, z] == 0, g[x, y, z] == 0, h[x, y, z] == 0}, {x,
y, z}]]

(*
==> {{-Sqrt[β] Sqrt[-1 + ρ], -Sqrt[β]
Sqrt[-1 + ρ], -1 + ρ}, {Sqrt[β]
Sqrt[-1 + ρ], Sqrt[β] Sqrt[-1 + ρ], -1 + ρ}}
*)

jac[{x_, y_, z_}] =
D[{f[x, y, z], g[x, y, z], h[x, y, z]}, {{x, y, z}}];

planeVectors =
Chop@Table[{vec[] + vec[], I (vec[] - vec[])}, {vec,
N[(Rest[Eigenvectors[jac[#]]] & /@ fixedPoints) /. set]}];

Show[
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol], {t, 0, tf},
PlotStyle -> Tube[.1], BoxRatios -> {1, 1, 1}, PlotRange -> All],
ParametricPlot3D[
MapThread[(#1 + {u, v}.#2) /. set &, {fixedPoints,
planeVectors}], {u, -10, 10}, {v, -10, 10}, Mesh -> None,
PlotStyle -> Opacity[.5], BoundaryStyle -> Directive[Thickness[.01], Black]]] The fixed points are obtained by setting the functions f, g, h to zero. I then calculate the Jacobian jac. All this can be done symbolically.

To draw the planes, I then insert the numerical parameters into the eigenvectors of jac, and retain only the ones corresponding to pairs of complex conjugate eigenvalues. This is done using Rest because of the way the results are sorted.

The complex-conjugate pairs of eigenvectors that are left can then be combined to get real vectors that span the desired plane. The two linear combinations I choose are simply their sum and difference (where the difference must be multiplied by i to make it real).

Once you have a pair of real vectors attached to each fixed point, you can sweep out the desired surface by multiplying one by u and the other by v, and then varying u, v as independent variables. This is done in ParametricPlot3D, which is then combined with the plot of the trajectory that was already posted in the question.

To make the trajectory stand out more, I turned it into a Tube.