Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 46
Текст из файла (страница 46)
WHICH GENERATES THE ORTHOGONAL POLYNOMIALS.CC D COEFFICIENTS OF THE POLYNOMIAL APPROXIMANT TO THE GIVEN DATA WITHRESPECT TO THE SEQUENCE OF ORTHOGONAL POLYNOMIALS.CTHE VALUE OF THE APPROXIMANT AT A POINT Y MAY BE OBTAINED BY ACCREFERENCE TOORTVAL(Y) .C****** M E T H 0 D ******C THE SECQUENCE P0, Pl, . . . . PNTERMS-l OF ORTHOGONAL POLYNOMIALS WITHC RESPECT TO THE DISCRETE INNER PRODUCTC(P,Q) = SUM ( P(X(I))*Q(X(I))*W(I) , I=l,...,NPOINT)C IS GENERATED IN TERMS OF THEIR THREE-TERM RECURRENCEPJPl(X) = (X - B(J+l))*PJ(X) - C(J+l)*PJMl(X) ,CC AND THE COEFFICIENT D(J) OF THE WEIGHTED LEAST SQUARES APPROXIMATC ION TO THE GIVEN DATA IS OBTAINED CONCURRENTLY ASCD(Jt1) = (F,PJ)/(PJ,PJ) , J=0,...,NTERMS-1.264CCCCAPPROXIMATIONACTUALLY, IN ORDER TO REDUCE CANCELLATION, (F,PJ) IS CALCULATED AS(ERROR, PJ) , WITH ERROR = F INITIALLY, AND, FOR EACH J , ERROR REDUCED BY D(J+l) *PJ AS SOON AS D(J+l) BECOMES AVAILABLE.*INTEGER NPOINT ,NTERMS,I, JREAL B,C,D,ERROR(NPOINT) ,F(NPOINT) ,PJ(NPOINT) ,PJMl(NPOINT),W(NPOINT) ,X(NPOINT),P,S (20)COMMON /POLY/ NTERMS,B(20) ,C(20) ,D(20)C9101112DO 9 J=l ,NTERMSB(J) = 0.D(J) = 3.S(J) = 0.C(1) = 0.DO 10 I=l,NPOINTD(l) = D(l) + F(I)*W(I)B(l) = B(l) + X(I)*W(I)S(l) = S(l) + W(I)D(l) = D(l)/S(l)CO 11 I=l,NPOINTERROR(I) = F(I) - D(l)IF (NTERMS .EQ.
l)B(l) = B(l)/S(l)DO 12 I=l,NPOINTPJMl(I) = l.PJ(I) = X(I) - B(l)RETURNC21222730DO 30 J=2,NTERMSDO 21 I=l,NPOINTP = PJ(I)*W(I)D(J) = D(J) + ERROR(I)*PP = P*PJ(I)B(J) = B(J) + X(I)*PS(J) = S(J) + PD(J) = D(J)/S(J)DO 22 I=l,NPOINTERROR (I) = ERROR(I) - D(J)*PJ(I)RETURNIF (J .EQ. NTERMS)B(J) = B(J)/S(J)C(J) = S(J)/S(J-l)DO 27 I=l,NPOINTP = PJ(l)PJ(I) = (X(I) - B(J))*PJ(I) - C(J) *PJMl (I)PJMl(I) = PCONTINUERETURNENDThe calculation of the D(j) as carried out in this subprogram needsperhaps some clarification. Since D(j) = d*j-1 we get from (6.34) that(6.37)whereas in the program, D(j) is calculated as(6.38)withERROR( n ) = fn - D(l)P0(xn) - · · · -D(j - l)Pj-2 (x n )all n(6.39)*6.4LEAST-SQUARES APPROXIMATION BY POLYNOMIALS265If one substitutes (6.39) into (6.38), one getsD (j)since Pj-1 is orthogonal to P0(x), .
. . , Pj-2(x). Hence, in exact or infiniteprecision arithmetic, both (6.37) and (6.38) give the same value for D(j).But in finite-precision arithmetic, (6.38) can be expected to be moreaccurate for the following reason: SinceP* r (x) = D(l)P0 (x) + · · · + D(r + l)Pr(x)is the (weighted) least-squares approximation to f(x) by polynomials ofdegree < r, it follows that the numbersERROR(n ) = fn - P*j-1(xn)can1, .lossthanbe. .ofisn = 1 , . . .
, NPOINTexpected to be of smaller size than are the numbers f n , n =, NPOINT. Hence the calculation of (6.38) is less likely to producesignificance due to subtraction of quantities of nearly equal sizethe calculation of (6.37) (see Exercise 6.4-l).Example 6.13 Given the values fn of f(x) = ex at xn = ( n - 1)/10 - 1(n = 1, . . . ,21),rounded to two places after the decimal point. Try to recover the information about f(x)contained in these data.We attempt to solve this problem by calculating the polynomial p * 3 (x) whichminimizesover all polynomials p3(x) of degree < 3.
The following FORTRAN program calculatesp*3(x) with the aid of the subprogram ORTPOL mentioned earlier, then evaluates p*3(x)at the xn using the FUNCTION ORTVAL, which is based on Algorithm 6.1.CPROGRAM FOR EXAMPLE 6.13 .PARAMETER NPMAX=l00I,J,NPOINTINTEGER NTERMS,REAL B, C, D, ERROR (NPMAX), F(NPMAX), PJ(NPMAX), PJMl(NPMAX), W(NPMAX)*,X(NPMAX)COMMON /POLY/ NTERMS,B(20),C(20),D(20)NPOINT = 21DO 1 I=l,NPOINTW(I) = 1.X(I) = -1.
+ FLOAT(I-1)/10.266APPROXIMATIONF(I) = FLOAT(IFIX(EXP(X(I))*l00. + .5))/100.NTERMS = 4CALL ORTPOL( X, F, W, NPOINT, PJMl, PJ, ERROR )PRINT 601, (J,B(J),C(J),D(J),J=l,NTERMS)601 FORMAT(I2,3E16.8)DO 60 I=l,NPOINTPJMl(T) = EXP(X(I))PJ(I) = ORTVAL(X(I))60PRINT 660, (X(I),F(I),PJ(I),ERROR(I),PJMl(I),I=l,NPOINT)660 FORMAT(F5.1,F8.3,Fl0.5,E13.3,F10.5)STOPEND1REAL FUNCTION ORTVAL (X)RETURNS THE VALUE AT X OF THE POLYNOMIAL OF DEGREE .LT. NTERMSGIVEN BYD(l)*P0(X) + D(2)*Pl(X) + . . . + D(NTLRMS)*PNTERMS-l(X),WITH THE SEQUENCE P0, Pl, . .
. OF ORTHOGONAL POLYNOMIALS GENERATEDBY THE THREE-TERM RECURRENCEPJPl(X) = (X - B(J+l))*PJ(X) - C(J+l)*PJMl(X) , ALL J .CCCCCCCCOMMON /POLY/ NTERMS,B(20),C(20),D(20)PREV = 0.ORTVAL = D(NTERMS)IF (NTERMS .EQ. 1)RETURNDO 10 K=NTERMS-l,l,-lPREV2 = PREVPREV = ORTVALORTVAL = D(K) + (X - B(K))*PREV - C(K + I)*PREV210 CONTINUE.RETURNENDTable 6.4 Computer results for Example 6.13xnfn-1.0 0.370-0.9 0.410-0.8 0.450-0.7 0.500-0.6 0.550-0.5 0.610-0.4 0.670-0.3 0.740-0.2 0.820-0.1 0.9000.0 1.0000.1 1.1100.2 1.2200.3 1.3500.4 1.4900.5 1.6500.6 1.8200.7 2.0100.8 2.2300.9 2.4601.0 2.720p * 3 (xn )0.363870.408740.454810.503150.554840.610940.675240.740700.816500.901010.995301.100441.217511.347581.491721.651001.826502.019292.230442.461022.71211fn - p * 3 (xn )6.130E1.263E-4.806E-3.148E-4.836E-9.436E-2.542E-7.045E3.497E-1.010E4.710E9.558E2.49OE2.422E-1.717E-1.000E-6.499E-9.287E-4.368E-1.020E7.890E-030303030304030403030303030303030303040303p*4(xn)0.371150.408740.450970.498040.550210.607890.671560.741830.819400.905070.999761.104501.220401.348711.490741647951.821882.014182.226602.461022.71939fn - p * 4 (xn )-1.154E1.263E-9.719E1.964E-2.134E2.108E-1.565E-1.832E6.029E-5.070E2.358E5.499E-4.045E1.294E-7.399E2.052E-1.876E-4.176E0.397E-1.020E6.061E-03 0.3678803 0.406570 4 0.4493303 0.4965904 0.5488103 0.6065303 0.6703203 0.7408204 0.8187303 0.9048404 1.0000003 1.1051704 1.2214003 1.3498604 1.4918203 1.6487203 1.8221203 2.0137503 2.2255403 2.459600 4 2.71828*6.4LEAST-SQUARES APPROXIMATION BY POLYNOMIALS267Figure 6.6 The error in the least-squares approximation to the data of Example 6.13 bypolynomials of degree (a) zero, (b) one, (c ) two, (d ) three, (e ) four, (f ) five.Table 6.4 gives the results of the calculations which were carried out on a CDC6500.
We have plotted the error, fn - p*3(xn), in Fig. 6.6 d, which shows the error tobehave in a somewhat regular fashion, suggesting the p*3(x) does not represent all theinformation contained in the given data. We therefore calculate also the least-squaresapproximation p*4(x) to the given data by polynomials of degree < 4. The results arealso listed in Table 6.4. The error fn - p*4(x) is plotted in Fig. 6.6 e, and is seen to behavequite irregularly.
Hence p*4(x) can be assumed to represent all the information containedin the given data fn. Increasing the degree of the approximating polynomial any furtherwould only serve to give the approximating function the additional freedom to approximate the noise in the data, too.EXERCISES6.4-l If f(x) = 6,000 + x, then any least-squares approximation to f(x) by straight lines isf(x) itself. Calculate the polynomialwhich minimizesNote that 1 and x are already orthogonal, so that one merely has to calculate d*0 and d*1.Show the difference between (6.37) and (6.38) by calculating d*1 both ways, using four-decimal-digit floating-point arithmetic.268APPROXIMATION6.4-2 Calculate the polynomial of degree < 2 which minimizesover all polynomials p(x) of degree < 2. Use Legendre polynomials and carry out allcalculations to five decimal places.
(Note: π = 3.141593.)6.4-3 Implement the subroutine ORTPOL on your computer. Then use this subroutine tosolve the following problem. From a table of values of f(x) = sin πx, find fn = sin πxn atx n = (n - 1)/10 - 1 (n = 1,. . . , 21). rounded off to three decimal places. Then find thepolynomial p*4(x) which minimizesover all polynomials p4(x) of degree < 4.*6.5 APPROXIMATION BY TRIGONOMETRICPOLYNOMIALSMany physical phenomena, such as light and sound, have periodic character. They are described by functions f(x) which are periodic, i.e., whichsatisfyfor all x and some fixed number τ, the period of the function. Since theonly periodic polynomials are the constant functions, one has to use otherfunction classes for the effective approximation of periodic functions, andthe trigonometric polynomials offer themselves as an appropriate alternative.A trigonometric polynomial of order n is, by definition, any function ofthe form(6.40)with a0, .
. . , an and b1, . . . , bn real or complex constants. Such a trigonometric polynomial is 2π-periodic. We would therefore have to make someadjustment when approximating a τ-periodic function f(x) withWe agree to consider in such a case the 2π-periodic function g(x) =Then, having constructed a trigonometric polynomial approximation p(x) to g(x), we obtain from it a τ-periodic approximation forf(x) in the formWith this, we will assume from now on that thefunction f(x) to be approximated is already 2π-periodic.As it turns out, it is often more convenient to write trigonometricpolynomials of order n in the equivalent complex form(6.41)*6.5APPROXIMATION BY TRIGONOMETRIC POLYNOMIALS269Here, and for the remainder of this section and the next, the symbol istands for the imaginary unit,and the connection between (6.40) and (6.41) is provided by Euler’sformulaeix = cos x + i sin x(6.42)(a proof of which can be found in Exercise 1.7-9).
From Euler’s formula,we find [with cos (-jx) = cos jx, sin (-jx) = -sin jx] thatThis shows that (6.41) is of the form (6.40) witha j = c i + c -jb j = i(c j - c - j )j = 0, . . . , n(6.43 a)This relationship is easily inverted to give that (6.40) is of the form (6.41)withcj = (aj - ibj)/2c - j = (aj + ibj)/2j = 0, . . .
, n (6.436)Note that (6.41) represents a real function if and only if it is its owncomplex conjugate. But, sincethis means that (6.41) is a real function if and only ifall j(6.44)Thus, if (6.40) or (6.41) is a real function, then (6.43a ) simplifies tob j = -2 Im cj(6.45)aj = 2 Re cjApproximation by trigonometric polynomials is dominated by theFourier series(6.46)with the Fourier coefficientscalculated by(6.47)This series converges to f(x) under rather mild conditions [but not forevery f(x)].
For example, the series converges uniformly if f(x) is continuous with a piecewise-continuous first derivative.270 APPROXIMATIONThe Fourier series derives from the following fact:This shows that the functions 1, e±ix , e±i2x , . . . are orthonormal withrespect to the scalar or inner productIn other words,This provesTheorem 6.4 The partial sumof the Fourier series for f(x) is the best approximation to f(x) bytrigonometric polynomials of order n with respect to the normFurther, it can be shown that Parseval’s relation(6.48)holds.The Fourier coefficientsfor the function f(x) are used to “understand” the function f(x), as follows.
Suppose f(x) is a real 2π-periodicfunction. If we think of f(x) as the position at time x of some objectmoving on a line, then our 2π-periodic function f(x) describes a periodicmotion. If now[the polar form for the complex numberseries for f(x) asthen we can write the Fourier(see Exercise 6.5-7). In this way, we have represented our periodic motion*6.5APPROXIMATION BY TRIGONOMETRIC POLYNOMIALS271described by f(x) as a sum or superposition of simple harmonic oscillations. The jth such motion,has amplitudefrequency j/(2π), angular frequency j, period orwavelength 2 π/j, and phase angle θj.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
