Solution of ODE's with a Chebyshev polynomial spectral method and a standard Fourier Spectral method

This notebook has been written in Mathematica by 

Mark J. McCready
Professor and Chair of Chemical Engineering
University of Notre Dame
Notre Dame IN 46556
USA

Mark.J.McCready.1@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by whatever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, that if you use it, that this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version: 2/14/14

This notebook demonstrates the Orszag-tau spectral method for a simple ODE and an ODE eigenvalue problem.  It is intended as a first introduction to solving these simple problems with a spectral method.

Reference:  S. A. Orszag (1971) "Accurate solution of the Orr-Sommerfeld stability equation", Journal of Fluid Mechanics, 50 pp 689-703.

The coefficients of the algebraic equations are computed directly from the orthogonality property of the Chebyshev polynomials using integration.

Mathematica aside

I have used a mix of input notation to show how it can be done.   Also, some of the major equations are shown in input form and also in typeset form for clarity.  The typset input equations are usually marked as inactive and so are not read in by Mathematica .. but they could be.

I would enter the dynamic boundary condition from a key pad in input form as

dc=- γ D[η[x,t],{x,2}] -
		(Subscript[ρ, 1])(-D[(Subscript[ϕ, 1])[x,y,t],t]-U1 D[(Subscript[ϕ, 1])[x,y,t],x]-g η[x,t])+
			(Subscript[ρ, 2])(-D[(Subscript[ϕ, 2])[x,y,t],t]-
					U2 D[(Subscript[ϕ, 2])[x,y,t],x]-g η[x,t])

In Standard form, which is a little shorter, you would need the type set window to make this preactical.

dc=-(γ (∂_({x,2}))η[x,t])-(Subscript[ρ, 1]) (-(∂_t)(Subscript[ϕ, 1])[x,y,t]-U1 (∂_x)(Subscript[ϕ, 1])[x,y,t]-g η[x,t])+(Subscript[ρ, 2]) (-(∂_t)(Subscript[ϕ, 2])[x,y,t]-U2 (∂_x)(Subscript[ϕ, 2])[x,y,t]-g η[x,t])

Finally it may be more easily read in Traditional form as

dc=-(γ (((∂^2)η(x,t))/(∂(x^2))))-(Subscript[ρ, 1]) (-((∂(Subscript[ϕ, 1])[x,y,t])/(∂t))-U1 ((∂(Subscript[ϕ, 1])[x,y,t])/(∂x))-g η(x,t))+(Subscript[ρ, 2]) (-((∂(Subscript[ϕ, 2])[x,y,t])/(∂t))-U2 ((∂(Subscript[ϕ, 2])[x,y,t])/(∂x))-g η(x,t))

Note that you can convert any "cryptic" input expression by selecting the cell, going up to the cell menu and selecting "convert to traditional form".  This is "shift curly t" on the keypad.

Some properties of Chebyshev polynomials

For this example, we will be using the ChebyshevT polynomials which are orthogonal and 
non singular on the interval [-1,1].  They also have the value of T[i,1]=T[i,-1]=±1.  
Let's see what these are

somechebys=Table[ChebyshevT[i,y],{i,0,10}]

If we wish, we can plot all of these.  
The Evaluate is needed to make the function be recognized before it is numerically evaluated.

Plot[Evaluate[somechebys],{y,-1,1}]

This is too complicated to see much.  But there are even polynomials and odd ones with the expected symmetry.
Also, it can be seen that the interval is spanned effectively with these so they look to be good exapnsion functions.

somechebys=Table[ChebyshevT[i,y],{i,0,10,2}]

If we wish, we can plot all of these.  
The Evaluate is needed to make the function be recognized before it is numerically evaluated.

Plot[Evaluate[somechebys],{y,-1,1}]

These are orthogonal, let's see how

(∫_(-1)^1)ChebyshevT[2,y] ChebyshevT[1,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y

(((∫_(-1)^1)ChebyshevT[0,y] ChebyshevT[0,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/π)

((((∫_(-1)^1)ChebyshevT[1,y] ChebyshevT[1,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y) 2)/π)

((((∫_(-1)^1)ChebyshevT[2,y] ChebyshevT[2,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y) 2)/π)

((((∫_(-1)^1)ChebyshevT[0,y] ChebyshevT[4,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y) 2)/π)

((((∫_(-1)^1)ChebyshevT[3,y] ChebyshevT[6,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y) 2)/π)

((((∫_(-1)^1)ChebyshevT[5,y] ChebyshevT[5,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y) 2)/π)

I hope that we see that pattern.  The weight function is 1/sqrt[1-y^2] 2/πi/c[i].
  The c[i]=1 for all i except i=0 when it equals 2.

ChebyshevT[0,y]

ChebyshevT[1,y]

ChebyshevT[2,y]

ChebyshevT[4,y]

((2 (∫_(-1)^1)(1-(y^2)) ChebyshevT[0,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/π)

((2 (∫_(-1)^1)(1-(y^2)) ChebyshevT[2,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/π)

((2 (∫_(-1)^1)(1-(y^2)) ChebyshevT[3,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/π)

(∫_(-1)^1)((ChebyshevT[2,y])^2) (1/(Sqrt[1-(y^2)]))\[DifferentialD]y

ChebyshevT[0,y]-ChebyshevT[2,y]

((2 (∫_(-1)^1)((ChebyshevT[3,y])^2) (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/π)

(∫_(-1)^1)((ChebyshevT[0,y])^2) (1/(Sqrt[1-(y^2)]))\[DifferentialD]y

What about some derivatives.

D[ChebyshevT[i,y],{y,1}]

D[ChebyshevT[i,y],{y,2}]

D[ChebyshevT[i,y],{y,3}]

D[ChebyshevT[i,y],{y,4}]

Basic idea of how to use a spectral numerical method to solve an ODE.

We have seen that we can fit an arbitrary function using a set of orthogonal polynomials by forming the inner product to evaluate the coefficients.  We thus conjecture that we could solve an ode, which must have some function form (albeit) unknown, with an appropriate set of orthogonal polynomials, or other orthogonal functions.  The obvious questions are which functions and how do we get the cofficients.

The starting point is an assumed solution form of:

u[x]  =(∑_(j=1)^∞)(Subscript[c, j]) (Subscript[ϕ, j])(x)

We need the (Subscript[ϕ, j])'s to be orthogonal on the domain of interest and capable of producing a boundary conditions of the needed form.  

We will substitute the assumed form for the solution into the ODE and boundary conditions, form an inner product with every term (not the boundary conditions however) (why not), and then see if we can evaluate the coefficients.  If we can, we have a solution.

Note however, that we will expect to use a finite number of terms -- as opposed to the series solution of an ode where we expect to find a recursion relation valid for all terms.  We then hope that the solution will get more accurate as we increase the number of terms.

We first start with a very simple ode that has an analytical solution with which to compare.

Solution of a simple ode

Here is a simple ODE that has an easy analytical solution.  We will then use a spectral technique and then compare the answers. 

This is a case where the viscosity of a fluid in a pressure driven flow increases or decreases away from the centerline.  This model for viscosity is not in any way physically defensible. but it shows the solution technique.

((d (α0 + α2 y^2)((d u[y])/dy))/dy) - (dP/dx)=0, 
 u[1]=u[-1]=0

Analytical solution

ode=((∂((α0 + α2 y^2) ((∂ u(y))/(∂y))))/(∂y)) -dpdx

ans1=DSolve[{ode==0,u[1]==0,u[-1]==0},u[y],y]

soln1=Simplify[u[y]/.ans1[[1]]]

p1=Plot[soln1/.{α0->.01,α2->.0001,dpdx->-.01},{y,-1,1}]

Galerkin Fourier mode solution

Now solve the equation using Cosines as the basis functions.

ode=((∂((α0 + α2 y^2) ((∂ vv(y))/(∂y))))/(∂y)) -dpdx

The basic assumption of a Galerkin spectral method is that the numerical solution, can be written as a sum of orthogonal functions.  This trial solution has a finite number of terms and is valid over the entire domain.  If the functions are chosen to automatically fit the boundary conditions, then we just need to worry about transforming the differential equation.

 Here we try Cosines

Plot[{Cos[π/2 y],Cos[π y],Cos[3/2 π y]},{y,-1,1}]

Integrate[Cos[π/2 y]^2,{y,-1,1}]

Integrate[Cos[π/2 y]  Cos[3 π/2 y],{y,-1,1}]

If we choose Cos[ i/2 y] for i = 1,n (i=odd), we will have the boundary conditions fit and you can see from the “1”, the ability to produce a non-zero, always positive solution

Expand the homogeneous part of the differential equation.

We need to keep the inhomogeneous part separate and add it in after the expansions...

ode1=((∂((α0 + α2 y^2) ((∂ vv(y))/(∂y))))/(∂y))

odeinhomg = -dpdx

temp1=ode1/.{vv[y]->a[i] Cos[π/2 i y],(∂_({y,a1_}))vv[y]->(∂_({y,a1}))(a[i] Cos[π/2 i y])}

temp2=ExpandAll[temp1]

Now get a finite size sum to solve.   We will be doing only the odd terms, so the number of “modes” will be n/2.

n=12;

The Δi=2, means to skip every other term.

temp3=(∑_({{i=1}, {Δ i=2}})^n)temp2

This might go quicker if we simplify it.

temp4=Simplify[temp3]

We need to find algebraic equations for the a[i]'s.  We will use the orthonality property of the functions to pick out equations for each a[i]. 

(∫_(-1)^1)  ρ(y)*(Subscript[ϕ, j])(y)*(substituted ODE with a number of (Subscript[ϕ, i])(y)'sdy . 

For Sines and Cosines, the weighting functio,  ρ(y) =1.  
  
Now set up integrals in the  n/2 equations for the n/2 coefficients.  We note that the integration will take of itself

temp5=(temp4+odeinhomg)  Cos[π/2 j y]

Expand[temp5]

Integrate[ Cos[ π/2 y] Cos[π/2 y],{y,-1,1}]

Integrate[ Cos[ π/2 y] Cos[ 5π/2 y],{y,-1,1}]

Here we generate the n/2 equations in a table form.  The original equation was equal to 0, now each of these pieces is also equal to 0.

temp6=Table[temp5,{j,1,n,2}]

Now we complete the calculation by integrating "temp6" from -1 to 1.

temp7=Integrate[temp6,{y,-1,1}]

We can solve at this point for the coefficients.  If the functions did not “automatically” fit the boundary conditions we would be doing additional work at this point using the “Tau” method.

temp8=Table[temp7[[i]]==0,{i,1,(n/2)}]

We solve for the coefficients.

coef=Solve[temp8,Table[a[i],{i,1,n,2}]]

We want to reconstruct the solution.  First make the series.

expand1=(∑_({{i=1}, {Δ i=2}})^n)a[i]  Cos[π/2 i y]

Then substitute the coefficients that we calculated above to get the answer.

anscos=expand1/.coef[[1]]

Limit[anscos,α2->0]

Here is a plot of the numerical solution.

p2=Plot[anscos/.{α0->.01,α2->-.007,dpdx->-.01},{y,-1,1}]

Check with the exact solution.

p1=Plot[soln1/.{α0->.01,α2->-.007,dpdx->-.01},{y,-1,1}]

Show[p1,p2]

These are almost the same!

Look at the series expansions.  First the exact solution:

analytseries=Chop[N[Series[soln1/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the numerical solution:

numseries=Chop[N[Series[anscos/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the difference.  The numerical solution works, but only by cancelation of +/- terms near to the endpoints.

comp1=analytseries-numseries

Normal[(comp1)]

Here is an integral measure of the error.  If we use enough terms, it is difficult for Mathematica to compute the error.

e[n]=NIntegrate[Normal[(comp1)],{y,-1,1}]

This approach will work for any n, but will get real slow as we increase n.  Thus the next section employs the same technique but with a little more simplifications before doing the integrals.  This will allow reasonable processing times for larger n's.

Chebyshev Numerical solution

Now we use Chebyshev polynomials as the basis functions.  We use direct substitution and don't attempt to reduce the derivatives with analytical manipulations (If you want high accuracy, you will want to find recursion relations for the coefficients.)

ode=((∂((α0 + α2 y^2) ((∂ vv(y))/(∂y))))/(∂y)) -dpdx

We need to keep the inhomogeneous part separate and add it in after the expansions...

ode1=((∂((α0 + α2 y^2) ((∂ vv(y))/(∂y))))/(∂y))

odeinhomg = -dpdx

The basic assumption of a Galerkin spectral method is that the numerical solution, can be written as a sum of orthogonal functions.  This trial solution has a finite number of terms and is valid over the entire domain.  If the functions are chosen to automatically fit the boundary conditions (e.g., perhaps by some clever choices of sums of groups of polynomials, or shear luck).

 If the functions do not assure automatic fitting of the boundary conditions, we will have to make some adjustments.  The tau method, used below, replaces one of the equations for a[i] generated by the ODE, for each algebraic equation that is needed to assure boundary condition fitting.  Stated again, if we have 4 boundary conditions, we will replace 4 of the equations generated from the ode with the 4 equations that must be fit on the boundary.  To minimize the error, we think of this as starting with more terms in the sum than we need, then dropping of the highest "frequency" modes that are assumed (this is confirmed) to be less important to getting an accurate solution than the lowest modes.

Here is where we substitute Chebyshev polynomials for vv[y].

temp1=ode1/.{vv[y]->a[i] ChebyshevT[i,y],(∂_({y,a1_}))vv[y]->(∂_({y,a1}))(a[i] ChebyshevT[i,y])}

temp2=ExpandAll[temp1]

Now get a finite size sum to solve.

n=8;

Because of the symmetry of the boundary conditions and the occurence of only even derivatives, the solution will be an even function so we need only even polynomials.  The Δi=2, means to skip every other term.

temp3=(∑_({{i=0}, {Δ i=2}})^n)temp2

This might go quicker if we simplify it.

temp4=Simplify[temp3]

We need to find algebraic equations for the a[i]'s.  We will use the orthonality property of the polynomials to pick out equations for each a[i]. 

(∫_(-1)^1)  ρ(y)*(Subscript[T, j])(y)*(substituted ODE with a number of (Subscript[T, i])(y)'sdy . 

One other "trick" needs to be recognized.  The orthogonality property will take care of "choosing" the nonzero pieces of the equation for each coefficient even if we multiply everthing out ahead of time as we have done above.

  
Now set up integrals noting the complete tricks needed for orthogonality.  We multiply by the weighting function and an arbitrary CheybyshevT.  We will let the "i" take on all possible values between 0 and n to form n/2 equations for the n/2 coefficients.  We note that the integration will take of itself

temp5=(((temp4+odeinhomg) (1/(Sqrt[1-(y^2)])) ChebyshevT[j,y] 2)/(π c[j]))

Here we generate the n/2 equations in a table form.  The original equation was equal to 0, now each of these pieces is also equal to 0.

temp6=Table[temp5,{j,0,n,2}]

Integrate[ dpdx/Sqrt[1-y^2] 1,{y,-1,1}]

Integrate[ dpdx/Sqrt[1-y^2] (-1+2 (y^2)),{y,-1,1}]

Let's get rid of the c[j]'s

temp7=(temp6/.c[0]->2)/.Table[c[i]->1,{i,2,n,2}]

temp7[[1]]

temp7[[2]]

temp7[[3]]

Now we complete the calcultion by integrating "temp7" from -1 to 1.

temp8=Integrate[temp7,{y,-1,1}]

temp8[[1]]

Each of the entries in this table is an equation that equals 0.  We can solve this at this point at get the trivial solution since there are 4 unknown a(j)'s and 4 equations.  Of course, since we have thus far just substituted into a homogeneous equation and not yet used the nonhomogeneous boundary conditions, what other answer could we expect!  

So we need to enforce the nonhomogeneous boundary conditions both to get a nonzero answer and also to get the correct answer.  

As mentioned above, we will replace one of these n/2+1 equations with the boundary condition.  You might ask which one should we replace.  The obvious one(s) to replace are the equations generated by the highest values of i.  We justify this two ways.  One is that we employ extra terms, which are necessarily at high values of "i" in the expansion so that there are enough to get a good answer.  Thus dropping a few at high i does not seem to hurt.  A second way to think of this is that the functions have an intrinsic "frequency" and we are just dropping off some of the higher frequencies that are not important.

What is the boundary condition in terms of the Chebyshev's?

We need u[1]=u[-1]=0

bc1= (∑_({{i=0}, {Δ i=2}})^n)a[i]

So we will have n/2 equations that =0 and one (bc1), that equals 1 in our set to solve.

Here is where we take the first n/2 rows of the coefficients from the ODE and Append the one boundary condition.

temp9=Append[Table[temp8[[i]]==0,{i,1,(n/2)}],bc1==0]

We solve for the coefficients.

coef=Solve[temp9,Table[a[i],{i,0,n,2}]]

We want to reconstruct the solution.  First make the series.

expand1=(∑_({{i=0}, {Δ i=2}})^n)a[i] ChebyshevT[i,y]

Then substitute the coefficients that we calculated above to get the answer.

anscheby=expand1/.coef[[1]]

FullSimplify[anscheby]

%/.α2->0

FullSimplify[soln1]

Limit[soln1,α2->0]

Here is a plot of the numerical solution.

p1=Plot[soln1/.{α0->.01,α2->.0001,dpdx->-.01},{y,-1,1}]

p2=Plot[anscheby/.{α0->.01,α2->.0001,dpdx->-.01},{y,-1,1}]

Check with the exact solution, we can't see any difference.

p1=Plot[soln1/.{α0->.01,α2->-.007,dpdx->-.01},{y,-1,1}]

Show[p1,p2]

Look at the series expansions.  First the exact solution:

analytseries=Chop[N[Series[soln1/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the numerical solution:

numseries=Chop[N[Series[anscheby/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the difference.

comp1=analytseries-numseries

Normal[(comp1)]

Here is an integral measure of the error.  If we use enough terms, it is difficult for Mathematica to compute the error.

e[n]=NIntegrate[Normal[(comp1)],{y,-1,1}]

This approach will work for any n, but will get real slow as we increase n.  Thus the next section employs the same technique but with a little more simplifications before doing the integrals.  This will allow reasonable processing times for larger n's.

Analytical solution

Plot[ (α0 + α2 y^2)/.{α0->.01, α2->-.007},{y,-1,1}]

Plot[ (α0 Cosh[π/2] - α2 Cosh[π/2 y])/.{α0->.006, α2->.005},{y,-1,1}]

Show[%,%517,PlotRange->All]

ode=((∂((α0 Cosh[π/2] - α2 Cosh[π/2 y])((∂ u(y))/(∂y))))/(∂y)) -dpdx

ans1=DSolve[{ode==0,u[1]==0,u[-1]==0},u[y],y]

soln1=Simplify[u[y]/.ans1[[1]]]

p1=Plot[soln1/.{α0->.006, α2->.005,dpdx->-.01},{y,-1,1}]

Numerical solution

Now solve the equation with a Galerkin (name of the person associated with this kind of spectral approach) spectral technique that involves substituting a series of Chebyshevs into the differential equation, employing orthogonality and then solving the algebraic equations. We use direct substitution and don't attempt to reduce the derivatives with analytical manipulations.

ode=((∂((α0 Cosh[π/2] - α2 Cosh[π/2 y]) ((∂ vv(y))/(∂y))))/(∂y)) -dpdx

We need to keep the inhomogeneous part separate and add it in after the expansions...

ode1=((∂((α0 Cosh[π/2] - α2 Cosh[π/2 y])((∂ vv(y))/(∂y))))/(∂y))

odeinhomg = -dpdx

The basic assumption of a Galerkin spectral method is that the numerical solution, can be written as a sum of orthogonal functions.  This trial solution has a finite number of terms and is valid over the entire domain.  If the functions are chosen to automatically fit the boundary conditions (e.g., perhaps by some clever choices of sums of groups of polynomials, or shear luck).

 If the functions do not assure automatic fitting of the boundary conditions, we will have to make some adjustments.  The tau method, used below, replaces one of the equations for a[i] generated by the ODE, for each algebraic equation that is needed to assure boundary condition fitting.  Stated again, if we have 4 boundary conditions, we will replace 4 of the equations generated from the ode with the 4 equations that must be fit on the boundary.  To minimize the error, we think of this as starting with more terms in the sum than we need, then dropping of the highest "frequency" modes that are assumed (this is confirmed) to be less important to getting an accurate solution than the lowest modes.

Here is where we substitute Chebyshev polynomials for vv[y].

temp1=ode1/.{vv[y]->a[i] ChebyshevT[i,y],(∂_({y,a1_}))vv[y]->(∂_({y,a1}))(a[i] ChebyshevT[i,y])}

temp2=ExpandAll[temp1]

Now get a finite size sum to solve.

n=8;

Because of the symmetry of the boundary conditions and the occurence of only even derivatives, the solution will be an even function so we need only even polynomials.  The Δi=2, means to skip every other term.

temp3=(∑_({{i=0}, {Δ i=2}})^n)temp2

This might go quicker if we simplify it.

temp4=Simplify[temp3]

We need to find algebraic equations for the a[i]'s.  We will use the orthonality property of the polynomials to pick out equations for each a[i]. 

(∫_(-1)^1)  ρ(y)*(Subscript[T, j])(y)*(substituted ODE with a number of (Subscript[T, i])(y)'sdy . 

One other "trick" needs to be recognized.  The orthogonality property will take care of "choosing" the nonzero pieces of the equation for each coefficient even if we multiply everthing out ahead of time as we have done above.

  
Now set up integrals noting the complete tricks needed for orthogonality.  We multiply by the weighting function and an arbitrary CheybyshevT.  We will let the "i" take on all possible values between 0 and n to form n/2 equations for the n/2 coefficients.  We note that the integration will take of itself

temp5=(((temp4+odeinhomg) (1/(Sqrt[1-(y^2)])) ChebyshevT[j,y] 2)/(π c[j]))

Here we generate the n/2 equations in a table form.  The original equation was equal to 0, now each of these pieces is also equal to 0.

temp6=Table[temp5,{j,0,n,2}]

Let's get rid of the c[j]'s

temp7=(temp6/.c[0]->2)/.Table[c[i]->1,{i,2,n,2}]

Now we complete the calcultion by integrating "temp7" from -1 to 1.

temp8=Integrate[temp7,{y,-1,1}]

Each of the entries in this table is an equation that equals 0.  We can solve this at this point at get the trivial solution since there are 4 unknown a(j)'s and 4 equations.  Of course, since we have thus far just substituted into a homogeneous equation and not yet used the nonhomogeneous boundary conditions, what other answer could we expect!  

So we need to enforce the nonhomogeneous boundary conditions both to get a nonzero answer and also to get the correct answer.  

As mentioned above, we will replace one of these n/2+1 equations with the boundary condition.  You might ask which one should we replace.  The obvious one(s) to replace are the equations generated by the highest values of i.  We justify this two ways.  One is that we employ extra terms, which are necessarily at high values of "i" in the expansion so that there are enough to get a good answer.  Thus dropping a few at high i does not seem to hurt.  A second way to think of this is that the functions have an intrinsic "frequency" and we are just dropping off some of the higher frequencies that are not important.

What is the boundary condition in terms of the Chebyshev's?

We need u[1]=u[-1]=0

bc1= (∑_({{i=0}, {Δ i=2}})^n)a[i]

So we will have n/2 equations that =0 and one (bc1), that equals 1 in our set to solve.

Here is where we take the first n/2 rows of the coefficients from the ODE and Append the one boundary condition.

temp9=Append[Table[temp8[[i]]==0,{i,1,(n/2)}],bc1==0]

We solve for the coefficients.

coef=Solve[temp9,Table[a[i],{i,0,n,2}]]

We want to reconstruct the solution.  First make the series.

expand1=(∑_({{i=0}, {Δ i=2}})^n)a[i] ChebyshevT[i,y]

Then substitute the coefficients that we calculated above to get the answer.

anscheby=expand1/.coef[[1]]

FullSimplify[anscheby]

%/.α2->0

FullSimplify[soln1]

Limit[soln1,α2->0]

Here is a plot of the numerical solution.

p2=Plot[anscheby/.{α0->.006, α2->.005,dpdx->-.01},{y,-1,1}]

Check with the exact solution, we can't see any difference.

p1=Plot[soln1/.{α0->.006, α2->.005,dpdx->-.01},{y,-1,1}]

Show[p1,p2]

Look at the series expansions.  First the exact solution:

analytseries=Chop[N[Series[soln1/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the numerical solution:

numseries=Chop[N[Series[anscheby/.{α0->.01,α2->.0001,dpdx->-.01},{y,0,n}]]]

Here is the difference.

comp1=analytseries-numseries

Normal[(comp1)]

Here is an integral measure of the error.  If we use enough terms, it is difficult for Mathematica to compute the error.

e[n]=NIntegrate[Normal[(comp1)],{y,-1,1}]

This approach will work for any n, but will get real slow as we increase n.  Thus the next section employs the same technique but with a little more simplifications before doing the integrals.  This will allow reasonable processing times for larger n's.

Solution of the same simple ode using a slightly more sophisticated approach

Spectral numerical solution

The same ODE.

ode1=(((∂^2)vv(y))/(∂(y^2)))+vv(y)

The basic assumption of a Galerkin spectral method is that the numerical solution, can be written as a sum of orthogonal functions.  This trial solution has a finite number of terms and is valid over the entire domain.  If the functions are chosen to automatically fit the boundary conditions (e.g., perhaps by some clever choices of sums of groups of polynomials, or shear luck).

 If the functions do not assure automatic fitting of the boundary conditions, we will have to make some adjustments.  The tau method, used below, replaces one of the equations for a[i] generated by the ODE, for each algebraic equation that is needed to assure boundary condition fitting.  Stated again, if we have 4 boundary conditions, we will replace 4 of the equations generated from the ode with the 4 equations that must be fit on the boundary.  To minimize the error, we think of this as starting with more terms in the sum than we need, then dropping of the highest "frequency" modes that are assumed (this is confirmed) to be less important to getting an accurate solution than the lowest modes.

Here is where we substitute Chebyshev polynomials for vv[y].

temp1=ode1/.{vv[y]->a[i] ChebyshevT[i,y],(∂_({y,a1_}))vv[y]->(∂_({y,a1}))(a[i] ChebyshevT[i,y])}

temp2=ExpandAll[temp1]

Now get a finite size sum to solve.

n=20;

Because of the symmetry of the boundary conditions and the occurence of only even derivatives, the solution will be an even function so we need only even polynomials.  The Δi=2, means to skip every other term.

temp3=(∑_({{i=0}, {Δ i=2}})^n)temp2

This might go quicker if we simplify it.

temp4=Simplify[temp3]

We need to find algebraic equations for the a[i]'s.  We will use the orthonality property of the polynomials to pick out equations for each a[i]. 

(∫_(-1)^1)  ρ(y)*(Subscript[T, j])(y)*(substituted ODE with a number of (Subscript[T, i])(y)'sdy . 

One other "trick" needs to be recognized.  The orthogonality property will take care of "choosing" the nonzero pieces of the equation for each coefficient even if we multiply everthing out ahead of time as we have done above.

  
Now set up integrals noting the complete tricks needed for orthogonality.  We multiply by the weighting function and an arbitrary CheybyshevT.  We will let the "i" take on all possible values between 0 and n to form n/2 equations for the n/2 coefficients.  We note that the integration will take of itself

temp5=((temp4 (1/(Sqrt[1-(y^2)])) ChebyshevT[j,y] 2)/(π c[j]))

Here we generate the n/2 equations in a table form.  The original equation was equal to 0, now each of these pieces is also equal to 0.

temp6=Table[temp5,{j,0,n,2}]

Let's get rid of the c[j]'s

temp7=(temp6/.c[0]->2)/.Table[c[j]->1,{i,2,n,2}]

At this point, we could just integrate "temp7" from -1 to 1 and we would have the desired algebraic equations.  However the integration is very slow for a lot of terms and we need to speed it up.  To do this, we factor out the common terms.  These are all of the form  (y^(2n))/ √(1-(y^2)).  We first collect the terms in powers of y.  We use only the numerator for this so that the collect works.

temp8=Collect[Expand[Numerator[temp7]],y]

Here is where we factor them out.  The Transpose is needed to put the matrix back in the form that would multiply a column vector of the   (y^(2n))/ √(1-(y^2))'s.

temp9=Transpose[Table[Coefficient[temp8,y,i],{i,0,2 n,2}]]

Here is what we integrate:

temp10=Table[((y^i)/(π (Sqrt[1-(y^2)]))),{i,0,2 n,2}]

Here is the integration showing the time it takes.

Timing[temp11=(∫_(-1)^1)temp10\[DifferentialD]y]

Now multiply the coefficients times the integrated polynomial pieces.

temp12= Simplify[Expand[temp9.temp11]]

We are using the Tau method.  We  have to drop off a number of the highest modes equal to the number of boundary conditions.  Since the solution is symmetric we have already implicitly made the two boundary conditions equal by using only even polynomials.  Thus need to drop off only one the last one.   Note that on the boundaries, each polynomial has a value of 1 thus we write the boundary condition as

bc1= (∑_({{i=0}, {Δ i=2}})^n)a[i]

Here is where we take the first n/2 rows of the coefficients from the ODE and Append the one boundary condition.  Note how the value of the boundary condition was included while allowing "LogicalExpand" to do its trick.  This simple command produces a set of equations that can be solved.

temp13=LogicalExpand[Append[Table[temp12[[i]],{i,1,(n/2)}],bc1-1]==0]

We solve for the coefficients.

coef=Solve[temp13,Table[a[i],{i,0,n,2}]]

We want to reconstruct the solution.  First make the series.

expand1=(∑_({{i=0}, {Δ i=2}})^n)a[i] ChebyshevT[i,y]

Then substitute the coefficients that we calculated above to get the answer.

anscheby=expand1/.coef[[1]]

Here is a plot of the numerical solution.

p2=Plot[anscheby,{y,-1,1}]

Check with the exact solution, we can't see any difference.

Show[p1,p2]

Look at the series expansions.  First the exact solution:

analytseries=Chop[N[Series[soln1,{y,0,n}]]]

Here is the numerical solution:

numseries=N[Series[anscheby,{y,0,n}]]

Here is the difference.

comp1=analytseries-numseries

Here is an integral measure of the error.  If we use enough terms, it is difficult for Mathematica to compute the error.

e[n]=NIntegrate[(soln1-anscheby),{y,-1,1}]

Here are some answers from previous expansions and we see that the answer is getting better really quick

e[2]

e[4]

e[6]

e[8]

e[10]

e[12]

e[14]

e[16]

Evaluation of convergence

Here we have some saved results.

testtable=Table[{i,e[i]},{i,2,16,2}]

We need to load the graphics package to do log plots.

<<Graphics`Graphics`

errorplot=LinearLogListPlot[Abs[testtable],PlotJoined->True,
	AxesLabel->{"number of terms","error"}]

We see that we are getting exponential convergence with an increase in the number of terms, at least until we hit machine precision.

Numerical solution using just a power series

Now solve the equation by assuming that the solution is just a power series.  This will be substituted into the ODE and we will generate algebraic equations.  We will equate the coefficients for each power of y to 0 to generate algebraic equations.

Here, again, is the ODE.

ode1=(∂_({y,2}))vv[y] + vv[y]

n=12;

Here is our presumed solution.  Note that it is only even powers because we expect it to be symmetric.

seriesexpand=(∑_({{i=0}, {Δ i=2}})^n)a[i] (y^i)

Here we substitute.

temp1=ode1/.{vv[y]->seriesexpand,(∂_({y,a1_}))vv[y]->(∂_({y,a1}))seriesexpand}

Expand so that we can collect the individual powers of y.

temp2=ExpandAll[temp1]

Here we collect.

temp3=Collect[Expand[temp2],y]

This generates a table of the coefficients.

eqs=Table[Coefficient[temp3,y,i],{i,0,n,2}]

We need one boundary condition as above.

bcx1=seriesexpand/.y->1

This is how to make the equations for the coefficients.  Note that we replace the highest power with the boundary condition.

temp10=LogicalExpand[Append[Table[eqs[[i]],{i,1,n/2}],bcx1-1]==0]

We solve for the coefficients.

coef=Solve[temp10,Table[a[i],{i,0,n,2}]]

We want to reconstruct the solution.

ansx=seriesexpand/.coef[[1]]

Look at the series expansions.

analytseries=N[Series[soln1,{y,0,n}]]

numseries=N[Series[ansx,{y,0,n}]]

comp1=analytseries-numseries

Here is an integral measure of the error.

ep[n]=NIntegrate[(soln1-ansx),{y,-1,1}]

Here are some comparison of error for the Chebyshev and power series solutions.

Here are some answers from the Chebyshev spectral solutions

e[2]

e[4]

e[6]

e[8]

e[10]

Here are the power series solutions

ep[2]

ep[4]

ep[6]

ep[8]

ep[10]

ep[12]

testpowertable=Table[{i,ep[i]},{i,2,12,2}]

errorpowerplot=LinearLogListPlot[Abs[testpowertable],PlotJoined->True,
	AxesLabel->{"number of terms","error"}]

Show[errorplot,errorpowerplot]

We see that the power series does not converge nearly as fast at the Chebyshev expansion.  However, since it does not involve an integration, it runs faster than the Chebyshev version.  This is why people use formulas to obtained from various Trig identities to generate the algebraic equations for all of the derivatives.

For complex problems you would need exponential convergence and the power series method would not be a valid approach.

Eigenvalue problem

We will solve the ODE, with an explicit parameter,λ , for homogeneous boundary conditions.

u''[y]+(1+λ) u[y]=0,u[1]=0, u[-1]=0

Analytical solution to eigenvalue problem

Here we get the answer analytically.

When we solve this in Mathematica, we don't evaluate both boundary conditions or it will return the homogeneous solution.

DSolve[{u''[y]+(1+λ) u[y]==0,u[1]==0},u[y],y]

We can find a simpler even function answer:

anseig=A Cos[(Sqrt[1+λ]) y]

Check this answer

(∂_({y,2}))anseig+(1+λ) anseig

anseig/.y->-1

anseig/.y->1

We can see the eigenvalues from a plot as these are the 0's of the function, A Cos[(Sqrt[1+λ]) y].

Plot[anseig/.{A->1,y->1},{λ,0,30}]

Plot[anseig/.{A->1,y->1},{λ,0,240}]

The roots are easy to figure out:

FindRoot[(anseig/.{A->1,y->1})==0,{λ,2}]

We should be able to figure out the roots, as it is a simple equation!!

(Sqrt[2.4674])

N[(π/2)]

FindRoot[(anseig/.{A->1,y->1})==0,{λ,20}]

Here are some eigenvalues

eigsanlyt = Table[N[((i*Pi)/2)^2 - 1], {i, 1, 13, 2}]

Here are some plots of eigenfunctions

λ = 1.4674

afunc1=(anseig/.A->1)/.λ->eigsanlyt[[1]]

aplot1=Plot[afunc1/.A->1,{y,-1,1}]

λ = 21.2066

afunc2=(anseig/.A->1)/.λ->eigsanlyt[[2]]

aplot2=Plot[afunc2/.A->1,{y,-1,1}]

λ = 60.7158

afunc3=(anseig/.A->1)/.λ->eigsanlyt[[3]]

aplot3=Plot[afunc3/.A->1,{y,-1,1}]

λ = 119.903

afunc4=(anseig/.A->1)/.λ->eigsanlyt[[4]]

aplot4=Plot[afunc4/.A->1,{y,-1,1}]

Show[aplot1,aplot2,aplot3,aplot4]

Here is the series expansion of the answer

Series[anseig/.{A->1},{y,0,8}]

Numerical Solution to the eigenvalue problem

Numerical solution

Here is our equation again.

odeeig1=(((∂^2)vv(y))/(∂(y^2)))+(λ+1) vv(y)

Here is the standard substitution.

t1 = odeeig1 /. 
   {vv[y] -> a[i]*ChebyshevT[i, y], 
    D[vv[y], {y, a1_}] -> D[a[i]*ChebyshevT[i, y], {y, a1}]}

t2=ExpandAll[t1]

Now get a finite size sum to solve.  Note that we will need only even terms because the solution is clearly even.

n=12;

t3=(∑_({{i=0}, {Δ i=2}})^n)t2

This might go quicker if we simplify it.

t4=Simplify[t3]

We need to find equations for the a[i]'s.  We can use the orthonality property of the polynomials to pick out equations for each a[i].   
Set up integrals

t5=((t4 (1/(Sqrt[1-(y^2)])) ChebyshevT[i,y] 2)/(π c[i]))

t6=Table[t5,{i,0,n,2}]

t7=(t6/.c[0]->2)/.Table[c[i]->1,{i,2,n}]

t8=Collect[Expand[Numerator[t7]],y]

t9=Transpose[Table[Coefficient[t8,y,i],{i,0,2 n,2}]]

t10=Table[((y^i)/(π (Sqrt[1-(y^2)]))),{i,0,2 n,2}]

Timing[t11=(∫_(-1)^1)t10\[DifferentialD]y]

t12=Simplify[Expand[t9.t11]]

Here is the eigenvalue problem from the ODE.  At this point, we have not enforced the boundary conditions so we cannot just solve this matrix for λ.

t13=Transpose[Table[Coefficient[t12,a[i]],{i,0,n,2}]]

Here is the homogeneous boundary condition at y=1.

b1=Table[1,{i,0,n,2}]

To enforce the boundary condition, we drop off the last equation (highest mode) and replace it with the boundary condition.  If course, this makes the eigenvalue problem singular.

t14=Append[Table[t13[[i]],{i,1,(n/2)}],b1]

Because eigenvalue problem is singular and we need to algebraically remove the last equation.  We will solve for the highest coefficient, a[n], and then use result this to eliminate a[n] everywhere in the array.

Here we just take the last equation (by taking the last row, multiplying by the a[i]'s) and then solve for a[n].

s1=Solve[t14[[(n/2)+1]].Table[a[i],{i,0,n,2}]==0,a[n]]

Let's be inelegant.  Remake the set of equations and then do the elimination of a[n] in the rest of the system.  (You can redo this another way!)

t15=Table[t14[[i]].Table[a[i],{i,0,n,2}],{i,1,(n/2)}]

Here is the substitution for a[n] which runs off of the page.

t16=Expand[t15/.s1]

Now remake the eigenvalue matrix, which is now of lower order, but which has all finite eigenvalues.

t17=Transpose[Table[Table[Coefficient[t16[[1,i]],a[j]],{i,1,(n/2)}],{j,0,n-2,2}]]

Here we compute the eigenvalues.  Note the insertion of the - sign.   

We expect them to decrease in accuracy.  Thus if we want n eigenvalues, we might need n+4 or more modes in the original expansion.

Mathematica sorts these for us and does not seem to complain about computing them.  I would not try to get them algebraically however!!

neigs= N[Eigenvalues[-t17/.λ->0]]

At this point, for many problems, you might be finished.  You have computed the eigenvalues.

Here we compute the eigenfunctions for each eigenvalue

If you need the eigenfunctions, that correspond to each of the eigenvalues, you will have to solve the matrix, t17, which enforces the boundary condition, for the a[i]'s.  Of course, there cannot be a unique solution because this is an eigenvalue problem, matrix t17 is (it better be) singular.

(Note that the eigenfunctions are probably less accurate than the eigenvalues!!  Thus you will need even more modes to get the eigenfunctions correct.)

 Here we again remake the equations

t18=Table[t17[[i]].Table[a[i],{i,0,n-2,2}],{i,1,(n/2)-1}]

Set them up to solve.

t19=LogicalExpand[t18==0]

Solve them.

t20=Solve[t19,Table[a[i],{i,2,n-2,2}]]

Now we can compute the eigenfunctions for any of the eigenvalues.  Here is the basic expansion.  Note that I have left off the last term.

expand2=(∑_({{i=0}, {Δ i=2}})^(n-2))a[i] ChebyshevT[i,y]

To get the eigenfunction, you just substitute the a[i]'s and the desired λ, like this.  Remember that the solution is not unique so you have an arbitrary scale factor, a[0], that you can pick.  For numerical problems you usually need to do some sort of normalization.  We can pick the a[0]'s so that each eigenfunction has a value of 1 at y=0.

efunc1=((expand2/.t20)/.λ->neigs[[1]])/.a[0]->1

Here are some plots of eigenfunctions from previous calculations.

neigs

λ = 1.4674

efunc1=(expand2/.t20)/.λ->neigs[[1]]

Here is where we normalize the eigenfunction to have a value of 1 at y=0.

a01=a[0]/. Solve[((efunc1/.y->0)[[1]])==1,a[0]][[1]]

eplot1=Plot[efunc1/.a[0]->a01,{y,-1,1},
PlotStyle->Dashing[{.02,.02}]]

Here we compare the first  analytical and numerical eigenfunctions

Show[aplot1,eplot1]

λ = 21.2066

efunc2=(expand2/.t20)/.λ->neigs[[2]]

Here is where we normalize the eigenfunction to have a value of 1 at y=0.

a02=a[0]/.Solve[((efunc2/.y->0)[[1]])==1,
a[0]][[1]]

eplot2=Plot[efunc2/.a[0]->a02,{y,-1,1},
PlotStyle->Dashing[{.02,.02}]]

Here we compare the second  analytical and numerical eigenfunctions

Show[aplot2,eplot2]

λ = 60.7158

efunc3=(expand2/.t20)/.λ->neigs[[3]]

Here is where we normalize the eigenfunction to have a value of 1 at y=0.

a03=a[0]/.Solve[((efunc3/.y->0)[[1]])==1,
a[0]][[1]]

eplot3=Plot[efunc3/.a[0]->a03,{y,-1,1},
PlotStyle->Dashing[{.02,.02}]]

Here we compare the first  analytical and numerical eigenfunctions

Show[aplot3,eplot3]

For n=12, this one is looking a little off because the λ is inaccurate.  It should be

N[(7 Pi/2)^2-1]

However if we use the computed value:

λ = 126.83

efunc4=(expand2/.t20)/.λ->neigs[[4]]

Here is where we normalize the eigenfunction to have a value of 1 at y=0.

a04=a[0]/.Solve[((efunc4/.y->0)[[1]])==1,
a[0]][[1]]

This one does not match well to the analytical solution.

eplot4=Plot[efunc4/.a[0]->a04,{y,-1,1},
PlotStyle->Dashing[{.02,.02}]]

We see that at this point, the error is quite noticeable.  This is the fourth eigenfunction.  Compare this to the number of modes that were used, = n/2+1.

Show[aplot4,eplot4]

Here are the numerical ones all together.

Show[eplot1,eplot2,eplot3,eplot4]

Here are the analytical eigenfunctions.

Show[aplot1,aplot2,aplot3,aplot4]

We see that this does not match real well with the same analytical modes as the mode frequency increases.  The numerical ones vary in amplitude across the interval!!.

Now you can get the eigenvalues and the eigen functions for a differential eigenvalue problem using a Chebyshev spectral method!

Orr-Sommerfeld Equation

This try with u[x,y,t] = f[x,t] psi[y]

The linearized x-momentum equations for a nearly-parallel flow are
x direction (first in traditional form)

xmom=((∂p(x,y,t))/(∂x))-(((((∂^2)u(x,y,t))/(∂(x^2)))+(((∂^2)u(x,y,t))/(∂(y^2))))/Re)+((∂u(x,y,t))/(∂t))+((∂u(x,y,t))/(∂x)) u0(y)+((∂u0(y))/(∂y)) v(x,y,t)

xmom=(∂_x)p[x,y,t]-(((∂_({x,2}))u[x,y,t]+(∂_({y,2}))u[x,y,t])/Re)+(∂_t)u[x,y,t]+(∂_x)u[x,y,t] u0[y]+(∂_y)u0[y] v[x,y,t]

Back to conclusions

y direction (first in traditional form)

ymom=((∂p(x,y,t))/(∂y))-(((((∂^2)v(x,y,t))/(∂(x^2)))+(((∂^2)v(x,y,t))/(∂(y^2))))/Re)+((∂v(x,y,t))/(∂t))+((∂v(x,y,t))/(∂x)) u0(y)

ymom=(∂_y)p[x,y,t]-(((∂_({x,2}))v[x,y,t]+(∂_({y,2}))v[x,y,t])/Re)+(∂_t)v[x,y,t]+(∂_x)v[x,y,t] u0[y]

The continuity equation is

cont=((∂u(x,y,t))/(∂x))+((∂v(x,y,t))/(∂y))

cont=(∂_x)u[x,y,t]+(∂_y)v[x,y,t]

xos1 = D[xmom,y] - D[ymom,x]

We need to expand these in terms of normal modes where the variable will be assumed to have the form ξ[x,y,t]=(OverHat[ξ])[y] Exp[I(α x -α c t)] based on the idea that we need functions that can describe spatially and temporally varying disturbances

Back to conclusions

xos2=xos1/.{u[x,y,t]->uh[y] ff[x,t] ,v[x,y,t]->vh[y]  ff[x,t],Derivative[a1_,a2_,a3_][u][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(uh[y]  ff[x,t]),Derivative[a1_,a2_,a3_][v][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(vh[y]  ff[x,t])}

which is the result for the y equation

os5=cont/.{u[x,y,t]->uh[y] ff[x,t] ,v[x,y,t]->vh[y]  ff[x,t],Derivative[a1_,a2_,a3_][u][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(uh[y]  ff[x,t]),Derivative[a1_,a2_,a3_][v][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(vh[y]  ff[x,t])}

which is the continuity equation

Now we use the continuity equation and the disturbance stream function definition
which is uh[y]=((∂ψ[y])/(∂y)) and vh[y]=-I αψ[y].

Back to conclusions

xos3=xos2/.{uh[y]->(∂_y)ψ[y],vh[y]->-D[ff[x,t],x] ψ[y]/ff[x,t],Derivative[a1_][uh][y]:>(∂_({y,a1}))((∂_y)ψ[y]),Derivative[a1_][vh][y]:>(∂_({y,a1}))(-D[ff[x,t],x] ψ[y]/ff[x,t])}

xos4 =xos3/.{ff[x,t]-> Exp[I (α x-α c t)],Derivative[a1_,a2_][ff][x,t]:>(∂_({x,a1},{t,a2}))( Exp[I (α x-α c t)])}

xos5= ExpandAll[xos4/Exp[I (α x-α c t)]]

Cancel[Expand[(xos5/(I α))]]

(1/(Re α))(ψ(y) (I (α^4)+c Re (α^3)-Re u0(y) (α^3)-Re u0''(y) α)+α (-c Re+u0(y) Re-2 I α) ψ''(y)+I Derivative[4][ψ](y))

Expand[%]

This, with little effort, is the Orr-Sommerfeld equation.

This try with u[x,y,t] = f[x] g[t] psi[y]

The linearized x-momentum equations for a nearly-parallel flow are
x direction (first in traditional form)

xmom=((∂p(x,y,t))/(∂x))-(((((∂^2)u(x,y,t))/(∂(x^2)))+(((∂^2)u(x,y,t))/(∂(y^2))))/Re)+((∂u(x,y,t))/(∂t))+((∂u(x,y,t))/(∂x)) u0(y)+((∂u0(y))/(∂y)) v(x,y,t)

xmom=(∂_x)p[x,y,t]-(((∂_({x,2}))u[x,y,t]+(∂_({y,2}))u[x,y,t])/Re)+(∂_t)u[x,y,t]+(∂_x)u[x,y,t] u0[y]+(∂_y)u0[y] v[x,y,t]

Back to conclusions

y direction (first in traditional form)

ymom=((∂p(x,y,t))/(∂y))-(((((∂^2)v(x,y,t))/(∂(x^2)))+(((∂^2)v(x,y,t))/(∂(y^2))))/Re)+((∂v(x,y,t))/(∂t))+((∂v(x,y,t))/(∂x)) u0(y)

ymom=(∂_y)p[x,y,t]-(((∂_({x,2}))v[x,y,t]+(∂_({y,2}))v[x,y,t])/Re)+(∂_t)v[x,y,t]+(∂_x)v[x,y,t] u0[y]

The continuity equation is

cont=((∂u(x,y,t))/(∂x))+((∂v(x,y,t))/(∂y))

cont=(∂_x)u[x,y,t]+(∂_y)v[x,y,t]

xos1 = D[xmom,y] - D[ymom,x]

We need to expand these in terms of normal modes where the variable will be assumed to have the form ξ[x,y,t]=(OverHat[ξ])[y] Exp[I(α x -α c t)] based on the idea that we need functions that can describe spatially and temporally varying disturbances

Back to conclusions

xos2=xos1/.{u[x,y,t]->uh[y] ff[x] g[t] ,v[x,y,t]->vh[y] ff[x] g[t],Derivative[a1_,a2_,a3_][u][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(uh[y]  ff[x] g[t]),Derivative[a1_,a2_,a3_][v][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(vh[y]  ff[x] g[t])}

which is the result for the y equation

os5=cont/.{u[x,y,t]->uh[y] ff[x] g[t] ,v[x,y,t]->vh[y] ff[x] g[t],Derivative[a1_,a2_,a3_][u][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(uh[y]  ff[x] g[t]),Derivative[a1_,a2_,a3_][v][x,y,t]:>(∂_({x,a1},{y,a2},{t,a3}))(vh[y]  ff[x] g[t])}

which is the continuity equation

Now we use the continuity equation and the disturbance stream function definition
which is uh[y]=((∂ψ[y])/(∂y)) and vh[y]=-I αψ[y].

Back to conclusions

xos3=xos2/.{uh[y]->(∂_y)ψ[y],vh[y]->-ff'[x] ψ[y]/ff[x],Derivative[a1_][uh][y]:>(∂_({y,a1}))((∂_y)ψ[y]),Derivative[a1_][vh][y]:>(∂_({y,a1}))(-ff'[x] ψ[y]/ff[x])}

xos4 =xos3/.{ff[x]-> Exp[I (α x)],Derivative[a1_][ff][x]:>(∂_({x,a1}))( Exp[I (α x)])}

xos5= ExpandAll[xos4/Exp[I (α x)]]

Cancel[Expand[(xos5/(I α g[t]))]]

Solve[%476==0,g'[t]]

g'[t]/.%[[1]]

Expand[%]

NDSolve[{

This, with little effort, is the Orr-Sommerfeld equation.

Let’s do a spatial spectral expansion and then integrate in time

xos5

This section expands the terms of the Orr-Sommerfeld equation in terms of the Chebyshev polynomials and gets the coefficients using the orthogonality.

We need use only even modes because the boundary conditions are symmetric.  There are then two boundary conditions, no slip and no flow through the wall, thus we will drop off the two highest equations produced from substitution of the Chebyshev's in the ODE.

Here is the start of the calculations.  
In this example we use nn/2 +1 modes.

nn=12;

nx=(nn/2)+1;

Here is the rest of the Orr-Sommerfeld equation

XD=xos5

We can substitute the average profile and the (Overscript[U, ─])''(y)

xd1=XD/.{u0[y]->1-(y^2),u0''[y]->(∂_({y,2}))(1-(y^2))}

xd2 = xd1 /. 
   {ψ[y] -> a[i]*ChebyshevT[i, y], 
    D[ψ[y], {y, a1_}] -> D[a[i]*ChebyshevT[i, y], {y, a1}]}

tempxd=(∑_({{i=0}, {Δ i=2}})^nn)xd2;

tempxda=Simplify[tempxd]

z0d=Table[a0d[i]=((2 (∫_(-1)^1)tempxda ChebyshevT[i,y] (1/(Sqrt[1-(y^2)]))\[DifferentialD]y)/(π c[i])),{i,0,nn,2}]

Here are the boundary conditions.

vexpand=(∑_({{i=0}, {Δ i=2}})^nn)a[i] ChebyshevT[i,y]

bc1=vexpand/.{y->-1}

bc2=(∂_y)vexpand/.{y->-1}

Tau solution method

We solve for the two highest coefficients because the two highest modes are dropped from the matrix of the differential equation mode coefficients (the Tau method).  The boundary conditions replace these modes assure that the BC's are satisfied.  However, the BC equations do not contain the eigen value and thus must be eliminated algebraically before solving the eigenvalue problem.

bcsolve=Solve[{bc1==0,bc2==0},{a[nn],a[nn-2]}]

Here is the table of coefficient equations from the O-S equation

dcoefs=z0d

Here is where we use the boundary conditions to replace the two highest coefficients.

dcoefsub=dcoefs/.bcsolve[[1]]

Here we pick out the number of valid equations

coefs=Table[dcoefsub[[i]],{i,1,nx-2}]

temp1=Expand[(coefs/.c[0]->2)/.Table[c[i]->1,{i,2,nn,2}]]

%[[1]]

Solve[%==0,g'[t]]

Here we construct the matrix that comprises the generalized eigenvalue problem 
AB = (AA  -   B  λ)

AB=Table[Table[Coefficient[temp1[[j]],a[i]],{i,0,nn-4,2}],{j,1,nx-2}]

The left side matrix will be AA

AA=AB/.g'[t]->0

Now we get BB

B=((AB-AA)/(g'[t]))

At this point we start doing things numerically.
First choose α and Reynolds number

alfx=1.02056;
rrx=5772.22;

Here is the calculation of Binverse to solve the GEP.

Binv=Inverse[-B];

Here is the matrix of the now regular eigenvalue problem

Binv.AA = AZ

AZ=Binv.AA;

Now make it numerical for sure

N[AZ/.{α->alfx,rr->rrx}]

Get the eigenvalues.

Eigenvalues[N[AZ/.{α->alfx,rr->rrx}]]

These are not very accurate!!  (See below for the good ones)

Try a different Reynolds number and wavenumber.

alfx=1;
rrx=10000;

Here is the calculation of Binverse to solve the GEP.

Binv=Inverse[-B];

Here is the matrix of the now regular eigenvalue problem

Binv.AA = AZ

AZ=Binv.AA;

Now make it numerical for sure

N[AZ/.{α->alfx,rr->rrx}]

Get the eigenvalues.

Eigenvalues[N[AZ/.{α->alfx,rr->rrx}]]