Thompson - Computing for Scientists and Engineers (523188), страница 21
Текст из файла (страница 21)
Wedefer discussion of the functions for numerical derivatives (function names startingwith Der) to Sections 4.4 and 4.5. Working Function is structured as follows.The main program outermost loop is over the order of the polynomial, Npoly.Within this loop are input and testing sections for Npoly, a section for the input ofthe polynomial coefficients, a section where derivative coefficients are calculated,then a loop for choosing x values at which to evaluate the function and its first andsecond derivatives.The polynomials are evaluated using the Horner algorithm and a modified version of the C function from the preceding subsection. Function Horner-Poly-2has a two-dimensional array (rather than one-dimensional, as in Horner-Poly), andan additional argument, i, to control the first index of the array.
For the functionand its derivatives the array elements are the wk in (4.2) as the zeroth elements of thearray, and the Wi,k in (4.3) as the elements for the ith derivative. At the expense of aslight awkwardness, this array is also two-dimensional, which allows the Hornerpolynomial program function to be used also to evaluate coefficients for the integral.Note that the labeling of arrays from zero, rather than unity, in the C language isa great advantage in Working Function.
Therefore, to maintain a closeness between the analysis and the numerics, we break with our general rule (Section 1.2)of starting arrays at unity (as is common practice among Fortran programmers).Exercise 4.3Code and test Working Function. For testing it may be useful to add an output section for the coefficient elements w[i][k] and Int eg [0][k] .
Test coefficients for all these are given in Table 4.1 for the working function (4.2),which is equivalent to (4.1). The latter provides x values at which the function,but not its derivatives, should evaluate to zero. nWe now have a working function and a program for its numerical properties touse in examples and tests of the numerical methods in this chapter. After listing theprogram, we consider some limitations of numerical computation.4.1 THE WORKING FUNCTION AND ITS PROPERTIESPROGRAM 4.2 The working function (4.1) or (4.2): its values and derivatives.#include <stdio.h>#include <math.h>#define MAXDER 7main (){/* Working Function: Value and Derivatives */double w[MAXDER][MAXDER];double x,ywDeriv_i,h;double ym2.ym1.y0,yp1,yp2,yf,yC,yCCD,y3CD,y5CD;int Npoly,k,i;char polyOK,choose_x;doubleHorner_Poly_2(),Der_1F(),De_r1C(),Der_2CCD(),Der_23CD(),Der25CD();printf(Working Function: Value and Derivatives\n");Npoly = 1;while ( Npoly >= 0 ){printf("\nNew polynomial;Input order,Npoly (<0 to end): ");scanf("%i",&Npoly);if(Npoly<O){printf("\nEnd Working Function Value and Derivatives");exit (0);}polyOK = 'y';if ( Npoly >= MAXDER ){printf("\n !!Order=%i>max_order=%i\n",Npoly,MAXDER-1);polyOK = 'n';elseprintf("\nInput %i coefficients:\n", Npoly+l);for ( k = 0; k <= Npoly; k++ )scanf("%lf",&w[O][k]);107108NUMERICAL DERIVATIVES AND INTEGRALS/* Generate all derivative coefficients */if (Npoly > O)for ( i = 1; i <= Npoly; i++ ) /* derivative order */{for ( k = 0; k <= Npoly-i; k++ )/* polynomial terms */{W [i ] [k] = (k+l)*w[i-1][k+l];}}}}if ( polyOK == 'y' )/* Evaluate working function at input x values */choose_x = 'y';while ( choose_x == 'y' )printf("\n\nWant an x ( y or n)? ");scanf("%s",&choose_x);if ( choose_x == 'y' )printf("Input x and stepsize h: ");scanf("%lf%lf",&x,&h);for ( i = 0; i <= Npoly; i++ )ywDeriv_i = Horner_Poly_2(w,Npoly-i,i,x);if ( i == 0 )printf("\nPolynomial value is %lf",ywDeriv_i);elseprintf("\ni=%i, derivative is %lf",i,ywDeriv_i);/* Numerical derivatives; function values first */ym2 = Horner_Poly_2(w,Npoly,0,x-2*h);ym1 = Horner_Poly_2(w,Npoly,0,x-h);y0 = Horner_Poly_2(w,Npoly,0,x);yp1 = Horner_Poly_2(w,Npoly,0,x+h);yp2 = Horner_Poly_2 (w,Npoly,0,x+2*h);/* First derivatives: */yF = Der_1F(y0,yp1,h); /* Forward *//* Central */yC = Der_1C(ym1,yp1,h);/* Second derivatives: */yCCD = Der_2CCD(ym2,y0,yp2,h);/* Central */y3CD = Der_23CD(ym1,y0,yp1,h);/* 3-point */y5CD = Der_25CD(ym2,ym1,y0,yp1,yp2,h);/* 5-point */4.1 THE WORKING FUNCTION AND ITS PROPERTIESprintf("\nNumerical first derivatives:");printf("yF = %lf, yC = %lf",yF,yC);printf("\nNumerical second derivatives:");printf("\nyCCD=%lf y3CD=%lf y5CD=%lf",yCCD,y3CD,y5CD);}} /* end while loop for x */}} /* end Npoly loop */}double Horner_Poly_2 (a,N,i,x)/* Homer Polynomial Expansion *//* Coefficients have two subscripts;first subscript is i */double a[MAXDER][MAXDER],x;int N,i;{double poly;int k;poly = a[i][N];for ( k = N-l; k >= 0; k-- ){poly = a[i] [k]+poly*x;}return poly;}double Der_1F(y0,yp1,h)/* Forward-difference first derivative */double y0,yp1,h;{double y1F0;y1F0 = (yp1-y0)/h;return y1F0;}double Der_1C(ym1,yp1,h)/* Central-difference first derivative */double ym1,yp1,h;{double y1C0;y1C0 = (yp1-ym1)/(2*h);return y1C0;}doubleDer_2CCD(ym2,y0,yp2,h)109110NUMERICAL, DERIVATIVES AND INTEGRALS/* Central-difference first derivative */double ym2,y0,yp2,h;double y2CCD;y2CCD = ( (yp2-y0)-(y0-ym2) ) / (4*h*h);return y2CCD;}double Der_23CD(ym1,y0,yp1,h)/* Three-point central-difference derivative */double ym1,y0,yp1,h;{double y23CD;y23CD = ( (yp1-y0)+(ym1-y0) ) / (h*h);return y23CD;}doubleDer_25CD(ym2,ym1,y0,yp1,yp2,h)/* Five-point central-difference derivative */doubleym2,ym1,y0,yp2,h;{double diff1,diff2,y25CD;diffl = (yp1-y0) + (ym1-y0);diff2 = (yp2-y0) + (ym2-y0);y25CD = ( 4*diff1 - diff2/4 )/(3*h*h);return y25CD;}4.2DISCRETE DATA AND NUMERICAL MATHEMATICSWe wish to emphasize briefly the importance of the discreteness of data and of numerical mathematics.
Both topics place fundamental limits on the completeness ofour interpretation of quantitative information.The discreteness of dataIn experimental work and also in numerical calculations, data are always obtained atdiscrete values of the independent variables, even though we may believe that thequantity being measured is a smooth function of the parameters (independent variables) of the problem.
Similarly, in numerical applications of mathematical analysisthe values are always obtained at finite parameter spacing. Such values are often,quite appropriately, called “numerical data.” Subsequent operations on data, such asapproximations to derivatives and integrals of the presumed smooth function underlying them, will always be limited in accuracy by the discreteness of the data.4.3NUMERICAL NOISE IN COMPUTING111Numerical mathematicsIn this chapter, and in much of this book, we have a different emphasis than in puremathematics. The latter emphasizes the existence of mathematical entities, the logicalconsistency of results, and the analytic solution and properties of quantities such asderivatives and integrals.
Numerical analysis, such as we are developing, must always be guided by pure mathematics, especially on questions of the existence andconvergence of quantities that we attempt to evaluate numerically. This viewpointexplains why in Figure 1.1 we showed mathematical analysis as the foundationupon which are built numerical methods and computer programming for various scientific applications.When evaluating a formula numerically, discrete values of the variables are always required, rather than the continuous range of values often allowed by the formula definition. For practical reasons, such as minimizing computer storage and execution time, one often characterizes a function by calculating accurately as few values as possible with large stepsizes for the independent variables.
Related propertiesof the function, such as values at intermediate values of the variables (interpolation),slopes (differentiation), and area under the function between two limits (integration),are then estimated by numerical procedures.One often refers to mathematics that emphasizes numerical properties as numerical analysis, and the process of replacing continuous functions by discrete values iscalled discretization. They are to be distinguished from two other distinct branchesof mathematics, number theory and discrete mathematics.4.3 NUMERICAL NOISE IN COMPUTINGIn a digital computer, most numbers can be represented to only a finite number ofsignificant figures, unless they can be expressed exactly in binary (base-2) arithmeticrequiring fewer than the number of bits used for each number.
For example, in a32-bit computer word the largest integer is 232- 1 (allowing for a sign bit), andthere can be at most a precision of 1 part in 232, that is, less than 1 part in 1010. Thedecimal fraction 0.5 is exactly represented, since 0.5 = 2-1, so that in base-2 arithmetic 0.510 = 0.12. On the other hand, however, the precise decimal fraction0.410 = 0.0110011001...2.
which is an imprecise binary fraction.The term “computer arithmetic,” in which the number of digits is finite, is oftenused to make a distinction from formal or exact arithmetic. Note that errors in numerical calculations introduced by computer arithmetic must be clearly distinguishedfrom random errors introduced in computations by computer hardware (or software)malfunctions. The type of error that we discuss should be completely reproducibleover time, provided that you execute exactly the same program as previously. (Inthe worst possible case of program bugs, this may be difficult to do unless the computer is performing exactly the same from run to run.) Note that the random errorrate in the arithmetic hardware of a computer is typically less than 1 in 1012, whichis less than one error per year of continuous running.112NUMERICAL DERIVATIVES AND INTEGRALSTwo main methods of handling the least-significant digits after an arithmetic operation are roundoff and truncation.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
