Solution of ODE's and eigenvalue problems with a Chebyshev polynomial 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, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version:  5/18/99
More recent versions of this notebook should be available at the web site:
http://www.nd.edu/~mjm/spectral.ode.eigens.nb

Summary

This notebook demonstrates the Orszag-tau ( a modification of the Galerkin) spectral method for a simple ODE and an ODE eigenvalue problem.  It is intended as a first introduction to solving these problems with a spectral numerical 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 properties of the Chebyshev polynomials using orthogonality.

Reference:  R. Miesen and B. J. Boersma (1995) "Hydrodynamic stability of a sheared liquid film", Journal of Fluid Mechanics, 301 pp 175-202.

Mathematica aside

In Mathematica, it is convenient to give all expressions a "name".  I try to pick ones that are consistent with what is being done (but sometimes "temp" is used).  This assignment is done with an "=" sign.  To make an equation, a "==" is used.  This distinction is very useful in computer algebra and is employed in all of the packages with which I am familiar .

Input notation

I have used a mix of input notation to show how it can be done.

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

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

In Standard form, which is a little shorter, you would need the typeset window to make this practical.

\(dc = 
      \(-\((γ\ ∂\_{x, 2}η[x, t])\)\) - 
        ρ\_1\ 
          \((\(-∂\_t ϕ\_1[x, y, t]\) - 
              U1\ ∂\_x ϕ\_1[x, y, t] - g\ η[x, t])\) + 
        ρ\_2\ 
          \((\(-∂\_t ϕ\_2[x, y, t]\) - 
              U2\ ∂\_x ϕ\_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 endpoint values @y=1 of 1, and @ y=-1, of 1 or -1 depending if they have even or odd index numbers.
  
Let's see what these are

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

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 expansion functions.

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 the 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.

Here are some more examples.

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

Here are 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.

You probably know from that it is possible to fit an arbitrary function using a set of orthogonal polynomials by forming the inner product to evaluate the coefficients.  (If you need a basic introduction to spectral theory you my want to check out http://www.nd.edu/~mjm/spectral.solutions.nb .)  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]  =\(TraditionalForm
      \`∑\+\(j = 1\)\%∞ c\_j\ \(\(ϕ\_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.

Return to conclusions

Solution of a simple ode

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

Analytical solution

Here is a simple differential equation that we can get an analytic solution to. 
 It will be useful for comparison when we solve it with a spectral method.
 
 Note it is first written in "traditional form" and the cell is inactive.

ode=(((∂^2) u(y))/(∂(y^2))) + u(y)

Now we solve the ODE.

\(ode = ∂\_{y, 2}u[y] + u[y]\)

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

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

p1 = Plot[soln1, {y, -1, 1}];

Numerical solution

Now solve the equation with a Galerkin 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 expressions.
 
 Note it is first written in "traditional form" and the cell is inactive.

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

Here we start the numerical solution.

\(ode1 = ∂\_{y, 2}vv[y] + 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.  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 = 12;

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.

This sum, with appropriate values for the a(i)'s would have a value of approximately 0 and thus solve the ODE.

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

Return to conclusions

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

Return to conclusions

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

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, {i, 0, n, 2}];

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

temp7=(temp6/.c[0]->2)/.Table[c[i]->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^(2  n))/ √(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^(2  n))/ √(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]]

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 n/2 unknown a(j)'s and n/2 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. 

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

Return to conclusions

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]

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]

Return to conclusions

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 is a 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]

Here are the two methods compared, first the Spectral approach:

testtable = {{2, 0.225926560420915478`}, {4, 
          0.000951142525143355577`}, {6, -6.18435756960661642`*^-7}, {8, 
          8.9471794231341768`*^-10}, {10, -9.81661958634065534`*^-13}, {
          12, 5.91919878400052912`*^-16}, {14, 9.66027008208424398`*^-17}, {
          16, 1.04611309920646089`*^-17}}

Here is the power series approach,

Return to conclusions

testpowertable = {{2, -0.218517884023529118`}, {4, 
          0.00712314161749632557`}, {6, -0.00013130059838550192`}, {8, 
          1.4845532182493013`*^-6}, {10, -1.13778653577798039`*^-8}, {12, 
          6.30329954060986352`*^-11}}

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}];

Here are the first four.

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

Computation of the eigenvalues

Here is our equation again.

\(odeeig1 = ∂\_{y, 2}vv[y] + \((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}]]

Here are the factors that we pulled out.

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

Now we just need to integrate each of them once.  If we had left them interspersed in the equations, we would have done about (((2  n))^3) integrations!

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, 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}]]

Return to conclusions

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]]

Here are the analytically determined eigenvalues from above.  Note how the accuracy decreases as they get bigger.  If you need them very accurately, you will need several more modes than the number of eigenvalues you need.

eigsanlyt

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

Computation of the eigenfunctions

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.

Return to conclusions

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

Here are some plots of eigenfunctions for different values of λ.

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 analytical and numerical eigenfunctions

Show[aplot3,eplot3];

The first three (for n=12) were not too bad but the fourth one will look worse.  Note that the eigenvalue is not real accurate.    The analytical value is:

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

The computed value is λ = 126.83

We will see below it is not just this inaccuracy in  λ that is responsible.

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}]];

efunc4x=(expand2/.t20)/.λ->eigsanlyt[[4]]

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

a04x=a[0]/.Solve[((efunc4x/.y->0)[[1]])==1,
a[0]][[1]]

This one does not match well to the analytical solution.

eplot4x=Plot[efunc4x/.a[0]->a04x,{y,-1,1},
PlotStyle->Dashing[{.02,.02}]];

Here we show the eigen function with the correct λ and the numerical one.  There is not much difference.

Show[eplot4, eplot4x];

Now compare the numerical one to the analytical one.  The error is quite noticeable.  This is the fourth eigenfunction and we generated a total of 6.  This gives some idea of how many modes can be expected to be accurate in one of these calculations.

Show[aplot4,eplot4];

Here are the numerical ones all together.

Show[eplot1,eplot2,eplot3,eplot4];

Return to conclusions

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!

Conclusions

1.  A spectral numerical method assumes that the approximate solution to (in this case) an ODE can be represented as a sum of Orthogonal Polynomials with appropriate coefficients.  To find this solution, you must find these coefficients.  

2.  The best way to generate equations for the unknown coefficients is to take advantage of the orthogonality property of the polynomials.  That is, form the inner product of the substituted expansion and each Poynomial, then integrate on the interval (or use some other appropriate inner product.).    

3.  Once you have a set of algebraic equations, obtained from the ODE for each mode, the boundary conditions are enforced by dropping the necessary number of equations for the highest modes and replacing these with the equations for the boundary conditions.  

4.  An alternative solution strategy would be to just use a power series and equate each power of y to 0.  This works for a simple problem but not as well as the tau-spectral method.  

5.  The eigenvalue problem works much the same way except that you need to solve for the eigenvalues at somepoint.  This can be done without evaluating all of the coefficients.  However, you need to evaluate the coefficients to determine the eigenfunctions.  

6.  This method works well but it is important to note that the accuracy of the eigenvalues and eigenfunctions decreases as the mode number increases.  Thus you will need several extra modes to get the desired number accurately.