Linear Algebra Powertool

Polynomial Approximation - III

Jacobi Polynomials

Worksheet by Mike May, S.J., maymk@slu.edu

> restart;

Overview

We have been looking at polynomial approximations to continuous functions defined on the domain [-1, 1] using orthogonal polynomials of various kinds. In this worksheet we use yet another norm, produing the Jacobi polynomials, paying attention to how this changes the definition of best fit and the polynomials that can be approximated.

We begin by loading Maple's package for othogonal polynomials. We also load the plots package so that we can graph our results. We repeat the technical commands from the last worksheet.

> with(orthopoly); with(plots):
assume('i',integer): assume('j',integer):
interface(showassumed=0):

Warning, the name changecoords has been redefined

For approximation by orthogonal polynomials, we define the norm by

<f(x), g(x)> = .,

where wtfunc(x) is a weight function. For Legandre polynomials, the weight function is the constant 1, and all points in the domain have equal weight in defining best fit approximation. For Chebyshev polynomials of the first kind, the weight function is (1-x^2)^(-1/2), so extra weight is added to the ends of the domain.

With Jacobi polynomials, with parameters a and b, the weight function is (1-x)^a*(1+x)^b, where a and b are both > -1. When a = b = 0 this reduces to the case of Legandre polynomials. If a and b are both positive the weight of the ends of the intervals is decreased. If a and b are negative, the weight of the ends of the intervals is increased.

Consider the plot of the weight functions for Jacobi polynomials with a variety of values for the parameters a and b.

> plot({1, (1-x)*(1+x), (1-x)^2*(1+x)^3, (1-x)^6*(1+x)^8,
(1-x)^(-.4)*(1+x)^(-.6), (1-x)^(-.9)*(1+x)^(-.9)}, x=-1..1, y=0..3);

We would also like to give a numeric measure to how the weighting shifts with different parameter values. One easy intuitive way is to compute the percent of the weight that is in the middle 60 % of the domain.

> "w(0, 0)" = int(1,x=-0.6..0.6)/int(1,x=-1..1);
"w(1, 1)" = int((1-x)*(1+x),x=-0.6..0.6)/int((1-x)*(1+x),x=-1..1);
"w(2, 3)" = int((1-x)^2*(1+x)^3,x=-0.6..0.6)/int((1-x)^2*(1+x)^3,x=-1..1);
"w(6, 8)" = int((1-x)^6*(1+x)^8,x=-0.6..0.6)/int((1-x)^6*(1+x)^8,x=-1..1);
"w(-.4, -.6)" = evalf(Int((1-x)^(-0.4)*(1+x)^(-0.6),x=-0.6..0.6))/
evalf(Int((1-x)^(-0.4)*(1+x)^(-0.6),x=-1..1));
"w(-.9, -.9)" = evalf(Int((1-x)^(-0.9)*(1+x)^(-0.9),x=-0.6..0.6))/
evalf(Int((1-x)^(-0.9)*(1+x)^(-0.9),x=-1..1));

>

As mentioned above, one reason to change the parameter values is to change the way closeness is defined. A second reason is to change the subspace of functions over which the inner product makes sense. In particular,using positive values for a and b produces an inner product under which functions with poles at the end of the domain have finite norm. Consider what happens to the function f(x) = 1/(1-x) with the Legandre norm and with the (3,3) Jacobi norm.

> "<f,f>_L" = int(1*(1-x)^(-2), x=-1..1);
"<f,f>_J(3,3)" = int((1-x^2)^3*(1-x)^(-2), x=-1..1);

>

Exe rcises:

1) Describe the rational functions that have a finite norm under the (a, b) Jacobi norm. (To simplify the problem, assume that the rational functions are in a normal form with partial fraction decomposition, i.e., each function is a sum of proper fractions with the denominator a power of an irreducible polynomial with integer coefficients.)

>

2) For the (2,3) Jacobi norm, what percentage of the weight is on the interval [0, 1]? Repeat the question with the (2, 5) Jacobi norm.

>

Jacobi polynomials

As we did in previous worksheets we check the help page for the orthogonal polynomials we are looking at in this worksheet.

> ?orthopoly,P

Once we have loaded the orthopoly package, we obtain the nth (a,b) Jacobi polynomial with P(n,a,b,x);

> "P(3,1,1,x)"=P(3,1,1,x);
"P(3,1,2,x)"=P(3,1,2,x);
"P(3,2,1,x)"=P(3,2,1,x);
"P(3,-.4,-.6,x)"=P(3,-.4,-.6,x);

E xercise:

3) Find the 2nd and 5th Jacobi polynomials with parameters (1,1), (2,3), and (6,6). Verify that pair associated with each parameter is orthogonal under the oppropriate inner product.

>

An extended example

Time to look at an example to see how this works out in a particular case. Once again we will use sin(3¹x) as our function.

> func := x -> sin(3*Pi*x);

Computing coefficients and approximating polynomials

Next we start computing coefficients and defining the approximating polynomials up to degree 12 for the inner products of type (0,0), (2,3), and (-.4, -.6). Recall that if P(n,a,b,x) is the nth (a,b) Jacobi polynomial, the projection coeffient of f(x) onto P(n,a,b,x) should be divided by .

> for i from 0 to 12 do
normP00[i] := evalf(Int((P(i,0,0,x))^2, x=-1..1)):
coeffP00[i] := evalf(Int(simplify(P(i,0,0,x)*func(x)*1.0), x=-1..1));
od:
Jacobi00approx := proc(m,x)
local count:
sum(P(count,0,0,x)*coeffP00[count]/normP00[count],count=0..m):
end;
for i from 0 to 12 do
normP23[i] := evalf(Int((1-x)^2*(1+x)^3*(P(i,2,3,x))^2, x=-1..1)):
coeffP23[i] := evalf(Int(simplify((1-x)^2*(1+x)^3*P(i,2,3,x)*func(x)*1.0), x=-1..1)):
od:
Jacobi23approx := proc(m,x)
local count:
sum(P(count,2,3,x)*coeffP23[count]/normP23[count],count=0..m):
end;
for i from 0 to 12 do
normP46[i] := evalf(Int((1-x)^(-.4)*(1+x)^(-.6)*(P(i,-.4,-.6,x))^2, x=-1..1));
coeffP46[i] := evalf(Int(simplify((1-x)^(-.4)*(1+x)^(-.6)*
P(i,-.4,-.6,x)*func(x)*1.0), x=-1..1));
od:
Jacobi46approx := proc(m,x)
local count:
sum(P(count,-.4,-.6,x)*coeffP46[count]/normP46[count],count=0..m):
end;

Comparing the approximations

Our past experience with this function is that the approximation start to be interesting with the ninth degree. Note first that the polynomial approximations are significantly different.

> n:= 9:
sort(Jacobi00approx(n,x));
sort(Jacobi23approx(n,x));
sort(Jacobi46approx(n,x));

Next compare the graphs of the approximations and of the errors.

> lowy := -1.2: highy := 1.6: n:= 7:
pl1 := plot(func(x), x=-1..1,y=lowy..highy, color=red):
pl2 := plot(Jacobi00approx(n,x), x=-1..1,y=lowy..highy, color=green):
pl3 := plot(Jacobi23approx(n,x), x=-1..1,y=lowy..highy, color=blue):
pl4 := plot(Jacobi46approx(n,x), x=-1..1,y=lowy..highy, color=black):
B := textplot({[-0.8,highy-(highy-lowy)/20,
`The function and the `||n||` term Jacobi approximations`],
[-0.8,highy-(highy-lowy)/9,
`(0,0) (green) vs (2,3) (blue), vs (-.4, -.6) (black)`]},
align=RIGHT, font = [TIMES, BOLD, 12] ):
display({pl1, pl2, pl3, pl4, B},view=[-1..1, lowy..highy]);

As expected, the (2, 3) approximation tends to be best in the center of the graph, while the (-.4, -.6) approximation tends to be best on the ends. This becomes even clearer when we look at the graphs of the errors.

> lowy := -1.0: highy := 1.4: n:= 10:
pl1 := plot(func(x) - Jacobi00approx(n,x), x=-1..1,y=lowy..highy, color=green):
pl2 := plot(func(x) - Jacobi23approx(n,x), x=-1..1,y=lowy..highy, color=blue):
pl3 := plot(func(x) - Jacobi46approx(n,x), x=-1..1,y=lowy..highy, color=black):
B := textplot({[-0.8,highy-(highy-lowy)/20,
`The `||n||` term Jacobi error `],
[-0.8,highy-(highy-lowy)/9,
`(0,0) (green) vs (2,3) (blue), vs (-.4, -.6) (black)`]},
align=RIGHT, font = [HELVETICA, BOLD, 12] ):
display({pl1, pl2, pl3, B});

It is clear that the (2,3) Jacobi approximation is best in the middle and the (-.4, -.6) Jacobi approximation is best on the ends.

>

Exercises:

>

4) Describe the kind kind of error you would expect with a (-.6, 5) Jacobi approximation.

>

5) Plot the function y=cos(2*Pi*x) along with the degree 10 (-.6, 5) Jacobi approximation. In a separate graph plot the error of the approxiamtion.

>

>

>