Quasi-periodic solution, coupled oscillations (Subscript[K, 1]),(Subscript[K, 2])≠0

Clear[tmax]
k1=2.0; (*coupling constants*)
k2=3.0;
ω1=5+(Sqrt[π]); (*ω natural frequency*)
ω2=1;(*2 then periodic (default 1)*)
tmax=3 π;

Style[StringForm["Ratio of fundamental frequencies ``",ω1/ω2] ,Red,17]
Style[StringForm["Is ω1/ω2 rational?: ``",ω1/ω2∈Rationals] ,Red,17]
kk=k1+k2;
Style[StringForm["k1 + k2 = ``",kk] ,Red,17]
ww=Abs[ω1-ω2]; (*if this is bigger than k1+k2 and k1,k2=/=0, then sol is quasiperiodic. Othervise periodic or periodic phase locking solution*)
Style[StringForm["Abs[ω1-ω2] = ``",ww] ,Red,17]
If[ww≥kk,
Print[Style["Solution is quasi-periodic (only if ω1/ω2 is irrational)",Red,17]],
Print[Style["Solution is periodic or long term sol'n is.",Blue,17]]]

If[ww≤ kk,
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20,Epilog->{Red,AbsolutePointSize[10],Point[{ϕ,0}/.Quiet[Solve[ω1-ω2-(k1+k2) Sin[ϕ]==0&&ϕ≥0&&ϕ<2 π,ϕ,Reals]]]}],
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20]]
dθ1=ω1+k1 Sin[θ2-θ1]; (*θ phase*)
dθ2=ω2+k2 Sin[θ1-θ2];
sp=StreamPlot[{dθ1,dθ2},{θ1,0,2 π},{θ2,0,2 π},FrameLabel->{θ1[t],θ2[t]},LabelStyle->20];
integ=NDSolve[{
θ1'[t]==ω1+k1 Sin[θ2[t]-θ1[t]],
θ2'[t]==ω2+k2 Sin[θ1[t]-θ2[t]],
θ1[0]==0,θ2[0]==1},
{θ1,θ2},{t,0,tmax}];
plot1=ParametricPlot[Evaluate[{θ1[t],θ2[t]}/.integ],{t,0,tmax},AxesLabel->{θ1[t],θ2[t]},ImageSize->260,PlotStyle->Red,LabelStyle->20];
plot2=Plot[Evaluate[θ1[t]/.integ],{t,0,tmax},AxesLabel-> {t,θ1[t]},PlotRange-> All,ImageSize->260,LabelStyle->20];
plot3=Plot[Evaluate[θ2[t]/.integ],{t,0,tmax},AxesLabel-> {t,θ2[t]},PlotRange-> All,ImageSize->260,LabelStyle->20];
Show[sp,plot1]
Row[{plot1,plot2,plot3}]
sinp1=Plot[Evaluate[Cos[θ1[t]]/.integ],{t,0,tmax},AxesLabel-> {t,Cos[θ1]},PlotRange-> All,ImageSize->420,LabelStyle->20];
sinp2=Plot[Evaluate[Sin[θ2[t]]/.integ],{t,0,tmax},AxesLabel-> {t,Sin[θ2]},PlotRange-> All,ImageSize->420,LabelStyle->20];
Row[{sinp1,sinp2}]
pplo1=ParametricPlot[Evaluate[{
Cos[θ1[t]],
D[Cos[θ1[t]],t]}/.integ],
{t,0,tmax},
AxesLabel-> {"cos θ1","(∂_t)cos θ1"},PlotLabel->"Phase portrait of cos θ1",PlotRange->{{-1.01,1.01},All},ImageSize->360,AspectRatio->1,LabelStyle->20];
pplo2=ParametricPlot[Evaluate[{
Sin[θ2[t]],
D[Sin[θ2[t]],t]}/.integ],
{t,0,tmax},
AxesLabel-> {"sin θ2","(∂_t)sin θ2"},PlotLabel->"Phase portrait of sin θ2",PlotRange->{{-1.01,1.01},All},ImageSize->360,AspectRatio->1,LabelStyle->20];
Row[{pplo1,pplo2}]

Power spectrum of the signals

Clear[tmax]
k1=2.0; (*coupling constants*)
k2=3.0;
ω1=5+(Sqrt[π]); (*ω natural frequency*)
ω2=1;(*2 then periodic (default 1)*)
tmax=10 π 2;

Style[StringForm["Ratio of fundamental frequencies ``",ω1/ω2] ,Red,17]
Style[StringForm["Is ω1/ω2 rational?: ``",ω1/ω2∈Rationals] ,Red,17]
kk=k1+k2;
Style[StringForm["k1 + k2 = ``",kk] ,Red,17]
ww=Abs[ω1-ω2]; (*if this is bigger than k1+k2 and k1,k2=/=0, then sol is quasiperiodic. Othervise periodic or periodic phase locking solution*)
Style[StringForm["Abs[ω1-ω2] = ``",ww] ,Red,17]
If[ww>=kk,
Print[Style["Solution is quasi-periodic (only if ω1/ω2 is irrational)",Red,17]],
Print[Style["Solution is periodic or long term sol'n is",Blue,17]]]

If[ww≤ kk,
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20,Epilog->{Red,AbsolutePointSize[10],Point[{ϕ,0}/.Quiet[Solve[ω1-ω2-(k1+k2) Sin[ϕ]==0&&ϕ≥0&&ϕ<2 π,ϕ,Reals]]]}],
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20]]

integ=NDSolve[{
θ1'[t]==ω1+k1 Sin[θ2[t]-θ1[t]],
θ2'[t]==ω2+k2 Sin[θ1[t]-θ2[t]],
θ1[0]==0,θ2[0]==0},
{θ1,θ2},{t,0,tmax}];
sinp1=Plot[Evaluate[Cos[θ1[t]]/.integ],{t,0,tmax},AxesLabel-> {t,Cos[θ1]},PlotRange-> All,ImageSize->380,LabelStyle->20];
sinp2=Plot[Evaluate[Sin[θ2[t]]/.integ],{t,0,tmax},AxesLabel-> {t,Sin[θ2]},PlotRange-> All,ImageSize->380,LabelStyle->20];

sr=11;(*Sampling Rate: keep Nyquist limit in mind*)
signal=Table[Evaluate[Cos[θ1[t]]/.integ[[1]]],{t,0,tmax,1/sr}];(*sample values*)
signal2=Table[Evaluate[Sin[θ2[t]]/.integ[[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"},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"},ImageSize->380,LabelStyle->20];
Row[{sinp1,specter1}]
Row[{sinp2,specter2}]

Phase differences, analysis of 1-D phase portrait

ϕ0=0.0; (*initial condition*)
k1=1.0; (*coupling constants*)
k2=2.0;
ω1=5.02; (*ω natural frequency*)
ω2=2.046766352;(*2.0467... periodic*)
tmax=15;

kk=k1+k2;
ww=Abs[ω1-ω2]; 

If[ww>=kk,
Print[Style["Solution is quasi-periodic (coupled cases only)",Red,17]],
Print[Style["Solution is periodic or long term sol'n is (coupled cases only)",Blue,17]]
]
integ2=NDSolve[{ϕ'[t]==ω1-ω2-(k1+k2) Sin[ϕ[t]],ϕ[0]==ϕ0},ϕ,{t,0,tmax}];
Plot[Evaluate[ϕ[t]/.integ2],{t,0,tmax},AxesLabel-> {t,ϕ[t]},PlotRange-> All,LabelStyle->20]

If[ww≤ kk,
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20,Epilog->{Red,AbsolutePointSize[10],Point[{ϕ,0}/.Quiet[Solve[ω1-ω2-(k1+k2) Sin[ϕ]==0&&ϕ≥0&&ϕ<2 π,ϕ,Reals]]]}],
Plot[ω1-ω2-(k1+k2) Sin[ϕ],{ϕ,0,2 π},AxesLabel->{ϕ,"(∂_t)ϕ"},LabelStyle->20]]

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