Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 45
Текст из файла (страница 45)
. . , k - 1, for the three-term recurrence relation (6.30) satisfied by the orthogonal polynomialsP0(x), . . . , Pk(x); given also the constant α0 = P0(x), the coefficientsd0, . . . , dk of p(x) in (6.31), and a pointIf k = 0, then EXITIf k = 1, then EXIT258 APPROXIMATIONThen, on EXIT,is given byFORTRAN implementations of this algorithm have to contend withthe minor difficulty that some FORTRAN dialects do not allow zerosubscripts. Also, storage requirements and the number of necessary calculations vary from one set of orthogonal polynomials to another.Example 6.10 Write a FORTRAN implementation of Algorithm 6.1 in case the orthogonal polynomials are the Chebyshev polynomials.In this case, the Ai, Bi, Ci need not be stored in arrays since they do not depend oni. Also, the calculation of requires onlyhence it is not necessary toandstore the full arrayThe FORTRAN FUNCTION CHEB below solves the given problem.
NTERMS isthe number of terms in p(x); that is, p(x) is of degree < NTERMS - 1. Both NTERMSand the coefficientsD(i) = d i - 1i - l, . . . , NTERMSare assumed to be in the labeled COMMON POLY.REAL FUNCTION CHEB (X)C RETURNS THE VALUE OF THE POLYNOMIAL OF DEGREE . LT. NTERMS WHOSEC CHEBYSHEV COEFFICIENTS ARE CONTAINED IN D .INTEGER NTERMS, KREAL D,X, PREV,PREV2,TWOXCOMMON /POLY/ NTERMS,D(30)IF (NTERMS .EQ. 1) THENCHEB = D(1)RETURNEND IFTWOX = 2.*XPREV2 = 0.PREV = D(NTERMS)IF (NTERMS .GT. 2) THENDO 10 K=NTERMS-1,2,-lCHEU = D(K) + TWOX*PREV - PREV2PREV2 = PREVPREV = CHEB10CONTINUEEND IFCHEB = D(1) + X*PREV - PREV2RETURNENDEXERCISES6.3-l Using the appropriate recurrence relation, generate the first five Laguerre polynomials(for α = 0).6.3-2 Find the zeros of the Legendre polynomials P2(x), P3(x), and P4(x).6.3-3 Find the zeros of the Hermite polynomials H2(x), H3(x), H4(x).6.3-4 Express the polynomial p(x) = x4 + 2x3 + x2 + 2x + 1 as a sum of Legendre polynomials.*6.4LEAST-SQUARES APPROXIMATION BY POLYNOMIALS2596.3-5 Verify directly that the Legendre polynomial P3(x) is orthogonal to any polynomial ofdegree 2.6.3-6 Prove that if Pk(x) is the Legendre polynomial of degree k, thenUse the three-term recurrence relation satisfied by Legendre polynomials.6.3-7 Let P0(x), P1(x), .
. . be a sequence of orthogonal polynomials and let x0, . . . , xk bethe k + 1 distinct zeros of Pk+1 (x). Prove that the Lagrange polynomialsx j )/( xi - xi), i = 0, . . . , k, for these points are orthogonal to each other. [ Hint: Show thatfor ij, li(x)lj(x) = Pk+1(x)g(x), where g(x) is some polynomial of degree < k. ]*6.4 LEAST-SQUARES APPROXIMATION BYPOLYNOMIALSIn this section, we discuss the use of sequences of orthogonal polynomialsfor the calculation of polynomial (weighted) least-squares approximations.Let f(x) be a function defined on some interval (a,b), and supposethat we wish to approximate f(x) on (a,b) by a polynomial of degree < k.If we measure the difference between f(x) and p(x) by(6.33)where the scalar product is given by (6.26) or (6.27), then it is natural toseek a polynomial of degree < k for which (6.33) is as small as possible.Such a polynomial is called a (weighted) least-squares approximation tof(x) by polynomials of degree < k.The problem of finding such a polynomial is solved in Sec.
6.2 for theparticular case that the scalar product is given by (6.27) with the weightfunction w(x) = 1. In the general case, one proceeds as follows: Supposethat we can find, for the chosen scalar product, a sequenceP0 (x), .
. . , Pk (x) of orthogonal polynomials. By Property 1 of suchsequences (see Sec. 6.3), every polynomial p(x) of degree < k can bewritten in the formp(x) = d0 P0 (x) + · · · + dkPk(x)for suitable coefficients d0, . . . , dk. Substituting this into (6.33) it followsthat we want to minimizeE(d0 , .
. . , dk) = <f(x) - d0 P0 (x) - · · · - d kPk(x),f(x) - d 0 P 0 (x) - · · · - d k P k (x)>260APPROXIMATIONover all possible choices of d0, . . . , dk. Proceeding as in Sec. 6.2, oneshows that “best” coefficients d0*, . . . , dk* must satisfy the normal equationsi = 0. . . . .kd0*<P0, Pi > + d1*<P1, Pi > + · · · + dk*<Pk, Pi > = <f, Pi >which, because of the orthogonality of the Pj(x), reduce todi *<Pi , Pi > = <f, Pi >i = 0, . . .
, kHence, ifSi = <Pi , Pi >i = 0, . . . , kare all nonzero, then the best coefficients are simply given byi = 0, . . . , k(6.34)Example 6.11 Calculate the polynomial of degree < 3 which minimizesover all polynomials p(x) of degree < 3.In this case, f(x) = ex, and the scalar product is given byFrom Example 6.6, we find the orthogonal polynomials for this scalar product to be theLegendre polynomials. Using Table 6.2 of these polynomials, we calculateOne can show that, for the Legendre polynomials (see Exercise 6.3-6),so that S 0 - 2,Using (6.34) to calculate the d i * and usinge = 2.71828183, we find that the least-squares approximation to ex on (-1, 1) by cubicpolynomials isp * (x) = 1.175201194P0(x) + 1.103638324P1(x) + 0.3578143506P2 (x)+ 0.07045563367P3 (x)If we replace Pi(x) by their equivalent expressions in powers of x using Table 6.2and rearrange, we obtainp*(x) = 0.9962940183 + 09979548730x + 0.5367215260x2 + 0.1761390842x 3On (-1, 1), this polynomial has a maximum deviation from ex of about 0.011.*6.4LEAST-SQUARES APPROXIMATION BY POLYNOMIALS261If the appropriate orthogonal polynomials cannot be found in tables,one has to generate them.
This can be done with the aid of the three-termrecurrence relation (6.30). We now give an algorithmic description of thistechnique for the practically important case when the scalar product is(6.35)with x1, . . . , xN Certain fixed points in (a,b).Algorithm 6.2: Generation of orthogonal polynomials For simplicity,we elect to get all orthogonal polynomials with leading coefficient 1, sothatAi = αi = 1all iStep 0 Set P0(x) = 1. Further, calculateIf N > 1 and w(x) > 0, then S 0 is not zero, and we can go on tocalculateP 1 ( x ) = (x - B0)P0(x) = x - B0where, by Property 4 of orthogonal polynomials (see Sec.
6.3),With P0 (x), . . . , Pj(x) already constructed, the general, or jth,step proceeds as follows:Step j CalculateSince Pj(x) is a polynomial of exact degree j, Sj can be zero only if nomore than j of the points x1, . . . , xN are distinct. Hence, if there aremore than j distinct points among the xn ‘s, we can calculateand get the next orthogonal polynomial asP j + 1 ( x ) = (x - Bj)Pj(x) - CjPj-1(x)(6.36)262APPROXIMATIONExample 6.12 Solve the least-squares approximation problem of Example 6.5 usingorthogonal polynomials.For this example, f(x) = 10 - 2x + x2 /10,n = 1, .
. . ,6and we seek the polynomial of degree < 2 which minimizesi.e., we are dealing with the scalar product (6.27) with w(x) = 1. Following the Algorithm 6.2. we calculateP0(x) = 1henceThereforeP1 (x) = x - 10.5and, as S10, we can go on to calculate P2(x). We getif we carry seven decimal places and round. This givesP2 (x) - (x - 10.5)2 - 0.1166667S2 = 0.05973332Next, we calculate the best coefficients d0*, d1*, d2* for the least-squares approximationp * (x) = d 0 * P0 (x) + d 1 * P1 (x) + d 2 * P2 (x)using (6.34) and continuing with seven-decimal-digit floating-point arithmetic.
This givesTo compare this with the results computed in Example 6.5, we write p*(x)in terms of 1, x, x2. We getp*(x) = 0.03666667 + 0.1(x - 10.5)2+ 0.0999999[( x - 10.5) - 0.11666671= 0.03666667 - 1.05 + 0.0999999(110.25 - 0.1166667)+ [0.1 + 0.0999999(-21)]x + 0.0999999x 2*6.4LEAST-SQUARES APPROXIMATION BY POLYNOMIALS263Hence, computed this way, the ci * of Example 6.5 becomec1* = 9.99998 · · ·c2* = -1.9999998 · · ·c3* = 0.0999999 · · ·By contrast, we obtained in Example 6.5cl* = 8.492 · · ·c2* = -1.712 · · ·c3* = 0.0863 · · ·when we solved the normal equations (6.23) for the ci * ‘s directly, usingseven-decimal-digit floating-point arithmetic. The results using orthogonalpolynomials thus show an impressive improvement in this example.Incidentally, one would normally not go to the trouble of expressingp*(x) in terms of the powers of x. Rather, one would use Algorithm 6.1together with the computed d i * whenever p*(x) is to be evaluated, sinceone has the coefficients Bi and Ci of the recurrence relation available.In a FORTRAN implementation, the generation of the orthogonalpolynomials and the calculation of the best coefficients d i * are bestcombined into one operation to save storage.
For the calculation of dj* andof Pj+1(x), we only need the numbersP j (x n )Pj-1(xn )n=1,...,N*Hence, if d j is calculated as soon as Pj(xn ), n = 1, . . . , N, becomeavailable, then Pj(xn), n = 1, . . . , N, can safely be forgotten once Pj+1(x)and Pj+2(x) have been calculated. Again, there is no need to construct thePj(x) explicitly in terms of the powers of x, say, since we need only theirvalues at the xn, n = 1, . . . , N.SUBROUTINE ORTPOL ( X, F, W, NPOINT, PJMl, PJ, ERROR )C CONSTRUCTS: THE DISCRETE WEIGHTED LEAST SQUARES APPROXIMATION BY POLYC NOMIALS OF DEGREE LT. NTERMS TO GIVEN DATA.C****** I N P U T ******C (X(I), F(I)), I=l,...,NPOINT GIVES THE ABSCISSAE AND ORDINATES OFCTHE GIVEN DATA POINTS TO BE FITTED.C W NPOINT-VECTOR CONTAINING THE POSITIVE WEIGHTS TO BE USED.C NPOINT NUMBER OF DATA POINTS.C****** I N P U T VIA COMMON BLOCK P 0 L Y ******C NTERMS GIVES THE ORDER (= DEGREE + 1) OF THE POLYNOMIAL APPROXIMANT.C****** W O R K A R E A S ******C PJMl, PJ ARRAYS OF LENGTH NPOINT TO CONTAIN THE VALUES AT THE X'SCOF THE TWO MOST RECENT ORTHOGONAL POLYNOMIALS.C****** O U T P U T ******C ERROR NPOINT-VECTOR CONTAINING THE ERROR AT THE X'S OF THE POLYNOMIAL APPROXIMANT TO THE GIVEN DATA.C****** O U T P U T VIA COMMON BLOCK P 0 L Y ******C B, C ARRAYS CONTAINING THE COEFFICIENTS FOR THE THREE-TERM RECURRENCE.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
