Thompson - Computing for Scientists and Engineers (523188), страница 26
Текст из файла (страница 26)
Include program statements to compute the difference between these derivativesand those obtained from the five formulas for first and second derivatives.thatare already programmed. Test this new program, called Cosine FunctionNumerical Derivatives, for example against hand calculations.4.6 NUMERICAL INTEGRATION METHODS133(b) Run the program with x = 0 and successively smaller values of stepsize h,say h = 0.2, 0.1, 0.01, 0.001, calculating all five numerical derivatives andtheir errors for each h. Except for the central-difference first derivative (whichhas zero error at x = 0), make a table and a log-log plot of the absolute value ofthe error against h, as in Table 4.8 and Figure 4.7.(c) Assuming that the graph for each method on a log-log plot is a straight line,show that the error has a power-law dependence on stepsize h, and that the overall scale of the error estimates the factor preceding the power of h.
How well dothe powers estimated for the five derivative methods from your graph agree withthose predicted from our Taylor-series analyses? nTABLE 4.8 Central-difference first and second derivatives errors for the cosine function evaluatedat x = 0. The exact second derivative is -1.Notice in Figure 4.7 that the most accurate numerical method, the five-point centraldifference second derivative, does not continue to improve as rapidly as predictedwhen h becomes of order 10-3.
For this small an h value the error in the numericalapproximation approaches my computer’s roundoff error of about 10-14.Now that you have experience in numerical differentiating with the badly behaved working function, with a monotonic exponential, and with the oscillatory cosine,you should be confident to tackle many functions and the differentiation of numericaldata.There are several other practical methods for numerical estimating derivatives.For example, cubic-spline fitting through a set of points will also produce, as part ofthe process, first through third derivatives, as described in Section 5.4.
Anothervery effective method is to consider the derivatives as functions of the stepsize h,then to extrapolate the derivative values at finite h to the limit h 0. This method,called Richardson extrapolation, is described in Section 7.2 of Maron’s book.4.6 NUMERICAL INTEGRATION METHODSWe remarked at the beginning of this chapter that integrals of functions are muchless likely to be able to be calculated analytically than are derivatives, but that integration is usually a more stable problem (in the sense discussed in Section 4.3) than134NUMERICAL DERIVATIVES AND INTEGRALSis differentiation. We therefore seek efficient numerical-integration methods becausewe will probably have to do many such integrals (absent analytical results), but afairly simple numerical algorithm may perform surprisingly well.The differentiation-integration question can be approached from another viewpoint, namely that the integral of a function is formally the “antiderivative” of thatfunction.
Since we saw in Section 4.4 that successively higher-order derivativesbecome numerically more difficult, therefore (at least formally) successive integralsshould become numerically less difficult. As we discussed for derivatives at the beginning of Section 4.4, we also base our integration formulas on Taylor expansionsof the integrand. If we truncate the expansion to a given order, N, we are using anNth-order polynomial approximation. The first neglected term in the Taylor expansion serves to estimate the error in the polynomial approximation.The working function, introduced at the beginning of this chapter as (4.1), isused for tests and for comparisons of the methods we develop.
Recall that this function was designed as a highly oscillatory but bounded function near x = 0, asshown in Figures 4.1 and 4.2. The expansion of its integral is given by (4.4), using the coefficients Ik in Table 4.1. The coefficients of the derivatives of variousorders, also in this table, are useful for estimating errors in integral approximations.We derive the two most common integration formulas, the trapezoid and Simpson approximations.
The emphasis in our discussion is on understanding how toderive these formulas and how to estimate their errors, with a view to indicating howto generalize the procedure for such approximations. Our approach is also consistent with that used in Section 4.4 for obtaining the derivative estimates.Integration algorithms using constant stepsize h are readily derived by using aTaylor expansion of the integrand within the region to be integrated.
Suppose thatthe integration is centered on x0, so that the appropriate Taylor expansion of an integrand, y (x), through the nth power of (x - x0) can be obtained by using Taylor’stheorem, (3.6), as(4.39)in which the derivative is evaluated at x0 and the remainder, Rn, is formally given by(3.7), but it may be approximated by the next term in the power series, which iswhat we shall use. The indefinite integral over x of y (x) is obtained by term-byterm integration of the series (4.39), a procedure which is usually valid for integralsof interest.Exercise 4.19Show that integration of the series terms in (4.39) leads to(4.40)for the indefinite integral. n4.6 NUMERICAL INTEGRATION METHODS135The presence of n derivatives in this formula, ignoring the remainder integral, requires that y be known at (n + 1) values for each basic integral.
We now show theconstruction explicitly for n = 1 (trapezoid rule) and for n = 2 (Simpson rule).Trapezoid formula and program for integrationTrapezoid rule formula. The lowest-order formula obtained using (4.40) is forn = 1, but let us carry the expressions to n = 2 then neglect the integral over the remainder. The n = 2 term will then be an estimate of the error in the n = 1 approximation. Suppose that we wish to integrate from x = a to x = a + h. The mostrapid convergence of the series (4.40) will then probably be achieved by choosingx 0 =a + h/2.Exercise 4.20(a) Make the indicated substitutions of variables in (4.40) in order to show that(4.41)in which derivatives of fourth and higher order have been neglected.(b) Write out the Taylor expansions (4.39) of y (a) and y (a + h) aboutx0 = a + h /2, then solve for y (a + h /2) in terms of these y values and thesecond derivative at a + h/2.
Thus show that (4.41) becomes(4.42)which, with neglect of the derivative term, is the basic trapezoid formula.(c) Show that if y is linear in x, then the trapezoid rule is exact. nThis derivation might seem a bit long-winded, because all we have ended up doing is to average the y values at the beginning and end of the integration interval,then multiplying by the length of the interval, which is h. However, if we did notuse the more complete analysis that you have just made (haven’t you?), then wewould not be able to estimate the error in the approximation.Notice that this error is better than you might have guessed, since the effects ofthe first derivative at the midpoint of the interval has vanished, leaving the secondderivative as the leading term in the error estimate.
By using the midpoint of the integration interval as the expansion point, all odd-order derivatives vanish from theintegration formula. Therefore, the next neglected term in the integral, and thereforein our error estimate, involves y( 4 )h 5 .136NUMERICAL DERIVATIVES AND INTEGRALSFIGURE 4.8 Integration formulas, illustrated for the working function (4.1) with stepsizeh = 0.1 in [-O.l, 0.1]. The trapezoid rule uses a straight-line segment in each panel. The Simpson rule uses the parabolic section (dashed) over two panels.Graphical interpretation of the trapezoid rule clarifies the analysis.
Figure 4.8shows our working function, (4.1), near x = 0, the same region where we investigate derivatives in Sections 4.4 and 4.5. The geometric interpretation of Exercise 4.20 (c) is that we approximate the function between the tabulated values atx = a and x = a + h by a straight-line segment. For example, in Figure 4.8 wehave a trapezoid between x = -0.1 and x = 0, with another between x = 0 andx = 0.1.Starting from the basic trapezoid formula (4.42), we may use additivity of integration to find an approximation to the integral between x = a and x = a + nh.Thereby we produce the composite trapezoid formula(4.43)in which the error estimate is(4.44)Here the second derivative is not larger in magnitude than the maximum second derivative in the entire interval.
The error scales as the quantity ny(2) h 3 /y, so halving the stepsize h (while doubling n to cover the same x range) should produceabout a factor of 4 greater accuracy in the integral estimated by the trapezoid rule.Trapezoid rule program. Programming the trapezoid function is quite straightforward. As usual, the code for the driver program to make the function testable islonger for the code than for the function itself. Also, we anticipate including theSimpson rule in the next subsection, so some code is included for this. Here weemphasize parts relevant to the trapezoid rule.4.6 NUMERICAL, INTEGRATION METHODS137PROGRAM 4.5 Trapezoid and Simpson integrals, composite formulas, applied to the workingfunction (4.1).#include#include<stdio.h><Math.h>main(){/* Trapezoid and Simpson Integrals */double a,h,IntAnalyt,IntByTrap,error,IntBySimp;int n,NTrap,NSimp;doubleSimpInt(),TrapInt(),y(),YWInt();printf("Trapezoid and Simpson Integrals\n");n = 2;while(n!=O)printf("\n\nInput n, a, h (n=O to end): ");scanf("%i%lf%lf",&n,&a,&h);if ( n == 0 )printf("\nEnd Trapezoid and Simpson Integrals");exit (0) ;}NTrap = n;if ( n == 2*(n/2) )NSimp = n;else NSimp = O;/* to bypass Simpson for n odd */IntAnalyt = YWInt(a+n*h)-YWInt(a);printf("Analytica1 integral=%lg",IntAnalyt);InByTrap = TrapInt(a,NTrap,h);/* trapezoid rule */error = IntAnalyt-IntByTrap;printf("\nTrapezoid integral=%lg",IntByTrap);printf(" Error=%lg",error);if ( NSimp != 0 )IntBySimp = SimpInt(a,NSimp,h);/* Simpson rule */error = IntAnalyt-IntBySimp;printf("\nSimpsonintegral=%lg",IntBySimp);printf(" Error=%lg",error);} /*end n loop */138NUMERICAL DERIVATIVES AND INTEGRALSdouble SimpInt(a,n,h)/* Simpson integral */double a,h;int n;double sum;double y();int n2,j;n2 = n/2;if ( (2*n2 != n) || (n2 == 1) ){printf("\n !! SimpInt has n=%i; not allowed",n);return 0;}elsesum= 0;for ( j = 1; j < n2; j++ ){sum = sum+y(a+2*j*h)+2*y(a+(2*j-l)*h);return h*(y(a)+4*y(a+(n-l)*h)+y(a+n*h)+2*sum)/3;double TrapInt(a,n,h)/* Trapezoid integral */double a,h;int n;{double sum;double y();int j;sum = ( y(a) + y(a+n*h) )/2;for ( j = 1; j < n; j++ ){sum = sum+y(a+j*h);}return h*sum;4.6 NUMERICAL INTEGRATION METHODS139double y(x)/* Working function value*/double x;double prod;prod=12O*(x+0.5)*(x+0.25)*x*(x-(l.0/3.0))*(x-0.2)*(x-1);return prod;double YWInt (x)/* Working function integral */double x;double sum;sum=(-0.5+x+(23.0/4.0)*x*x-10.2*pow(x,3)-(47.0/3.0)*pow(x,4)+(120.0/7.0)*pow(x,5))*x*x;return sum;The structure of Trapezoid and Simpson Integrals can be summarized asfollows.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
