With interpolation we were given a formula or data about a function f(x), and we made a model p(x) that passed through a given set of data points. We now consider approximation, in which we are still trying to build a model, but specify other conditions that the model must satisfy. It is still likely that p(x) will be equal to the function f(x) at some points, but we will not know in advance which points.
As with interpolation, we build this model out of a set of basis functions. The model is then a "recipe", the linear coefficients c that specify how much of each basis function to use when building the model.
Polynomials will be our most typical set of basis function. However, we will want to use a different set of basis polynomials than the obvious set of 1, x, x2, and so on. We'll still be building the same kind of objects, just using a more convenient set of tools.
Once we have our basis set, we will consider how we can determine the approximating function p(x), and we will look at the behavior of the approximation error.
It is sometimes necessary to refer to part of a vector or matrix. Matlab provides this capability using a colon. Supposing that the vector v has 100 components and you want to square the 25th through 40th components. You can do this by writing
v(25:40)=v(25:40).^2
Matlab provides the capability of defining functions "in line" instead of writing m-files to do it. This feature is most convenient when the function to be defined is very simple--a line of code, say--or when you have a function that requires several arguments and you want to fix all but one. The way to do this is use the inline function. (The full help description of inline is available here.) Suppose, for example, you want to define a function sq(x)=x^2 . You could do this by writing
sq=inline('x.^2','x');
You would then use this function just as if you had written it as an m-file: y=sq(x), or you could plot this function on the interval from 0 to 2 using fplot(sq,[0,2])
The Legendre polyonomials are a basis for the set of polynomials, appropriate for use on the interval [-1,1]. They are discussed in Atkinson starting on page 210, and the first few Legendre polynomials are:
P0(x) = 1
P1(x) = x
P2(x) = ( 3 x2 - 1 ) / 2
P3(x) = ( 5 x3 - 3 x ) / 2
P4(x) = ( 35 x4 - 30 x2 + 3 ) / 8
The value at x of any Legendre polynomial Pi can be determined using the following recursion:
P0(x) = 1
P1(x) = x
...
Pi(x) = ( (2*i-1) * x * Pi-1(x) - (i-1) * Pi-2(x) ) / i
Exercise 1: write a routine legpoly.m to evaluate the i-th Legendre polynomial Pi(x). Your routine should have the form:
function pval = legpoly ( i, x )
and the routine should work for vector arguments x.
The first thing you should note is that if i is larger than 2, you only need to retain the values Pj-1(x) and Pj-2(x) in order to compute Pj(x). Instead of creating a subscripted variable, just use the three variables pj, pjm1, and pjm2, where the "jm1" stands for "j minus 1" and "jm2" stands for "j minus 2." Each time you increase the value of j you move pjm1 into pjm2 and pj into pjm1 to keep things straight. Your code will look something like the following:
if i is 0 or 1
set pval
else
pjm1 = ones(size(x));
pj = x;
for j=2:i
pjm2=pjm1;
pjm1=pj;
pj= ???
end
pval=pj;
end
The ones function makes a vector the same size as x but containing all ones.
As a test of your routine, verify the values in the following table by evaluating legpoly using the vector x=[-0.5, 0, 1].
P\x| -0.5 0.0 1.0
---+----------------------------
P0 | 1.0000 1.0000 1.0000
P1 | -0.5000 0.0000 1.0000
P2 | -0.1250 -0.5000 1.0000
P3 | 0.4375 0.0000 1.0000
P4 | -0.2891 0.3750 1.0000
P5 | -0.0898 0.0000 1.0000
The Legendre polynomials form a basis for the linear space of polynomials. One thing we like any set of basis vectors to do is be orthogonal. If we were working with regular geometric vectors, we could draw them and see this condition. In generalized vector spaces, we have to define a dot product, and say that two "vectors" f(x) and g(x) are orthogonal if their dot product (f,g) (the integral of f(x) times g(x)) is equal to zero.
How shall we compute the integrals? We will be learning how numerical integration is done in the next lab, but for now we can use the function supplied by Matlab called quad8 ("quad" for "quadrature" which is another word for integration). The quad8 function requires at least three variables, a function name, the lower integration bound, and the upper integration bound. The integral of the sin(x) function from 0 to pi would be given by quad8('sin',0,pi). To find the integral dot product of, say, P2 and P1 you might first use inline to define the intgrand,
integrand=inline('legpoly(2,x).*legpoly(1,x)','x');
(note the use of elementwise multiplication because quad8 uses vector operations). Then the Integral Dot Product (P2(x),P1(x)) could be found as
IDP=quad8(integrand,-1,1)
You could compute the same Integral Dot Product in a single line as
IDP=quad8(inline('legpoly(2,x).*legpoly(1,x)','x'),-1,1)
Exercise 2: Using your Legendre polynomial routine legpoly, and the quad8 integration routine, estimate the following dot products:
( P3(x), P5(x) )
( P4(x), P4(x) )
( P0(x), P8(x) )
What should you have gotten?
In linear algebra, the length of a vector v can be computed as
||v|| = sqrt ( v, v )
In the same way, in our generalized vector space, the "length" of a vector is its L2 norm, which is the square root of its integral dot product with itself.
It will be useful for our Legendre polynomials to have unit length. For the rest of our work, we will use normalized Legendre polynomials. The definition above gives that ( Pi(x), Pi(x) ) = 2 / ( 2 * i + 1 ), so the Legendre polynomials can be normalized by multiplying Pi(x) as computed above by sqrt((2*i+1)/2). Before going on, insert this normalization into your legpoly function.
Exercise 3: Using the normalized Legendre polynomials, estimate the L2-norms of:
P0(x)
P4(x)
P8(x)
If the results are far away from 1, you are doing something wrong.
Theorem 4.3 in Atkinson says that a function f(x) defined on [-1,1] can be approximated by the polynomial:
c0 * P0(x) + c1 * P1(x) + ...+ cn * Pn(x)
The coefficients c can be computed using the integral dot product:
ci = ( f, Pi )
The Corollary to Theorem 4.3 essentially says that this choice of coefficients results in a model p(x) of the function f(x) with the property that p(x) "casts the same shadow" as f(x) does. The dot product can be thought of as measuring the shadow or projection of the vector f onto each of the basis functions. The resulting approximation p(x) will indeed satisfy this equation:
( p, Pi ) = ( f, Pi )
for every basis function i. In other words, if the dot product was our only measurement device, we couldn't tell p(x) apart from f(x). But this means that the error p(x)-f(x) that we make is orthogonal to our space - we can't do any better! Atkinson discusses this further in Section 4.5.
Exercise 4: Compute the first n Legendre coefficients for the function ex with n = 5. Compute the dot products c(i+1) = ( f, Pi ), for i from 0 to n. You could use quad8 for integration as described above. Unfortunately, you cannot easily build a loop to do this task, so compute the coefficients c(i) for i=1,2,3,4,5 individually. Use c(1) to denote c0, c(2) to denote c1 and so forth. (Actually, you can do it in a loop using sprintf, but that is more effort than I wish to put you to.) Your coefficients should be close to these values:
i | 0 1 2 3 4 5
--|-----------------------------------------------------------
c | 1.661985 0.901117 0.226302 0.037660 0.004698 0.000469
These are the values in Atkinson's Table 4.2.
Now we have all the pieces in place. Given a function f defined on [-1,1], we can compute the Legendre coefficients c by computing dot products. Then, to evaluate the approximation at some point x, we simply add up c0 times the zero-th Legendre polynomial, c1 times the second and so on.
Exercise 5: Write a routine legval.m to approximate a function f(x) over [-1,1] given its Legendre coefficients. Your routine should have the form:
function pvec = legval ( c, x )
Recall that c(1) contains the coefficient we've called c0, and c(n+1) contains cn. The point or vector x is where we want to evaluate. Your code might look like this:
pvec = zeros ( size ( x ) )
for i = 0:length(c)
p = ?
pvec = pvec + ?
endTest your evaluation routine by evaluating the approximation to ex from the previous exercise. Plot the true and approximated values over the interval [-1,1]. Once you have computed the coefficients, you might like to see what the approximation looks like, using plotting. You could do this with the commands:
xtest = linspace ( -1, 1, 101 )
ytest = exp ( xtest )
ptest = legval ( c, xtest )
You could also plot the error ptest-ytest and see if it oscillates evenly in the way Atkinson discusses. (Because our coefficients are inaccurate, you may see some bad behavior near the endpoints.)
You should wonder what happens over a larger interval, say [-2,2]. You can evaluate the original and approximating functions over this interval, and plot them. Outside of [-1,1], your approximation has no guaranteed behavior whatsoever!
What we have done only "works" on the interval [-1,1]. Suppose you want to approximate sin(x) from 0 to 2*pi? Everything we've done can be adjusted to work in this case, but you have to be careful. The easiest approach is to essentially consider the function sin(x(t)), where x(-1)=0 and x(+1)=2*pi. Now look at this as a function of t...
Another basis that may be used is the Chebyshev polynomials. These are also defined over [-1,1], but use a different dot product involving a weight function. There are also Laguerre, Jacobi, and Hermite polynomials which are appropriate for other dot products or intervals.
Finally, why couldn't we just work with the basis 1, x, x2? There are two reasons, convenience and accuracy. Because this basis is not orthogonal, it's much more work to solve for the coefficients. And because the functions xn and xn+1 are very similar, it becomes quite difficult to compute the coefficients accurately.
Exercise 6: Compute the Legendre polynomial approximation to the function
y = atan ( 5 * x )
on the interval [-1,1], using the first 6 basis functions in the same manner you used for exp(x) above. Explain the pattern you see in the coefficients. Plot both y and its Legendre polynomial approximation on the same graph.
Back to the MATH2070 page.
Last revised on November 10, 2000.