Chaotic systems

The Lorenz mill, two initial conditions

Clear[a,b,ω,da0,db0,dω0]
tmax=50;
Manipulate[
a0=0.1;b0=0.15;ω0=-0.2;
K=0.129;q=0.11;ν=1.2;g=9.8;r=2.3;Ii=0.3;
integration1=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0,b[0]==b0,ω[0]== ω0},{a,b,ω},{t,0,tmax}];
integration2=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0+da0,b[0]==b0+db0,ω[0]== ω0+dω0},{a,b,ω},{t,0,tmax}];
solω1=Plot[Evaluate[ω[t]/.integration1],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,ω[t]}];
solω2=Plot[Evaluate[ω[t]/.integration2],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,ω[t]},PlotStyle->{Red}];
Show[solω1,solω2,PlotLabel->Style["Angular velocity of the wheel",20],FrameLabel-> {t,ω[t]},LabelStyle->20,Frame->True,ImageSize->Large],
Control[{{da0,0.00001,"Deviation in a direction, (Subscript[Δa, 0])"},-0.00010,0.00010,ImageSize->Large,Appearance->"Labeled"}],
Control[{{db0,0,"Deviation in b direction, (Subscript[Δb, 0])"},-0.00010,0.00010,ImageSize->Large,Appearance->"Labeled"}],
Control[{{dω0,0,"Deviation in ω direction, (Subscript[Δω, 0])"},-0.00010,0.00010,ImageSize->Large,Appearance->"Labeled"}],
SaveDefinitions->True,TrackedSymbols:>{da0,db0,dω0}]

The Lorenz mill (with chaotic parameters)

Clear[a,b,ω]
tmax=40;

a0=0.1;
b0=0.15;
ω0=-0.2;

K=0.129;
q=0.11;(*outflow distribution, fourier amplitude, first mode*)
ν=1.2;
g=9.8;
r=2.3;
Ii=0.3;

diffa=ω b-K a;diffb=-ω a-K b+q;diffω=((-ν)/Ii) ω+((π g r)/Ii) a;
fixedpoints=Quiet[Solve[{diffa==0,diffb==0,diffω==0},{a,b,ω}]];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={a,b,ω}/.fixedpoints
integration=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0,b[0]==b0,ω[0]== ω0},{a,b,ω},{t,0,tmax}];
Plot[Evaluate[ω[t]/.integration],{t,0,tmax},PlotRange-> All,PlotLabel->Style["Angular velocity of the wheel",20],FrameLabel-> {t,ω[t]},LabelStyle->20,Frame->True]

Spectra of the above signals

Clear[a,b,ω]
tmax=250;

a0=0.1;
b0=0.15;
ω0=-0.2;
da0=0.00001;

K=0.129;
q=0.11;
ν=1.2;
g=9.8;
r=2.3;
Ii=0.3;

integration1=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0,b[0]==b0,ω[0]== ω0},{a,b,ω},{t,0,tmax}];
integration2=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0+da0,b[0]==b0,ω[0]== ω0},{a,b,ω},{t,0,tmax}];
solω1=Plot[Evaluate[ω[t]/.integration1],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,ω[t]},ImageSize->330,LabelStyle->20];
solω2=Plot[Evaluate[ω[t]/.integration2],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,ω[t]},PlotStyle->{Red},ImageSize->330,LabelStyle->20];
sr=10;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[Evaluate[ω[t]/.integration1[[1]]],{t,0,tmax,1/sr}];(*sample values*)
signal2=Table[Evaluate[ω[t]/.integration2[[1]]],{t,0,tmax,1/sr}];(*sample values*)
ft=Fourier[signal,FourierParameters->{-1,-1}];(*num fourier transform*)
ft2=Fourier[signal2,FourierParameters->{-1,-1}];(*num fourier transform*)
ff=Table[2 π (n-1) sr/Length[signal],{n,Length[signal]}];(*frequency values x-axis*)
specter1=ListLinePlot[Transpose[{ff,Log[Abs[ft]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},PlotStyle->Thickness[0.002],ImageSize->380,LabelStyle->20];
specter2=ListLinePlot[Transpose[{ff,Log[Abs[ft2]^2]}][[1;;Round[0.5 Length[signal]]]],PlotRange->All,Frame->True,FrameLabel->{"Frequency ω","Spectral Power Level"},PlotStyle->{Thickness[0.001],Red},ImageSize->380,LabelStyle->20];
Row[{Show[solω1,solω2],Show[specter1,specter2]}]
Row[{solω1,specter1}]
Row[{solω2,specter2}]

Full 3-D solution, the Lorenz mill

Clear[a,b,ω]
tmax=400;

a0=0.1;
b0=0.15;
ω0=-0.2;

K=0.129;
q=0.11;
ν=1.2;
g=9.8;
r=2.3;
Ii=0.3;

integration=NDSolve[{
a'[t]==ω[t] b[t]-K a[t],
b'[t]==-ω[t] a[t]-K b[t]+q,
ω'[t]==((-ν)/Ii) ω[t]+((π g r)/Ii) a[t],
a[0]==a0,b[0]==b0,ω[0]== ω0},{a,b,ω},{t,0,tmax}];
diffa=ω b-K a;diffb=-ω a-K b+q;diffω=((-ν)/Ii) ω+((π g r)/Ii) a;
fixedpoints=Quiet[Solve[{diffa==0,diffb==0,diffω==0},{a,b,ω}]];
Print[Style["Coordinates of fixed points",18,Blue]]
coordinates={a,b,ω}/.fixedpoints
Gfp=Graphics3D[{PointSize[0.02],Red,Point[coordinates]}];
sola=Plot[Evaluate[a[t]/.integration],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,a[t]},ImageSize->240,LabelStyle->20];
solb=Plot[Evaluate[b[t]/.integration],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,b[t]},ImageSize->240,LabelStyle->20];
solw=Plot[Evaluate[ω[t]/.integration],{t,0,tmax},PlotRange-> All,AxesLabel-> {t,ω[t]},ImageSize->240,LabelStyle->20];
Row[{sola,solb,solw}]
pplot=ParametricPlot3D[Evaluate[{a[t],b[t],ω[t]}/.integration],{t,0,tmax},PlotRange-> All,BoxRatios->{1,1,1},AxesLabel-> {a[t],b[t],ω[t]},PlotStyle->Thin,PlotPoints->500,LabelStyle->Large];
Show[pplot,Gfp]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).