I am glad you asked this question because I have worked on this problem in my post grad years, as well as the 3D channel. There are two answers, the short one and the long one, I am going to give you the long one...
First of all there has to be a pressure gradient because it is equivalent to a body force exerted on the fluid and it is this what drives the flow. If there where no pressure gradient there are still solutions but they consist of vortices that decay exponentially in time and the flow eventually stops. Not an interesting situation.
What you have described is the "velocity profile" which is the steady-state of the system. Let me call it $\mathbf{U}$, ie
$$U_i=U(x_2) \delta_{i1}= \dfrac{F}{2\mu}(x_2^2-a^2) \delta_{i1}$$
where $F= (p_2-p_1)/L$ is the pressure gradient or body force. An easy calculation shows that
$$\mathbf{U}. \nabla \mathbf{U}=0$$
and
$$ -\dfrac{1}{\rho_0} \mathbf{F} + \frac{\mu}{\rho_0} \nabla^2\mathbf{U}=0$$
So it's a solution of the Navier Stokes and also $\nabla. \mathbf{U}=0$ which we need for incompressible flow.
So far so good, now we ask the question are there any other flows that satisfy incompressibility and N-S equations. Let me introduce a small deviation from the steady state, by substituting $\mathbf{U}+\mathbf{u}$ instead of $\mathbf{u}$ into the N-S equation (which is actually the momentum conservation equation):
$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}.\nabla \mathbf{U}+\mathbf{U}.\nabla \mathbf{u}+\mathbf{u}.\nabla \mathbf{u}=\nu \nabla^2 \mathbf{u}$$
Happily the force dropped out of the problem because it balanced out viscosity, $\nu\nabla^2\mathbf{U}$. But we have introduced two more convective terms:
$$\mathbf{u}.\nabla \mathbf{U}=u_2U' \delta_{i1}$$
$$\mathbf{U}.\nabla{u}=U \frac{\partial u_i}{ \partial x_1}$$
Both of these terms are linear with respect to $\mathbf{u}$ but we cannot get rid of the term $\mathbf{u}.\nabla \mathbf{u}$ which is nonlinear. The viscosity term is also linear in $\mathbf{u}$ so we can write the equation finally as
$$\frac{\partial \mathbf{u}}{\partial t}=-\mathbf{u}.\nabla \mathbf{u}+\mathbf{L} \mathbf{u}$$
where $\mathbf{L}$ is a linear operator acting on $\mathbf{u}$. If $\mathbf{u}$ is really tiny the quadratic term is one order of magnitude smaller than the linear terms, so we might as well solve
$$\frac{\partial \mathbf{u}}{\partial t}=\mathbf{L} \mathbf{u}$$
together with our usual incompressibility constraint
$$\nabla . \mathbf{u}=0$$
This is typically solved by the method of eigenvalues and eigenfunctions of the linear operator:
$$\mathbf{u}( \mathbf{x},t)=\sum c_n e^{\alpha_n t} \mathbf{\phi_n(\mathbf{x})}$$
The time dependence is the essential thing here. For small Reynolds numbers all the eigenvalues $\alpha_n$ have negative real part and so the flow soon enough returns to the steady state $\mathbf{U}$. But there is a threshold above which some eigenvalues have positive real part and so we have exponential growth! And soon enough the nonlinear term cannot be neglected any more. In this case we have chaotic turbulent flow. It's not like genuine 3D turbulence, but it's still a kind of turbulence. You cannot have that in the linear case.
Now, all this was the framework, I still haven't answered your question! And the answer applies to the majority of nonlinear systems which are chaotic: Suppose we Fourier expand the solution in terms of $y$:
$$\mathbf{u}(x,y,t)=\sum_{n=1}^N \mathbf{v}_n(x,t) \sin{( \frac{n \pi y}{a})}$$
respecting the boundary conditions at $y=\pm a$ and keeping a finite number of terms. Never mind about the $\cos$ for the moment, let's assume it's odd in $y$ and keeping only a finite number of modes $N$. In our equation
$$\frac{\partial \mathbf{u}}{\partial t}=\sum_{n=1}^N \frac{\partial \mathbf{v}_n}{\partial t} \sin{( \frac{n \pi y}{a})}$$
i.e. is also expressible as a sum of the same modes $n=1,2, \dots,N$. The same applies with the linear term $\mathbf{L} \mathbf{u}$. Now we come to the "juice", the nonlinear term $\mathbf{u}.\nabla \mathbf{u}$. That is going to produce modes greater than $N$ that cannot be balanced out since all the other terms in the equation involve modes up to $N$. For example if we had the first mode only say $f(y)=\sin(\frac{\pi y}{a})$, the quadratic term
$$f(y) f'(y)=\frac{\pi}{a} \sin(\frac{\pi y}{a}) \cos(\frac{\pi y}{a})=\frac{\pi}{2a}\sin(\frac{2 \pi y}{a}) $$
which involves mode 2. Similarly with the $\cos$. And in general any finite sum of modes will give you modes more than you had considered. This is the so called "energy cascade", the mechanism through which energy is transfered from larger scales to smaller scales. Big eddies break up and transfer their energy to smaller and smaller eddies until viscosity is the major factor, not the nonlinear terms, and then friction takes over and their energy becomes heat. You are doing work pumping the fluid down the channel, where is that work going, it becomes heat. The energy cascade is essential to energy balance. So there is an energy leak in any finite solution to the N-S equation. We introduce more energy than we can handle.
And this is another reason why the N-S equations are incomplete: you need to introduce the heat equation as well to get it all balanced properly, which is coupled with the momentum equation but that's another story... In theory you could consider the infinite sum, i.e. all modes to infinity but then how are you going to simulate an infinite system on the computer?
Hope this helps, I 've just summed up two year's hard work!
