Least Squares Approximation Hector D. Ceniceros 1 Least Squares Approximation Let f be a continuous function on [a, b]. We would like to find the best approximation to f by a polynomial of degree at most n in the L2 notm. We have already studied this problem for the approximation of periodic functions by trigonometric (complex exponential) polynomials. The problem is: find a poly Pn (x) of degree n such that Z (1) b [f (x) pn (x)]2 dx = min! a Such polynomial, is also called the Least Squares approximation to f . As an illustration let us consider n = 1. We look for P1 (x) = a0 + a1 x, for x [a, b] which minimizes Z b [f (x) P1 (x)]2 dx E(a0 , a1 ) = a (2) Z b Z b Z b 2 = f (x)dx 2 f (x)(a0 + a1 x)dx + (a0 + a1 x)2 a a a These are lecture notes for Math 104 A. These notes and all course materials are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 E(a0 , a1 ) is a quadratic function of a0 and a1 . Necessarily, the minimum is the critical point Z b Z b E(a0 , a1 ) f (x)dx + 2 (a0 + a1 x)dx = 0, = 2 a0 a a Z b Z b E(a0 , a1 ) xf (x)dx + 2 (a0 + a1 x)x dx = 0, = 2 a1 a a which yields the following linear 2 2 system for a0 and a1 : \u0012Z b \u0013 \u0012Z b \u0013 Z b 1dx a0 + (3) xdx a1 = f (x)dx, a a a \u0013 \u0013 \u0012Z b \u0012Z b Z b 2 (4) x dx a1 = xf (x)dx. xdx a0 + a a a These two equations are known as the Normal Equations for n = 1. Example 1. Let f (x) = ex for x [0, 1]. Then Z 1 (5) ex dx = e 1, Z 10 (6) xex dx = 1, 0 and the normal equations are 1 a0 + a1 = e 1, 2 1 1 a0 + a1 = 1, 2 3 (7) (8) whose solution is a0 = 4e 10, a1 = 6e + 18. Therefore the least squares approximation to f (x) = ex by a linear polynomial is P1 (x) = 4e 10 + (18 6e)x. P1 and f are plot in Fig. 1 The Least Squares Approximation to a function f on an interval [a, b] by a polynomial of degree at most n is the best approximation of f in the L2 norm by that class of polynomials. It is the polynomial Pn (x) = a0 + a1 x + + an xn 2 2.8 P1(x) x f(x)=e 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 1: The function f (x) = ex on [0, 1] and its Least Squares Approximation P1 (x) = 4e 10 + (18 6e)x. such that Z (9) b [f (x) Pn (x)]2 dx = min . a Defining this squared L2 error as E(a0 , a1 , ..., an ) we have Z b [f (x) (a0 + a1 x + + an xn )]2 dx E(a0 , a1 , ..., an ) = a Z b 2 f (x)dx 2 = a n X Z aj j x f (x)dx + a j=0 b n X n X j=0 k=0 Z aj ak b xj+k dx. a E is a quadratic function of the parameters a0 , ..., an . A necessary condition 3 is that the set of a0 , ..., an which minimizes E is the critical point. That is E(a0 , a1 , . . . , an ) am Z b Z b Z b n n X X m j+m x f (x)dx + aj x dx + xk+m dx = 2 ak 0 = a Z = 2 j=0 n X b xm f (x)dx + 2 a a k=0 xj+m dx, aj a j=0 and we get the Normal equations Z b Z b n X j+m (10) xm f (x)dx, aj x dx = j=0 a b Z m = 0, 1, ..., n. a a We can write the Normal Equations in matrix as Rb Rb n Rb Rb f (x)dx 1dx xdx x dx a 0 a a R ba 2 R ba n+1 a1 R b xf (x)dx R b xdx x dx x dx a a a (11) . . .. = a . . . . . . . . . . . . . . . Rb n Rb n R b n+1 R b 2n an x f (x)dx x dx a x dx a x dx a a The matrix in this system is symmetric. It is also positive definite and hence nonsingular so that there is a unique solution. Indeed, denoting this matrix by H and a = [a0 a1 an ]T any n + 1 row vector we have T a Ha = = = n X n X i=0 j=0 n X n X Z ai aj i=0 j=0 n X n Z b X i=0 j=0 = ai aj Hij a = i=0 j=0 n b X a ai xi aj xj dx ai xi aj xj dx aj xj )2 dx 0 ( a xi+j dx a Z bX n X n Z b j=0 4 Moreover, aT Ha = 0 if and only if a = 0, i.e. H is positive definite. Fur2E is equal to 2H so it is also positive definite and thermore, the Hessian, ai a j therefore the critical point is indeed a minimum. In the interval [0, 1], 1 1 1 1 2 n+1 1 1 1 n+2 2 3 (12) H = .. .. .. , . . . . . . 1 n+1 1 n+2 1 2n+1 which is known as the [(n + 1) (n + 1)] Hilbert matrix. In principle, the process to obtain the Least Squares Approximation Pn (x) to f is to solve the normal equations (10) to find the coefficients a0 , a1 , . . . , an and set Pn (x) = a0 + a1 x + . . . an xn . There are however two problems with this approach: 1. It is difficult to solve this linear system numerically for large n because the matrix H is very sensitive to small perturbations. 2. If we want to increase the degree of the approximating polynomial we need to start over again and solve a larger set of normal equations. That is, we cannot use the a0 , a1 , . . . , an we already found. It is more efficient and easier to solve the Least Squares Approximation problem using the idea of orthogonality, as we did with approximation by trigonometric polynomials. Suppose that we have a set of polynomials (in [a, b]) {0 (x), 1 (x), ..., n (x)} such that j (x) is a polynomial of degree j. Then we can write any polynomial of degree at most n as a linear combination of 0 (x), 1 (x), ..., n (x). In particular, the Least Square Approximating polynomial Pn (x) can be written as n X Pn (x) = a0 0 (x) + a1 1 (x) + ... + an n (x) = aj j (x). j=0 5 Then Z n X b [f (x) E(a0 , ..., an ) = a Z j=0 n X b 2 f (x)dx 2 = (13) aj j (x)]2 dx a + j (x)f (x)dx aj a j=0 n X n X b Z b Z j (x)k (x)dx aj ak a j=0 k=0 and E 0= = 2 am Z b m (x)f (x)dx + 2 a n X Z b aj j=0 j (x)m (x)dx. a for m = 0, 1, . . . , n, which gives the normal equations (14) n X Z b aj j=0 b Z m (x)f (x)dx, j (x)m (x)dx = m = 0, 1, ..., n. a a Now if the set of approximating functions {0 (x), ....n (x)} is orthogonal, i.e. \u001a Z b Z b 0 j 6= m (15) j (x)m (x)dx = , j = 2j (x)dx > 0 j = m j a a then (16) 1 am = m Z b Z m (x)f (x)dx, m = b 2m (x)dx, m = 0, 1, ..., n. a a and Pn (x) = a0 0 (x) + a1 1 (x) + ... + an n (x). Note that if the set {0 (x), ....n (x)} is orthogonal, (16) and (13) imply the Bessel inequality (17) n X j a2j Z a j=0 6 b f 2 (x)dx. This inequality shows that if f is square integrable, i.e. if Z b f 2 (x)dx < , a P 2 j=0 j aj converges. then the series We can consider the Least Squares approximation for a class of linear combinations of orthogonal functions {0 (x), ..., n (x)} not necessarily polynomials (we saw an example of this with Fourier approximations. Rb For complex-valued functions orthogonality means a j (x)m (x)dx = 0 if j 6= m, where the bar denotes the complex conjugate). It is convenient to define a weighted L2 norm associated with the Least Squares problem \u0013 21 \u0012Z b 2 f (x)w(x)dx , kf kw,2 = (18) a 1 where w(x) 0 for all x (a, b) . A corresponding inner product is defined by Z b (19) f (x)g(x)w(x)dx. hf, gi = a Definition 1. A set of functions {0 (x), ..., n (x)} is orthogonal, with respect to the weighted inner product (19), if hj , k i = 0 for j 6= k. 1.1 Linear Independence and Gram-Schmidt Orthogonalization Definition 2. A set of functions {0 (x), ..., n (x)} defined on an interval [a, b] is said to be linearly independent if (20) a0 0 (x) + a1 1 (x) + . . . an n (x) = 0, for all x [a, b] then a0 = a1 = . . . = an = 0. Otherwise, it is said to be linearly dependent. Example 2. The set of functions {0 (x), ..., n (x)}, where j (x) is a polynomial of degree j for j = 0, 1, . . . , n is linearly independent on any interval [a, b]. For a0 0 (x) + a1 1 (x) + . . . an n (x) is a polynomial of degree at most n and hence a0 0 (x) + a1 1 (x) + . . . an n (x) = 0 for all x in a given interval [a, b] implies a0 = a1 = . . . = an = 0. Rb Rb More precisely, we will assume w 0, a w(x)dx > 0, and a xk w(x)dx < + for k = 0, 1, . . .. We call such a w and admissible weight function. 1 7 Given a set of linearly independent functions {0 (x), ..., n (x)} we can produce an orthogonal set {0 (x), ..., n (x)} by doing the Gram-Schmidt procedure: 0 (x) = 0 (x) 1 (x) = 1 (x) c0 0 (x), h1 , 0 i h0 , 0 i 2 (x) = 2 (x) c0 0 (x) c1 1 (x), h2 , 0 i h2 , 0 i = 0 c0 = h0 , 0 i h2 , 1 i h2 , 1 i = 0 c1 = h1 , 1 i .. . h1 , 0 i = 0 c0 = We can write this procedure recursively as (21) (22) 0 (x) = 0 (x), j (x) = j (x) j1 X ck k (x), k=0 1.2 ck = hj , k i . hk , k i Orthogonal Polynomials Let us take the set {1, x, . . . , xn } on a interval [a, b]. We can use the GramSchmidt process to obtain an orthogonal set {0 (x), ..., n (x)} of polynomials with respect to the inner product (19). Each of the j is a polynomial of degree j, determined up to a multiplicative constant (orthogonality is not changed). Suppose we select the j (x), j = 0, 1, . . . , n to be monic, i.e. the coefficient of xj is 1. Then k+1 (x) xk (x) = rk (x), where rk (x) is a polynomial of degree at most k. So we can write (23) k+1 (x) xk (x) = k k (x) k k1 (x) + k2 X cj j (x). j=0 Then taking the inner product of this expression with k and using orthogonality we get hxk , k i = k hk , k i 8 and k = hxk , k i . hk , k i Similarly, taking the inner product with k1 we obtain hxk , k1 i = k hk1 , k1 i but hxk , k1 i = hk , xk1 i and xk1 (x) = k (x)+qk1 (x), where qk1 (x) is a polynomial of degree at most k 1 then hk , xk1 i = hk , k i + hk , qk1 i = hk , k i, where we have used orthogonality in the last equation. Therefore k = hk , k i . hk1 , k1 i Finally, taking the inner product of (23) with m for m = 0, . . . , k 2 we get hk , xm i = cm hm , m i m = 0, . . . , k 2 but the left hand side is zero because xm (x) is a polynomial of degree at most k 1 and hence it is orthogonal to k (x). Collecting the results we obtain a three-term recursion formula (24) 0 (x) = 1, (25) 1 (x) = x 0 , (26) (27) (28) 0 = hx0 , 0 i h0 , 0 i for k = 1, . . . n k+1 (x) = (x k )k (x) k k1 (x), hxk , k i k = , hk , k i hk , k i k = . hk1 , k1 i Example 3. Let [a, b] = [1, 1] and w(x) 1. The corresponding orthogonal polynomials are known as the Legendre Polynomials and are widely used in a variety of numerical methods. Because xk2 (x)w(x) is an odd function 9 it follows that k = 0 for all k. We have 0 (x) = 1 and 1 (x) = x. We can now use the three-term recursion (26) to obtain R1 2 x dx = 1/3 1 = 1 R1 dx 1 and 2 (x) = x2 31 . Now for k = 2 we get R1 1 2 = (x2 31 )2 dx = 4/15 R1 2 dx x 1 4 and 3 (x) = x(x2 13 ) 15 x = x3 53 . We now collect Legendre polynomials we found: 0 (x) = 1, 1 (x) = x, 1 2 (x) = x2 , 3 3 3 (x) = x3 x, 5 .. . Theorem 1. The zeros of orthogonal polynomials are real, simple, and they all lie in [a, b]. Proof. Indeed, k (x) is orthogonal to 0 (x) = 1 for each k 1, thus Z (29) b k (x)w(x)dx = 0 a i.e. k has to change sign in [a, b] so it has a zero, say x1 (a, b). Suppose x1 is not a simple root, then q(x) = k (x)/(x x1 )2 is a polynomial of degree k 2 and so Z b k2 (x) 0 = hk , qi = w(x)dx > 0, 2 a (x x1 ) which is of course impossible. Assume that k (x) has only l zeros in (a, b), x1 , . . . , xl . Then k (x)(x x1 ) (x xl ) = qkl (x)(x x1 )2 (x xl )2 , 10 where qkl (x) is a polynomial of degree k l which does not change sign in [a, b]. Then Z b qkj (x)(x x1 )2 (x xl )2 w(x)dx 6= 0 hk , (x x1 ) (x xl )i = a but hk , (x x1 ) (x xl )i = 0 for l < k. Therefore l = k. 1.3 Chebyshev Polynomials There is a class of orthogonal polynomials that are very important in applications (integration, spectral methods for PDE's, etc.) as well as in other areas of mathematics (number theory, algebra, etc.) We take the interval [a, b] = [1, 1]. A general interval [a, b] can be obtain from a simple change of variables = 21 (a + b) 21 (a b)x. For each x [1, 1] there is a unique [0, ] such that x = cos . Define (30) Tn (x) = cos(n cos1 (x)) for n = 0, 1, . . . Although not immediately apparent, this is a polynomial in x of degree n, as we will see shortly. From the definition (30) (31) Tn (cos ) = cos n, [0, ]. Using the trigonometry identity (32) cos[(n + 1)] + cos[(n 1)] = 2 cos n cos , we immediately get (33) Tn+1 (cos ) + Tn1 (cos ) = 2Tn (cos ) cos and going back to the x variable we obtain the recursion formula (34) T0 (x) = 1, Tn+1 (x) = 2xTn (x) Tn1 (x), n 1. It is clear from the recursion formula that the Tn (x) are indeed polynomials. Let us generate a few of them. T0 (x) = 1, T1 (x) = cos(1 cos1 (x)) = x, T2 (x) = cos(2 cos1 (x)) = 2x x 1 = 2x2 1, T3 (x) = 2x (2x2 1) x = 4x3 3x, T4 (x) = 2xT3 (x) T2 (x) = 8x4 8x2 + 1. 11 From these few Chebyshev polynomials a pattern is observed. The Tn (x) is a polynomial of exactly degree n and its leading order coefficient is 2n1 . Moreover, Tn (x) is an even (odd) function of x if n is even (odd). These polynomials are orthogonal with respect to the weight function 1 . (x) = 1 x2 (35) Indeed: Z 1 Z Tn (x)Tm (x)(x)dx = 1 1 cos(n cos1 (x)) cos(m cos1 (x)) 1 1 dx 1 x2 1 x = cos , d = dx 2 1 x Z cos(n) cos(m)d = 0 and because 1 cos(n) cos(m) = [cos(m + n) + cos(m n)], 2 we get \u0015\f \u0014 Z 1 1 1 1 Tn (x)Tm (x)(x)dx = sin((n + m)) + sin(n m) 2 n+m nm 1 0 = 0, if m 6= n. Consequently Z 1 (36) Z Tn (x)Tm (x)(x)dx = 1 cos(m) cos(n)d = 0 0 m 6= n m=n>0 m=n=0 2 Given that Tn (cos ) = cos(n) the zeros of Tn occurs when n is an odd multiple of /2. Since [0, ], we get the n zeros at k = (2k1) , for 2n k = 1, 2, . . . , n. This corresponds to the points \u0012 \u0013 (2k 1) (37) xk = cos(k ) = cos , k = 1, ..., n. 2n 12 Tn (x) is oscillating between 1 and 1, as Tn (cos ) = cos(n). Thus, the extrema of Tn correspond to the values k such that cos nk = 1, i.e. k = k/n for k = 0, 1, ..., n. Therefore, the max and min of Tn are attained at \u0012 \u0013 k (38) , k = 0, 1, 2..., n, xk = cos n and Tn (xk ) = (1)k . These xk are the Chebyshev (also known as GaussLobatto) nodes we have used in interpolation. Let us now consider the monic Chebyshev Polynomials, Tn (x), n = 0, 1, . . .. These can be obtain by dividing Tn by its leading order constant 2n1 : (39) T0 (x) = 1, (40) Tn (x) = 1 2n1 Tn (x), n 1. As n increases the extrema of Tn (x) decreases and the polynomial equioscillates between its max and min. These particular monic polynomials are beautifully connected with the error in interpolation. Recall that f (x) Pn (x) = 1 f (n+1) ((x))(x x0 )(x x1 ) (x xn ) (n + 1)! We have no control over the term f (n+1) ((x)) but we can select the interpolation nodes x0 , x1 , . . . , xn to reduce |(x x0 )(x x1 ) (x xn )|. Note that (x x0 )(x x1 ) (x xn ) is a monic polynomial of degree n + 1. The question is then: out of all the monic polynomials of degree n + 1 which one has the smallest max in absolute value. Remarkably, the answer is Tn+1 (x)! Theorem 2. For m 1 (41) max |Tm (x)| max | pm (x)| x[1,1] x[1,1] for all monic polynomials pm (x) and equality occurs only if pm = Tm . Proof. Let pm (x) be a monic polynomial of degree m and suppose (42) max | pm (x)| max |Tm (x)| = x[1,1] x[1,1] 13 1 2m1 . The polynomial p(x) = Tm (x) pm (x) is of degree at most m 1 and at the extremal points xk of Tm we have (43) p(xk ) = (1)k pm (xk ), 2m1 k = 0, 1, . . . , m. Then (42) implies p(xk ) 0 for k even and p(xk ) 0 for k odd. That is, p(x) changes sign or is equal to zero at least m times and so it has m zeros. But p(x) is of degree at most m1, thus p(x) 0 and consequently pm = Tm . Going back to the interpolation error question we have that the term |(xx0 )(xx1 ) (xxn )| would be minimized if we select the interpolation nodes to be the zeros of the n + 1 Chebyshev polynomials, i.e. \u0013 \u0012 2k + 1 (44) , k = 0, 1, . . . n. xk = cos 2(n + 1) The following theorem summarizes this observation. Theorem 3. Let PnT (x) be the interpolating polynomial of degree at most n of f with respect to the nodes (44) then (45) max |f PnT (x)| x[1,1] 2 1 max |f n+1 (x)|. 2n (n + 1)! x[1,1] Discrete Least Squares Approximation Suppose that we are given the data (x1 , f1 ), (x2 , f2 ), , (xN , fN ) obtained from an experiment. Can we find a simple function that can appropriately fit these data? Suppose that empirically we determine that there is an approximate linear behavior between fk and xk , k = 1, . . . , N . What is the straight line y = a0 + a1 x that best fits these data? The answer depends on how we measure the error, the deviations fk (a0 + a1 xk ), i.e. which norm we use for the error. The most convenient norm is the 2-norm (the sum of the squares) because we will end up with a linear system of equations to find a0 and a1 . Other norms will yield a nonlinear system for the unknown parameters. So the problem is: Find a0 , a1 which minimize (46) E(a0 , a1 ) = N X [fk (a0 + a1 xk )]2 . k=1 14 We can repeat all the Least Squares Approximation theory that we have studied at the continuum level except that integrals are replaced by sums! The conditions for the minimum (47) N X E =2 [fk (a0 + a1 xk )](1) = 0, a0 k=1 (48) N X E =2 [fk (a0 + a1 xk )](xk ) = 0 a1 k=1 produce the Normal Equations: (49) a0 N X 1 + a1 k=1 a0 (50) N X x k + a1 k=1 N X xk = N X k=1 k=1 N X N X x2k = k=1 fk , xk f k . k=1 Approximation by a higher order polynomial Pn (x) = a0 + a1 x + + an xn , n < N 1 is similarly done. In practice we would like n << N . We find the a0 , a1 , ..., an which minimize (51) E(a0 , a1 , ..., an ) = N X [fk (a0 + a1 xk + + an xnk )]2 . k=1 Proceeding as in the continuum case, we get the Normal Equations N X (52) aj j=0 That is a0 N X xkj+m x0k k=1 x1k xm k fk , N X x1k + + an k=1 + a1 m = 0, 1, ..., n. k=1 + a1 k=1 a0 = k=1 N X N X N X N X x2k N X xnk = k=1 + + an k=1 N X k=1 .. . 15 xn+1 k N X x0k fk k=1 = N X k=1 x1k fk a0 N X xnk + a1 N X k=1 xn+1 + + an k k=1 N X x2n k = k=1 N X xnk fk . k=1 The matrix of coefficients of this linear system is, as in the continuum case, symmetric, positive definite, and highly ill-conditioned (small perturbations in the data can produce large variations in the solution). And again, if we want to increase the degree of the approximating polynomial we cannot reuse the coefficients we already computed for the lower order polynomial. But now we know how to go around these two problems: we use orthogonality! Let us define the (weighted) discrete inner product as (53) hf, giN = N X f (xk )g(xk )(xk ), (x) > 0. k=1 Let {0 (x), ..., n (x)} be a set of polynomials such that j (x) is of degree exactly j. Then the solution to the discrete Least Squares Approximation problem by a polynomial of degree n can be written as (54) Pn (x) = a0 0 (x) + a1 1 (x) + + an n (x) and the square of the error is given by (55) E(a0 , , an ) = N X [fk (a0 0 (xk ) + + an n (xk ))]2 (xk ). k=1 Consequently, the normal equations are (56) N X aj hj , m iN = hm , f iN , m = 0, 1, ..., n. j=0 If {0 (x), , n (x)} are orthogonal with respect to the inner product h, iN , i.e. if hj , m iN = 0 for j 6= m, then the coefficients of the Least Squares Approximation are given by (57) aj = hj , f iN , hj , j iN j = 0, 1, ..., n and Pn (x) = a0 0 (x) + a1 1 (x) + + + an n (x). 16 If the {0 (x), , n (x)} are not orthogonal we can produce an orthogonal set {0 (x), , n (x)} using the 3 term recursion formula adapted to the discrete inner product. We have (58) 0 (x) 1, (59) 1 (x) = x 0 , (60) (61) 0 = hx0 , 0 iN , h0 , 0 iN for k = 1, ..., n k+1 (x) = (x k )k (x) k k1 (x), hxk , k iN hk , k iN k = , k = . hk , k iN hk1 , k1 iN Then, aj = (62) hj , f iN , hj , j iN j = 0, ..., n. and the Least Squares Approximation is Pn (x) = a0 0 (x) + a1 1 (x) + + an n (x). Example 4. Suppose we are given the data: xk : 0, 1, 2, 3, fk = 1.1, 3.2, 5.1, 6.9 and we would like to fit to a line. the normal equations are (63) (64) a0 a0 4 X 1 + a1 4 X k=1 4 X k=1 4 X k=1 k=1 x k + a1 xk = x2k = 4 X k=1 4 X fk xk f k k=1 and performing the sums we have (65) (66) 4a0 + 6a1 = 16.3, 6a0 + 14a1 = 34.1. Solving this 2 2 linear system we get a0 = 1.18 and a1 = 1.93. Thus, the Least Squares Approximation is P1 (x) = 1.18 + 1.93x and the square of the error is E= 4 X [fk (1.18 + 1.93xk )]2 = 0.023. k=1 17 Example 5. Fitting to a exponential y = beaxk . Defining E(a, b) = N X [fk beaxk ]2 k=1 we get the conditions for a and b (67) N X E =2 [fk beaxk ](bxk eaxk ) = 0, a k=1 (68) N X E =2 [fk beaxk ](eaxk ) = 0. b k=1 which is a nonlinear system of equations. However, if we take the natural log of y = beaxk we have ln y = ln b + ax. Defining B = ln b the problem becomes linear in B and a. Tabulating (xk , ln fk ) we can obtain the normal equations (69) B B (70) N X 1+a N X k=1 N X k=1 N X k=1 k=1 xk + a xk = x2k = N X k=1 N X ln fk , xk ln fk , k=1 and solve this linear system for B and a. Then b = eB and a = a. If a is given and we only need to determine b then the problem is linear. From (68) we have b N X k=1 e 2axk = N X fk e axk PN b = Pk=1 N fk eaxk k=1 k=1 e2axk Example 6. Discrete orthogonal polynomials. Let us construct the first few orthogonal polynomials with respect to the discrete inner product with 1 and xk = Nk , k = 1, ..., N . Here N = 10 (the points are equidistributed in [0, 1]). We have 0 (x) = 1 and 1 (x) = x 0 , where PN xk hx0 , 0 iN = Pk=1 = 0.55. 0 = N h0 , 0 iN k=1 1 18 and hence 1 (x) = x 0.55. Now (71) (72) (73) 2 (x) = (x 1 )1 (x) 1 0 (x), PN 2 hx1 , 1 iN k=1 xk (xk 0.55) 1 = = P = 0.55, N 2 h1 , 1 iN (x 0.55) k k=1 h1 , 1 iN 1 = = 0.0825. h0 , 0 iN Therefore 2 (x) = (x 0.55)2 0.0825. We can now use these orthogonal polynomials to find the Least Squares Approximation by polynomial of degree at most two of a given set of data. Let us take fk = f (xk ) = x2k + 2xk + 3. Clearly the Least Squares Approximation should be P2 (x) = x2 + 2x + 3 ! Let us confirm this by using the orthogonal polynomials 0 , 1 and 2 . The Least Squares Approximation coefficients are given by (74) (75) (76) hf, 0 iN = 4.485, h0 , 0 iN hf, 1 iN a1 = = 3.1, h1 , 1 iN hf, 2 iN = 1, a2 = h2 , 2 iN a0 = which gives, P2 (x) = (x0.55)2 0.0825+(3.1)(x0.55)+4.485 = x2 +2x+3. 19 Homework # 6: Least Squares Approximation and Orthogonal Polynomials 1 1. Prove that any polynomial Pn (x) of degree at most n can be written as Pn (x) = n X aj j (x), j=0 where j (x) is a polynomial of degree exactly j, for j = 0, 1, . . . , n. 2. Suppose {0 (x), 1 (x), . . . , n (x)} is an orthogonal set of functions with respect to the L2 inner product, i.e. Z b hj , k i = j (x)k (x)dx = 0, if j 6= k. a Prove the Pythagorean theorem k0 + 1 + . . . n k2 = k0 k2 + k1 k2 + . . . kn k2 , where kf k2 = hf, f i. 3. The solution Pn (x) to the Least Squares Approximation problem of f by a polynomial of degree at most n is given explicitly in terms of orthogonal polynomials 0 (x), 1 (x), . . . , n (x), where j is a polynomial of degree j, by Pn (x) = n X aj j (x), j=0 aj = hf, j i . hj , j i (a) Let Pn be the space of polynomials of degree at most n. Prove that the error f Pn is orthogonal to this space, i.e. hf Pn , qi = 0 for each q Pn . (b) Using the analogy of vectors interpret this result geometrically (recall the concept of orthogonal projection). 4. (a) Obtain the first 4 Legendre polynomials in [1, 1]. b) Find the least squares polynomial approximations of degrees 1, 2, and 3 for the function f (x) = ex on [1, 1]. c) What is the polynomial least squares approximation of degree 4 for f (x) = x3 on [1, 1]? Explain. 1 All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials) and the intellectual content of the course itself are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 5. Plot the monic Chebyshev polynomials T0 (x), T1 (x), T2 (x), T3 (x), and T4 (x). 6. Prove that 2 Z 1 1 [T (x)]2 n dx = 1. 1 x2 (1) 7. The concentration c of a radioactive material decays according to the law c(t) = beat where t represents time in seconds, a = 0.1 sec1 , and b is the initial concentration. a)Using the Least Squares method and the data table (Table 1) below find b. b) Find the error in the least squares approximation. ti (sec) 1 2 3 4 Ci 0.91 0.80 0.76 0.65 Table 1: 8. Given a collection of data points {(xi , yi )}m i=1 find the best least squares approximation of the form y = ax2 + bx3 . 9. (a) Given a collection of data points {(xi , yi )}m i=1 find the best least squares approx2 imation of the form y = ax + bx . (b) Use this approximation to fit the data in Table 2. (c) Find the error in the least squares approximation. xi 1 2 3 4 yi 3.1 9.8 21.2 36.1 Table 2: 2 \f\f