Lorenz Equations

σ=40;
b=8/3;
f[x_,y_,z_]:=σ*(y-x);
g[x_,y_,z_]:=r*x-y-x*z;
h[x_,y_,z_]:=x*y-b*z;

Fixed Point Analysis

Simplify[Solve[{f[x,y,z]==0,g[x,y,z]==0,h[x,y,z]==0},{x,y,z}]]

MatrixForm[D[{f[x,y,z],g[x,y,z],h[x,y,z]},{{x,y,z}}]]

J[x_,y_,z_]:=(({{-10, 10, 0}, {r-z, -1, -x}, {y, x, -(8/3)}}))

Plot[Re[Eigenvalues[J[0,0,0]]],{r,0,3}]

Plot[Re[Eigenvalues[J[-2 (Sqrt[(2/3)]) (Sqrt[-1+r]),-2 (Sqrt[(2/3)]) (Sqrt[-1+r]),-1+r]]],{r,0,60},PlotRange->{-15,5}]

Plot[Im[Eigenvalues[J[-2 (Sqrt[(2/3)]) (Sqrt[-1+r]),-2 (Sqrt[(2/3)]) (Sqrt[-1+r]),-1+r]]],{r,0,60}]

There is a saddle node bifurcation and a Hopf bifurcation in this system. What about the dynamics?

Numerical Experiments

r is the bifurcation parameter. T is the length of time for numerical integration. If the plot is jagged reduce T to reduce interpolation over long times. R is the domain size for plotting.

r=30;
T=100;
R=20;

x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.1;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-3*R,3*R},{-3*R,3*R},{0,3*R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot1z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot1x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.2;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot2z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot2x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.3;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot3=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot3z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot3x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.4;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot4=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot4z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot4x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.5;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot5=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot5z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot5x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.6;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot6=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot6z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot6x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.7;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot7=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot7z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot7x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
x0=RandomReal[{-R/2,R/2}];
y0=RandomReal[{-R/2,R/2}];
z0=RandomReal[{0,R/2}];
p=.8;
s1=NDSolve[{x'[t]==σ*(y[t]-x[t]),y'[t]==r*x[t]-y[t]-x[t]*z[t],z'[t]==x[t]*y[t]-b*z[t],x[0]==x0,y[0]==y0,z[0]==z0},{x,y,z},{t,0,T}];
plot8=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s1],{t,0,T},PlotRange->{{-R,R},{-R,R},{-R,R}},PlotStyle->Hue[p],ViewPoint->{6,3,2}];
plot8z=Plot[Evaluate[z[t]/.s1],{t,0,T},PlotRange->{0,R},PlotStyle->Hue[p],ImageSize->Large];
plot8x=Plot[Evaluate[x[t]/.s1],{t,0,T},PlotRange->{-R,R},PlotStyle->Hue[p],ImageSize->Large];
Show[plot1,plot2,plot3,plot4,plot5,plot6,plot7,plot8]
Show[plot1z,plot2z,plot3z,plot4z,plot5z,plot6z,plot7z,plot8z]
Show[plot1x,plot2x,plot3x,plot4x,plot5x,plot6x,plot7x,plot8x]
Show[plot1,plot2]
Show[plot1z,plot2z]
Show[plot1x,plot2x]