Simulation Lorenz83 Attractor

Question

I tried to simulate the Lorenz83 Attractor that is defined by the following system of nonlinear ordinary differential equations: \begin{eqnarray*} \frac{dx}{dt}&=&-a x - y^2 - z^2 + a f\\ \frac{dy}{dt}&=& -y + x y - b x z + g\\ \frac{dz}{dt}&=& -z + b x y + x z \end{eqnarray*} Given the following parameters and initial conditions: $a=0.95,b=7.91,f=4.83,g=4.66$ and $ x(0)=-0.2,y(0)=-2.82,z(0)=2.71$

(* Parameters for the system *)
parameters = {a -> 0.95, b -> 7.91, f -> 4.83, g -> 4.66};

(* System of differential equations *)
equations = {
   x'[t] == -a x[t] - y[t]^2 - z[t]^2 + a f,
   y'[t] == -y[t] + x[t] y[t] - b x[t] z[t] + g,
   z'[t] == -z[t] + b x[t] y[t] + x[t] z[t]
} /. parameters;

(* Initial conditions *)
initialConditions = {x[0] ==-0.2, y[0] == -2.82, z[0] == 2.71};

(* Numerical solution *)
solution = NDSolve[{equations, initialConditions}, {x, y, z}, {t, 0, 20},Method->"StiffnessSwitching"];

(* 3D plot of the attractor with color gradient *)
trajectoryPlot = ParametricPlot3D[
   Evaluate[{x[t], y[t], z[t]} /. solution], {t, 0, 20},
   PlotRange -> All, AxesLabel -> {"x", "y", "z"},
   PlotStyle -> {Thick,Dashed}, ColorFunction -> "DarkRainbow", ColorFunctionScaling -> True,Boxed -> True
,PlotPoints -> 20];
trajectoryPlot

but it doesn't seem right because the picture given by the source Strange Attractors I'm studying is different

Answer 1

When using PlotPoints -> 1000, the visualization improves. However, if the time exceeds 20, no chaotic behavior is observed with these parameter values ​​and initial condition, and the system reaches a limit cycle.

Plot[x[t] /. solution, {t, 0, 80}, PlotRange -> All, 
PlotStyle -> {Thickness[0.001], Blue}, PlotPoints -> 200, 
PlotHighlighting -> None]

If you change the value of b, there is chaotic behavior:

(*Parameters for the system*)
parameters = {a -> 0.95, b -> 5.71, f -> 4.83, g -> 4.66};

(*System of differential equations*)
equations = {x'[t] == -a  x[t] - y[t]^2 - z[t]^2 + a  f, 
y'[t] == -y[t] + x[t]  y[t] - b  x[t]  z[t] + g, 
z'[t] == -z[t] + b  x[t]  y[t] + x[t]  z[t]} /. parameters;

(*Initial conditions*)
initialConditions = {x[0] == -0.2, y[0] == -2.82, z[0] == 2.71};

(*Numerical solution*)
solution = 
NDSolve[{equations, initialConditions}, {x, y, z}, {t, 0, 1000}];

Plot[x[t] /. solution, {t, 900, 1000}, PlotRange -> All, 
PlotStyle -> {Thickness[0.001], Blue}, PlotPoints -> 200, 
PlotHighlighting -> None]

Following a comment from @Chris K, here is the link to the model reference:

Irregularity: a fundamental property of the atmosphere by Edward N. Lorenz (1984)