Thompson - Computing for Scientists and Engineers (523188), страница 20
Текст из файла (страница 20)
The numerical methods developedfor differentiation in this chapter are useful for differential-equation solution in Sections 7.4 and 7.5 (first-order equations) and in Sections 8.3 - 8.5 (second-orderequations). In Section 10.4 we compute the convolution of a Gaussian distributionwith a Lorentzian distribution (the Voigt profile) by using a combination of analyticalmethods and numerical integration.4.1THE WORKING FUNCTION AND ITS PROPERTIESFor purposes of comparing various methods, particularly those for numerical derivatives and integrals, it is useful to have a working function and to systematize its analytical properties.
A polynomial, which can readily be differentiated and integratedanalytically to provide benchmark values, is more appropriate than the functions thatwe investigated by power series in Section 3.2. Further, most of our approximations are based on the use of low-order polynomials, so by using such a polynomialwe can make direct analytical comparison between approximate and exact values.The working function will stress the methods greatest if the polynomial oscillateswithin the region where the derivatives and integrals are being performed, and suchoscillations can be designed into the polynomial by choosing its roots appropriately.In this section we first derive some analytical and numerical properties of theworking function, then we develop a C function for efficient numerical evaluation ofpolynomials by using Horner’s algorithm, then we include this in a program for theworking function.Properties of the working functionWe define the working function, yw(x), by a product expansion in terms of itsroots, defining it by the expression(4.1)which is a sixth-order polynomial having all real roots, and five of these are withinthe range -0.5O.5 over which we will explore the function properties.
Thisworking function, yw is shown in Figure 4.1 for a slightly smaller range than this.Outside the range of the graph, the function has another zero at x = -0.5 and one atx = 1.0.Although the product expansion (4.1) is convenient for controlling the positionsof roots of the working function, and thereby controlling its oscillations, it is notconvenient for analytical differentiation or integration. These can be performed byexpanding (4.1) into powers of x, then differentiating or integrating term by term.This involves a little algebra and some numerical work, so why don’t you try it?4.1 THE WORKING FUNCTION AND ITS PROPERTIES101FIGURE 4.1 The working function (4.1), a pathological polynomial having six real roots inthe interval [-0.5, 1].Exercise 4.1(a) Show that the polynomial for the working function, (4.1), can be expandedinto(4.2)where the coefficients, wk, are given in Table 4.1,(b) Starting with the expansion (4.2) with the given wk, evaluate the successivederivatives of yW in terms of the expansion coefficients appearing in the formulas(4.3)in which for the i th derivative, i = 1, 2, .
. . . 6, the coefficients are given inTable 4.1.(c) Perform a similar derivation for the integral of the working function yw ( x )in order to show that its indefinite integral is given by(4.4)for which the coefficients Ik are given in the bottom row of Table 4.1. n102NUMERICAL DERIVATIVES AND INTEGRALSTABLE 4.1 Expansion coefficients of working function (4.1) for the value, wk, for derivatives,wi,k, and for the indefinite integral, Ik. Omitted values are zero.0123467066001This table, transformed into an array, is used in the program at the end of this section. It may also be used for the programs in Sections 4.4, 4.5, and 4.6.Derivatives of the working function, appropriately scaled, are shown in Figure 4.2.
Notice that the derivatives grow steadily in magnitude, but steadily becomesmoother with respect to x as the order of the derivative increases. Clearly, for ourpolynomial of sixth order, derivatives past the sixth are zero, as the constancy ofy(6) in Figure 4.2 indicates.-0.50FIGURE 4.2 The working polynomial, (4.1) or (4.2), and its derivatives. Note that the derivatives become larger but smoother as the order of the derivative increases.It is useful to have available a program for computing the working function (orother polynomial) and its derivatives. Since this function, its derivatives, and its integral are all polynomials, it is practical to have an efficient means of evaluatingthem.
Horner’s polynomial algorithm, which we now describe and for which weprovide a C function for its implementation, is efficient, accurate, and simple.4.1 THE WORKING FUNCTION AND ITS PROPERTIES103A C function for Horner’s algorithmPolynomials occur very frequently in numerical mathematics, for example in anyTaylor expansion approximation to a function, as considered in Section 3.2. It istherefore worthwhile to find an accurate and efficient way to evaluate polynomialsnumerically.
Horner’s method is one approach. Suppose that the polynomial to beevaluated, y ( x ), is of order N in the variable x , so that(4.5)We can avoid the inefficient (and sometimes inaccurate) direct evaluation of powersof x by starting with the highest-order term, then adding in the next-lowest coefficient plus the current value of the polynomial multiplied by x. For example, supposethat we have the quadratic polynomial(4.6)We may evaluate this by successive stages as(4.7)(4.8)(4.9)In each of these stages the symbolmeans “assign the right side to the left side.”After the last assignment, (4.9), we have completely evaluated the polynomial (4.5)without computing powers of x explicitly. Clearly, this recursive procedure can begeneralized to any order of polynomial, N. The Homer method provides an example of the effective use of repetition in computing, as we discuss in the Diversion inSection 3.4.Program 4.1, Horner Polynomials, implements the Horner algorithm for apolynomial of order N, where N is restricted only by the size of the array, a[MAX ] ,used to hold the values of the ak, k = 0, 1 , .
. . . N. In this program the constantterm in the polynomial, a0, is stored in a[0] , and so on for the other coefficients.Thus, we have an exception to our practice, stated in Section 1.2, of starting arraysat element [1]. By breaking this rule, however, we have a direct correspondencebetween the formula (4.5) and the coding.The main program consists mainly of bookkeeping to input and test the order ofthe polynomial N, to test its bounds (a negative value produces a graceful exit fromthe program), then to input the polynomial coefficients. Once these are chosen bythe program user, the polynomial, y, can be evaluated for as many x values as desired, except that input of x = 0 is used to signal the end of this polynomial.104NUMERICAL DERIVATIVES AND INTEGRALSPROGRAM 4.1 A test program and function for evaluating an Nth-order polynomial by theHorner algorithm.#include <stdio.h>#include <math.h>#define MAX 10/* Homer Polynomals */double a[MAX];double x,y;int N,k;double Horner_Poly();printf("Horner Polynomials\n");N = 1;while(N>=O){printf("\nNew polynomial; Input N (N<0 to end):");scanf("%i",&N);if(N<O){ printf("\nEnd Homer Polynomials"); exit(O); }if(N>MAX-1)printf("\n\n!! N=%i > maximum N=%i",N,MAX-1);}elseprintf("\nInput a[k], k=O,...,%i (%i values)\n",N,N+l);for ( k = 0; k <= N; k++ )scanf("%lf",&a[k]);x = 1;while(x!=0){printf("\n\nInput x (x=0 to end this polynomial):");scanf("%lf",&x);if ( x .!= 0 ){y = Horner_Poly(a,N,x);printf("y(%lf) = %lf",x,y);}} /*end x loop*/} /*end N loop*/4.1 THE WORKING FUNCTION AND ITS PROPERTIES105double Horner_Poly(a,N,x)/* Horner Polynomial Expansion */double a[],x;int N;double poly;int k;for ( k = N-l; k >= 0; k-- )poly = a[k]+poly*x;return poly;The Horner polynomial is in function Horner-Poly, which is very compact.Note that the for loop is not executed if N = 0, because the value of k is checked before any instructions within the loop is executed.
The function is therefore correct even for a zero-order polynomial.Exercise 4.2(a) Check the correctness of the Horner polynomial algorithm as coded in function Horner-Poly in Program 4.1 by using it to check the case N = 2 workedin formulas (4.7) through (4.9).(b) code program Horner Polynomials, then test it for a range of N values,x values, and correctness for some low-order polynomials. nNote that in the Horner algorithm the order of the polynomial, N, must beknown before the polynomial is evaluated.
The algorithm is therefore of limited usefor computing power-series expansions (as in Chapter 3), unless one has decidedbeforehand at what order to terminate the expansion.Program 3.3 in Section 3.2 gives an example of this latter usage for the functions CosPoly and SinPoly. On the other hand, if a convergence criterion (suchas a given fractional change in the series partial sum) is used, then the method ofgenerating successive terms by recurrence (Sections 3.4, 3.5) is appropriate. Program 3.4 in Project 3 (Section 3.5) shows this usage in the functions for the sixpower series.The Horner-poly function should be of general use to you in future programdevelopment. We apply it immediately to program the numerical evaluation of theworking function presented above.106NUMERICAL DERIVATIVES AND INTEGRALSProgramming the working functionThe working function, described in Section 4.1, is used in Sections 4.4 and 4.5 totest numerical derivatives and may be used in Section 4.6 for integrals.
It would bepossible, but tedious (especially if you use a different polynomial), to evaluate thenumerical values by hand or with help from a pocket calculator. Also, if one wantsthe derivatives of various orders or the integral this is also tedious and error-prone inhand calculations. Here we provide a general-purpose program that, given the orderand coefficients of a polynomial, determines the polynomials that represent its nonvanishing derivatives. Then the program allows the user to choose specific x valuesfor evaluating these quantities numerically.Here we give a preliminary discussion of the program Working Function:Value and Derivatives (henceforth abbreviated to Working Function).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
