Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Chebyshev Approximation

Generating Chebyshev Polynomials from a Recursion Relation

From the definition (Subscript[T, n])(x)=cos(n arccos x) we can get the recursion relation:

(Subscript[T, 0])(x) = 1
(Subscript[T, 1])(x)=x
(Subscript[T, n+1])(x)=2x (Subscript[T, n])(x)-(Subscript[T, n-1])(x), n=1,2,3,...

First Kind (Subscript[T, n])

Clear[Cheby]
num=5;
Cheby[0]=1;
Cheby[1]=x;
Cheby[n_]:=Cheby[n]=Collect[2x Cheby[n-1]-Cheby[n-2],x]

TableForm[Table[{n,Cheby[n],Cheby[n]/.x->1},{n,0,num}]]

Compare with the actual Chebyshev polynomials, using Mathematica's built-in definition of the Chebyshev polynomial, below.

TableForm[Table[{n,ChebyshevT[n,x],ChebyshevT[n,1]},{n,0,num}]]

Here is a plot of the first few Chebyshev polynomials:

len=11;
Do[plot[n]=Plot[{ChebyshevT[n,x]},{x,-1,1},PlotRange->{{-1,1},{-1,1}},PlotStyle->{Thick,Red},BaseStyle->20],{n,0,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,0,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

Second Kind (Subscript[U, n])

Although not used in the context of interpolation, there is a Chebyshev Polynomial of the Second Kind (Subscript[U, n])(x)=sin((n+1)arccos x)/sin (arccos x)

Clear[Cheby]
num=5;
Cheby[0]=1;
Cheby[1]=2x;
Cheby[n_]:=Cheby[n]=Collect[2x Cheby[n-1]-Cheby[n-2],x]

TableForm[Table[{n,Cheby[n],Cheby[n]/.x->1},{n,0,num}]]

Compare with the actual Chebyshev polynomials, using Mathematica's built-in definition of the Chebyshev polynomial, below.

TableForm[Table[{n,ChebyshevU[n,x],ChebyshevU[n,1]},{n,0,num}]]

Here is a plot of the first few Chebyshev polynomials:

len=11;
Do[plot[n]=Plot[{ChebyshevU[n,x]},{x,-1,1},PlotRange->{{-1,1},{-3,3}},PlotStyle->{Thick,Blue},BaseStyle->20],{n,0,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,0,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

Roots and Extrema of Chebyshev Polynomials of 1st Kind

There are two very nice formulas for the Chebyshev polynomials of the first and second kind, which provide the starting point for a variety of results.

(Subscript[T, n])(x)=cos(n arccos(x))
(Subscript[U, n])(x)=((sin(n arcsin(x)))/(Sqrt[1-(x^2)]))

These can be used to determine recursion relations (the ones we used above), normalization conditions, relations between the 1st and 2nd Kind Chebyshev polynomials, and a host of other properties. We are interested in the zeros and extremas of the Chebyshev polynomials of the 1st Kind. It is somewhat easy to show that the following properties are true.

(Subscript[(OverBar[x]), k]) = cos(((2k-1)/(2n))π), (Subscript[T, n])((Subscript[(OverBar[x]), k]))=0, k=1,2,...,n (zeros)
(Subscript[(OverBar[x])', k])= cos((k/n)π), (Subscript[T, n])((Subscript[(OverBar[x])', k]))=(((-1))^k), k=0,1,...,n (extrema)

Other relations include

(Subscript[T, n+1])(x)=2x (Subscript[T, n])(x)-(Subscript[T, n-1])(x)
(d/(d x))(Subscript[T, n])(x)=n (Subscript[U, n-1])(x)

Monic Chebyshev Polynomials of 1st Kind (OverTilde[(Subscript[T, n])])(x)

Monic polynomials have coefficient of 1 in front of the highest power of x.

num=5
TableForm[Table[{k,Simplify[ChebyshevT[k,x]/2^(k-1)],ChebyshevT[k,1]/2^(k-1)},{k,1,num}]]

len=11;
Do[plot[n]=Plot[{If[n==0,ChebyshevT[n,x],((ChebyshevT[n,x])/(2^(n-1)))]},{x,-1,1},PlotRange->{{-1,1},{-1,1}},PlotStyle->{Thick,Red},BaseStyle->20],{n,0,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,0,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

The zero and extremas are at the same values as for the regular Chebyshev polynomials, however the value of the extrema is scaled as

(Subscript[(OverBar[x]), k]) = cos(((2k-1)/(2n))π), (Subscript[(OverTilde[T]), n])((Subscript[(OverBar[x]), k]))=0, k=1,2,...,n (zeros)
(Subscript[(OverBar[x])', k])= cos((k/n)π), (Subscript[(OverTilde[T]), n])((Subscript[(OverBar[x])', k]))=((((-1))^k)/(2^(n-1))), k=0,1,...,n (extrema)

Polynomial Interpolation: Equispaced Nodes vs (Subscript[T, n])(x) Nodes on Interval [-1,1]

Equally spaced nodes

Clear[LP,fcn,n,x,t]

The following is valid only for n≥2. This is our familiar Lagrange interpolating polynomial on the interval [-1,1].

LP[fcn_,n_,t_]:=Module[{a,b,h,xp,i,k,res},
a=-1.0;
b=1.0;
h=(b-a)/(n-1);
Do[xp[i]=a+h*(i-1),{i,1,n}];
res=Sum[(fcn/.x->xp[k])*((∏_(i=1)^(k-1))((t-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((t-xp[i])/(xp[k]-xp[i]))),{k,1,n}];
Expand[res]
]

Here is a function we wish to approximate by using the Lagrange interpolating polynomial:

f[x_]=1/(1+15x^2)
LP[f[x],4,t]

Here are the Lagrange interpolating polynomial approximations to the function:

len=13;
Do[plot[n]=Plot[{LP[f[x],n,t],f[t]},{t,-1,1},PlotRange->{{-1,1},{0,1}},PlotStyle->{{Thick,Red},{Thick,Blue}},Ticks->None],{n,2,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,2,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

Here is the error of the approximations:

len=13;
Do[plot[n]=Plot[{Abs[LP[f[x],n,t]-f[t]]},{t,-1,1},PlotRange->{{-1,1},{0,1}},PlotStyle->{Thick,Blue}],{n,2,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,2,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

The approximation is not getting better at the edges. This is the wiggle problem we saw earlier.

Chebyshev roots spaced nodes

Clear[LP,fcn,n,x,t]

The following is valid only for n≥2. The change from equally spaced roots is the red bit of code. We know that the Chebyshev polynomials of the 1st kind put more weight on the edge of the region, so we can use that to our advantage.

CP[fcn_,n_,t_]:=Module[{h,xp,i,k,res},
Do[xp[i]=Cos[(2i-1.0)/(2n)*Pi],{i,1,n}];
res=Sum[(fcn/.x->xp[k])*((∏_(i=1)^(k-1))((t-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((t-xp[i])/(xp[k]-xp[i]))),{k,1,n}];
Expand[res]
]

Here is the function:

f[x_]=1/(1+15x^2)
CP[f[x],4,t]

len=13;
Do[plot[n]=Plot[{CP[f[x],n,t],f[t]},{t,-1,1},PlotRange->{{-1,1},{0,1}},PlotStyle->{{Red,Thick},{Blue,Thick}},Ticks->None],{n,2,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,2,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

Notice how the error is not increasing at the edge of the region:

len=13;
Do[plot[n]=Plot[{Abs[CP[f[x],n,t]-f[t]]},{t,-1,1},PlotRange->{{-1,1},{0,1}},PlotStyle->{Blue,Thick}],{n,2,len}]
Table[{plot[n],plot[n+1],plot[n+2]},{n,2,len,3}];
Show[GraphicsGrid[%],ImageSize->{650,Automatic}]

By adding more nodes to the edges, where the problem with the equally spaced nodes exists, we are able to get a better approximation.

Aside: More About Chebyshev Polynomials

Chebyshev polynomials belong to a class of functions called special functions, which inlcude Legendre polynomials, Hermite polynomials, Bessel functions. These special functions have a multitute of properties, and I will just list some of them here.

Orthogonal and Normalization and Weights: Chebyshev Polynomial of 1st and 2nd Kind

Weights

Here are the two different weight functions, which result in Chebyshev polynomial of First or Second Kind.

w1[x_]=(1/(Sqrt[1-(x^2)]));
w2[x_]=(Sqrt[1-(x^2)]);
Plot[{w1[x],w2[x]},{x,-1,1},PlotStyle->{{Red,Thick},{Blue,Thick}},PlotRange->{{-1,1},{0,2}},Frame->True,FrameLabel->{"x","Weight"},ImageSize->{500,Automatic},BaseStyle->20]

w1[x] gives more weight to the edge of the region, whereas w2[x] gives more weight to the center.

(Subscript[T, n])(x)has weight w(x)=(1/(Sqrt[1-(x^2)]))on interval [-1,1] (Chebyshev Polynomial of First Kind)
(Subscript[U, n])(x)has weight w(x)=(Sqrt[1-(x^2)])on interval [-1,1] (Chebyshev Polynomial of Second Kind)

Orthogonal and Normalization

Chebeshev polynomials are orthogonal, which means that on an interval [a,b] relative to a weight function w(x) the integral

(∫_a^b)w(x)(Subscript[ϕ, j])(x)(Subscript[ϕ, k])(x)dx=({{{, ({{0, j≠k}, {(Subscript[α, k]), j=k}})}})

Let's check what the normalization is, using the built-in Mathematica commands for the Chebyshev functions.

w[x]=(1/(Sqrt[1-(x^2)]));
Table[Integrate[ChebyshevT[k1,x]*ChebyshevT[k2,x]*w[x],{x,-1,1}],{k1,0,5},{k2,0,5}]//TableForm

w[x]=(Sqrt[1-(x^2)]);
Table[Integrate[ChebyshevU[k1,x]*ChebyshevU[k2,x]*w[x],{x,-1,1}],{k1,0,5},{k2,0,5}]//TableForm

Chebyshev Polynomials from Gram-Schmidt Process

A Gram Schmidt Module

Here is a way of generating the Chevyshev polynomials, based on the Gram Schmidt process with appropriate weight function and interval.

GramSchmidt[a_,b_,w_,num_]:=Module[{ϕ,B,A,k},

(Subscript[ϕ, 0])[x_]=1;
(Subscript[B, k_]):=(Subscript[B, k])=(((∫_a^b)x*w*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x)/((∫_a^b)w*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x)); 

(Subscript[A, k_]):=(Subscript[A, k])=(((∫_a^b)x*w*(Subscript[ϕ, k-1])[x]*(Subscript[ϕ, k-2])[x]\[DifferentialD]x)/((∫_a^b)w*(((Subscript[ϕ, k-2])[x])^2)\[DifferentialD]x));

(Subscript[ϕ, 1])[x_]=x-(Subscript[B, 1]);

Do[(Subscript[ϕ, k])[x_]=Expand[(x-(Subscript[B, k]))(Subscript[ϕ, k-1])[x]-(Subscript[A, k])(Subscript[ϕ, k-2])[x]],{k,2,num}];

Return[TableForm[Table[{n,(Subscript[ϕ, n])[x],(Subscript[ϕ, n])[1]},{n,0,num}]]]
]

First Kind: weight function is (1/(Sqrt[1-(x^2)]))and interval is x ∈(-1,1)

GramSchmidt[-1,1,(1/(Sqrt[1-(x^2)])),8]

TableForm[Table[{k,ChebyshevT[k,x],ChebyshevT[k,1]},{k,0,8}]]

We are missing the standardization (Subscript[T, n])(x)=1.

Second Kind: weight function is (Sqrt[1-(x^2)])and interval is x ∈(-1,1)

The only change from the First Kind is the weight function w(x).

GramSchmidt[-1,1,(Sqrt[1-(x^2)]),8]

TableForm[Table[{k,ChebyshevU[k,x],ChebyshevU[k,1]},{k,0,8}]]

We are missing the standardization (Subscript[U, n])(x)=1.