|Legendre Polynomials||Exercise 3|
|Orthogonality and Integration||Exercise 4|
|Least squares approximations in||Exercise 5|
|Legendre polynomial approximation||Extra Credit|
|Piecewise constant approximation|
|Piecewise linear approximation (Extra)|
With interpolation we were given a formula or data about a function , and we made a model 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 some condition of ``closeness'' that the model must satisfy. It is still likely that will be equal to the function 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 linear combination with coefficients that specify how much of each basis function to use when building the model.
In this lab we will consider four different selections of basis functions in the space . The first is the usual monomials , , , and so on. In this case, the coefficients are exactly the coefficients Matlab uses to specify a polynomial. The second is the set of Legendre polynomials, which will yield the same approximations but will turn out to have better numerical behavior. The third selection is the trigonometric functions, and the final selection is a set of piecewise constant functions. Advantages and disadvantages of each for numerical computations will be presented.
Once we have our basis set, we will consider how we can determine the approximating function as the ``best possible'' approximate for the given basis functions, and we will look at the behavior of the approximation error. Since we are working in , we will use the norm to measure error.
It turns out that approximation by monomials results in a matrix similar to the Hilbert matrix whose inversion can be quite inaccurate, even for small sizes. This inaccuracy translates into poor approximations. Use of orthogonal polynomials such as the Legendre polynomials, results in a diagonal matrix that can be inverted almost without error, but the right side can be difficult and slow to compute to high accuracy. In addition, roundoff accumulation in function reconstruction can be a serious problem. Fourier approximation substantially reduces the roundoff errors, but is slow to compute and evaluate, although fast methods (not discussed in this lab) can improve speed dramatically. Approximation by piecewise constants is not subject to these sources of error until ridiculously large numbers of pieces are employed, but can be slow to converge.
We will be attempting to approximate several functions in this lab, all on the interval [-1,1]. These functions include:
function y=partly_quadratic(x) % y=partly_quadratic(x) % input x (possibly a vector or matrix) % output y, where % for x<=0, y=0 % for x>0, y=x(1-x) y=(heaviside(x)-heaviside(x-1)).*x.*(1-x);Remark: The ``heaviside'' function (sometimes called the ``unit step function'') is named after Oliver Heaviside and is defined as
function y=sawshape9(x) % y=sawshape9(x) % input x (possibly a vector or matrix) % output y, where % y=(x+1) for -1<=x<0 % y=0 for x=0 % y=(x-1) for 0<x<=1 y= heaviside(-x).*(x+1) + heaviside(x).*(x-1);
Remark: In Lab 6, you saw a similar function defined using the Matlab find command, while heaviside is used here. There is no essential reason for this change.
This lab will take four sessions. If you print this lab, you may prefer to use the pdf version.
This kind of approximation requires evaluation of integrals. We will use the Matlab numerical integration function integral. Since we will be computing fairly small errors, will replace the default error tolerances with smaller ones. The Matlab command is
q=integral(func,-1,1,'AbsTol',1.e-14,'RelTol',1.e-12);where func is a function handle to a function written using vector (array) syntax. This command will result in an approximation, , satisfying
q=quadgk(func,-1,1,'AbsTol',1.e-14,'RelTol',1.e-12, ... 'MaxIntervalCount',15000);
Since these integration functions involve extra parameters and extra typing, You should create an abbreviation for them for use in this lab. Create an m-file named ntgr8 in the following way:
function q=ntgr8(func) % q=ntgr8(func) integrates func over [-1,1] with tight tolerances q=integral(func,-1,1,'AbsTol',1.e-14,'RelTol',1.e-12);with a similar function if you are using quadgk. This function will save you some typing for this lab. The name looks strange in print, but it is easy to remember because it is a pun. The letter n can be pronounced ``en,'' the letter t can be pronounced ``tee'' and gr8 can be pronounced ``grate.'' Put them together and it sounds like ``integrate.''
The problem of approximation can be described in the following way.
Problem: Given a function , and a (complete) set of functions , , then for a given , find the set of values so that
Quarteroni, Sacco, and Saleri, in Section 10.7, discuss least-squares approximation in function spaces such as . The idea is to minimize the norm of the difference between the given function and the approximation. Given a function and a set of approximating functions (such as the monomials ), for each vector of numbers define a functional
Consider the functional
Equation (5) is related to Equations (10.1) and (10.2) in Quarteroni, Sacco, and Saleri, but their presentation focusses on orthogonal polynomials. For an arbitrary value of , Equation (5) can be written in the following way, where the indexing corresponds with Matlab indexing (starting with 1 instead of 0) and with the Matlab convention for coefficients of polynomials (first coefficient is for the highest power of .)
Examining Equation (6), there are two sets of integrals that need to be evaluated in order to compute the coefficients . On the left side, the integrands involve the products . On the right side, the integrands involve the products . Please note that these expressions are made more complicated by the fact that they are indexed ``backwards.'' This is done to be consistent with Matlab's numbering scheme for coefficients. Later in the lab when we switch to Legendre polynomials and are free to number the coefficients as we wish, we will switch to a simpler numbering scheme.
Once the coefficients have been found, the Matlab polyval function can be used to evaluate the resulting polynomials.
function c=coef_mon(func,n) % c=coef_mon(func,n) % func is a function handle % ... more comments ... % your name and the dateand be called coef_mon.m. It will solve the system (6) by constructing the matrix and right side and solving the resulting system for . In coef_mon.m, func refers to a function handle.
for k=1:n % force b to be a column vector with second index b(k,1)=ntgr8(@(x) func(x).*x.^(n-k)); endWarning: You should already have created the abbreviation function ntgr8.m according to instructions in Section 1.2 above on Matlab remarks contained in the introduction to this lab.
for k=1:n for ell=1:n H(k,ell)=ntgr8( ??? ) end endNote 1: I try not to use the letter ``l'' as a variable because it looks so much like the number 1. Instead, I use ell.
c=H\b;Do not be surprised if Matlab warns you that H is poorly conditioned for larger values of n.
func=@partly_quadratic; c=coef_mon(func,n); xval=linspace(-1,1,10000); yval=polyval(c,xval); yexact=func(xval); plot(xval,yval,xval,yexact) % relative Euclidean norm is approximating % the relative integral least-squares (L2 norm) % using an approximate trapezoid rule relativeError=norm(yexact-yval)/norm(yexact)and use it to evaluate the approximation for n=1 and n=5. Look at the plot and estimate by eye if the area between the exact and approximate curves is divided equally between ``above'' and ``below.'' Further, the error for n=5 should be smaller than for n=1 and the plot should look much better, but still far from perfect.
Runge n relative error 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 20 _______________ 30 _______________ 40 _______________
partly_quadratic n relative error 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 20 _______________ 30 _______________ 40 _______________
sawshape9 n relative error 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 20 _______________ 30 _______________ 40 _______________
You should notice that much smaller errors are associated with the smooth Runge example function than with the non-differentiable partly quadratic function and with the discontinuous sawshape9 function.
You should see that the errors do not seem to be decreasing much for the final values of , and wonder why the method becomes poor as gets large. It may not be obvious, but the matrices (5) and (6) are related to the Hilbert matrix and are extremely difficult to invert. The reason inversion is difficult is because the monomials all start to look the same as gets larger, that is, they become almost parallel in the sense. Even if the integration could be performed without error, you would observe roundoff errors in evaluating the resulting high-order polynomial.
One way to make the approximation problem easier might be to pick a better set of functions than monomials. The following section discusses a good alternative choice of polynomials. These polynomials allow much larger values of .
The Legendre polynomials form an -orthogonal set of polynomials. You will see below why orthogonal polynomials make particularly good choices for approximation. In this section, we are going to write m-files to generate the Legendre polynomials and we are going to confirm that they form an orthogonal set in . Throughout this section, we will be representing polynomials as vectors of coefficients, in the usual way in Matlab.
The Legendre polyonomials are a basis for the set of all polynomials,
just as the usual monomial powers of are.
They are appropriate for use on the interval [-1,1] because they
are orthogonal when considered as members of
Polynomials that are orthogonal are discussed
by Quarteroni, Sacco, and Saleri in Chapter 10, with
Legendre polynomials discussed in Section 10.1.2.
The first few Legendre polynomials are:
The value at of any Legendre polynomial
can be determined using the following recursion:
The following recursive Matlab function computes the values of the Legendre polynomial.
function yval = recursive_legendre ( k , xval ) % yval = recursive_legendre ( k , xval ) % yval = values of the k-th Legendre polynomial % at values xval if k<0 error('recursive_legendre: k must be nonnegative.'); elseif k==0 % WARNING: no space between else and if! yval = ones(size(xval)); elseif k==1 % WARNING: no space between else and if! yval = xval; else yval = ((2*k-1)*xval.*recursive_legendre(k-1,xval) ... - (k-1)*recursive_legendre(k-2,xval) )/k; end
Unfortunately, this recursive function is too slow to be used in this lab. The alternative to recursive calculation of Legendre polynomials is one that uses loops. It is a general fact that any recursive algorithm can be implemented using a loop. The code for the loop is typically more complicated than the recursive formulation. In the following exercise, you will write an algorithm using loops for Legendre polynomials.
function yval = legen ( n , xval) % yval = legen ( n , xval) % more comments % your name and the dateand add appropriate comments.
xval=linspace(0,1,10); norm( legen(3,xval)-(5*xval.^3-3*xval)/2 )and showing that the difference is of roundoff size.
xval=linspace(0,1,20); norm( legen(10,xval) - recursive_legendre(10,xval) )The difference should be of roundoff size.
The Legendre polynomials form a basis for the linear space of polynomials. One good characteristic of any set of basis vectors is to be orthogonal. For functions, we use the standard dot product, and say that two functions and are orthogonal if their dot product
q=ntgr8(@(x) ffunc(x).*gfunc(x) );where ffunc gives the function and gfunc gives .
Legendre polynomial approximation in follows the same recipe as monomial approximation:
function d=coef_legen(func,n) % d=coef_legen(func,n) % comments % your name and the dateto compute the coefficients of the approximation as . You should use the ntgr8 abbreviation function in the same way as used in Exercise 2 Note: The factor comes from the inverse of the diagonal matrix .
function yval=eval_legen(d,xval) % yval=eval_legen(d,xval) % comments % your name and the date
Runge n relative error 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 20 _______________ 30 _______________ 40 _______________ 50 _______________
partly_quadratic n relative error elapsed time 5 _______________ 10 _______________ 20 _______________ 40 _______________ 80 _______________ 160 _______________ _______________ 320 _______________ _______________ 640 _______________ _______________
sawshape9 n relative error elapsed time 5 _______________ 10 _______________ 20 _______________ 40 _______________ 80 _______________ 160 _______________ _______________ 320 _______________ _______________ 640 _______________ _______________
You should find the same values as for approximation by monomials for small n, and you can accurately compute with larger values of n using Legendre polynomials than using monomials. However, using large values of n can result in computing times that grow more rapidly than expected as n increases because the integral function must compensate for roundoff error arising from rapid oscillations. In fact, the time does grow more rapidly than , where you just estimated .
There is another set of functions that is orthogonormal in the sense. This is the set of trigonometric functions
The sum of the first terms of the Fourier series for a function is given as
function [z,s,c]=coef_fourier(func,n) % [z,s,c]=coef_fourier(func,n) % more comments % your name and the dateto compute the first coefficients of the Fourier series using Equation (10).
function yval=eval_fourier(z,s,c,xval) % yval=eval_fourier(z,s,c,xval) % more comments % your name and the date
Runge n relative error elapsed time 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 50 _______________ 100 _______________ ______________ 200 _______________ ______________ 400 _______________ ______________ 800 _______________ ______________
partly_quadratic n relative error elapsed time 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 50 _______________ 100 _______________ ______________ 200 _______________ ______________ 400 _______________ ______________ 800 _______________ ______________
sawshape9 n relative error elapsed time 1 _______________ 2 _______________ 3 _______________ 4 _______________ 5 _______________ 6 _______________ 10 _______________ 50 _______________ 100 _______________ ______________ 200 _______________ ______________ 400 _______________ ______________ 800 _______________ ______________
We have learned that approximation is best done using matrices that are easy to invert accurately, like diagonal matrices. This is the reason for using sets of orthogonal basis functions. We would also like to be able to perform the right side integrals easily as well. A large part of the reason that the orders of approximations in the exercises above have been restricted is that the integrals are difficult to perform accurately because they ``wiggle'' a lot, a major source of inaccuracy in the approximation. Using the integral function hides the inaccuracy, but you pay for it because the integrations require substantial time for higher values of n, as you may have noticed.
In this section, we will look at approximation by piecewise constant functions. Approximation by piecewise linears or higher are also useful, but all the important steps are covered with piecewise constants. Furthermore, piecewise constants are easy to extend to higher dimensions.
Suppose that a number is given and that the interval is divided into equal subintervals and points , . For , a function can be defined as
If a vector of coefficients can be found to represent the piecewise constant approximation to a function , then the approximation can be evaluated as
In the following exercise, we will follow the same basic recipe as before to compute the coefficients and the approximation to .
function a=coef_pc(func,Npc) % a=coef_pc(func,Npc) % comments % your name and the dateto compute the coefficients of the approximation as
function yval=eval_pc(a,xval) % yval=eval_pc(a,xval) % comments % your name and the dateAs in coef_pc, you will need to use linspace to generate the points , and bracket to find the values of corresponding to the values of xval.
xval=[-0.95,-0.65,-0.45,-0.25,-0.05,0.15,0.35,0.55,0.75,0.95]Be sure your answers are correct before continuing.
Look carefully and critically at the plot for the Runge function with Npc=8. You should be able to justify to yourself visually that no other piecewise constant function would produce a better approximation. Send me the plot for the Runge function with Npc=8.
Runge example Npc relative error elapsed time 4 _______________ 8 _______________ 16 _______________ 64 _______________ 256 _______________ 1024 _______________ ______________ 4096 _______________ ______________ 16384 _______________ ______________
partly quadratic Npc relative error elapsed time 4 _______________ 8 _______________ 16 _______________ 64 _______________ 256 _______________ 1024 _______________ ______________ 4096 _______________ ______________ 16384 _______________ ______________
sawshape9 Npc relative error elapsed time 4 _______________ 8 _______________ 16 _______________ 64 _______________ 256 _______________ 1024 _______________ ______________ 4096 _______________ ______________ 16384 _______________ ______________
Remark: As you might imagine, approximation using piecewise linear functions will converge more rapidly than using piecewise constants. There are alternative approaches for using piecewise linears: piecewise linear functions on each interval with jumps at interval endpoints, as the piecewise constant functions have; and, piecewise linear functions that are continuous throughout whole interval. The first retains orthogonality and the diagonal form of the coefficient matrix . The second sacrifices the diagonal form for a banded form that is almost as easy to solve, as you may see in the next exercise. Continuity, however, can be worth the sacrifice, depending on the application. Even higher order piecewise polynomial approximation is possible, if the application can benefit.
Remark: In Exercises 3, 4, and 5 you were asked to estimate the growth of the running time with n. You found in Exercises 3 and 4 that the running time increased more rapdily than linearly (``superlinear'') as n increases, but in Exercise 5 it increased linearly. The distinction between superlinear and linear is important when choosing an algorithm to use. When running time increases too fast, algorithms can become too time-consuming to be useful. Since it is clear that approximation algorithms must scale at least linearly in n, Exercise 5 shows that piecewise constant approximation can be an attractive choice. In the extra credit problem below, you will find a similar result for piecewise linear approximation.
Piecewise linear approximations improve the rate of convergence over piecewise constant approximations, at the cost of increased work. In addition, piecewise linear approximations are commonly used in finite element approximations to differential equations. In this extra credit exercise, you will see how the same approach you saw above can be extended to piecewise linear approximations.
For this presentation, you again break the interval into equal subintervals. Denote the number of intervals by , although it is the same as above. Thus, there are points defined over the interval according to the Matlab function x=linspace(-1,1,Npl+1). For each , define a set of ``hat'' functions as
There is an equation analogous to (4), (8), (9), and (12):
Back to MATH2070 page.