Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 58
Текст из файла (страница 58)
Westart out with the composite trapezoid rule approximation (see Sec. 7.4)(7.72)to the numberHere N is a positive integer related to h byandf i = f i , N = f(a + ih)i = 0, . . . , NIf f(x) is four times continuously differentiable, we infer from (7.50) and(7.51) thatI(f) = T N (f) + C 1 h 2 +(7.73)where the constant C 1 = [f´(a) - f´(b)]/12 is independent of h. Henceextrapolation to the limit is applicable. We get thatis anapproximation to I(f), while in general, TN(f) has only anerror ofNote that the choice of q or N is restricted by the condition that N/qbe an integer. One usually chooses q = 2 (so that N must be even).
Thischoice for q has the computationally important advantage that all functionvalues used for the calculation of TN/q can also be used for the calculationof TN. Specifically, we prove that for even N,(7.74)For by (7.72),Here the first sum extends over the “odd” points and the second sum over*7.7ROMBERG INTEGRATION341the “even” points. The last two terms can be writtenHence, sincethese last two terms add up to TN/2 (f)/2. This proves (7.74). Note that(7.74) can be written more simplywith M denoting the composite midpoint rule (7.49a ) .If the integrand has 2k + 2 continuous derivatives, it can be shownthat, more explicitly than (7.73),where the constants C1, .
. . , Ck do not depend on h. Hence, withwe get thatwith the constants C21, . . . , Ck1 independent of h. Further extrapolation istherefore meaningful. Settingwe get thatMore generally, it is seen that, for m = 1, . . . , k,is anapproximation to I(f).Note that the calculation of T N m involvesand finally,must therefore be an integer,henceand T N . N/2msay, for TNm to be defined.
It is convenient to visualize these various342DIFFERENTIATION AND INTEGRATIONapproximations to I(f) as entries of a triangular array, the so-called Ttable:Here we have written TN0 for TN.Algorithm 7.2: Romberg integration Given a function f(x) defined on[a,b] and a positive integer M (usually, M = 1).h := (b - a) /MIf f(x) has 2m + 2 continuous derivatives, thenk = m, m + l, · · ·Also, if k is sufficiently large, thenBut before putting any faith in this inequality, check that at leastExample 7.9 Use Romberg integration for Example 7.1.The integral in question isThe FORTRAN program below has been set up to produce the first six rows of theT table and the corresponding table of ratios Rkm, as follows:Romberg T table0.7313700E0.7429838E0.7458653E0.7465842E0.7467639E0.7468069E0000 0.7468551E 0000 0.7468258E 00 0.7468238E 0000 0.7468238E 00 0.7468237E 0000 0.7468237E 00 0.7468237E 0000 0.7468212E 00 0.7468210E 000.7468237E 000.7468237E 00 0.7468237E 000.7468210E 00 0.7468209E 00*7.7ROMBERG INTEGRATION343Table of ratios4.034.014.004.1714.8816.500.050.00.00.0M was chosen to be 2, so that the first entry in the T table is T2(f).
Note that thefirst column of ratios converges very nicely to 4, but then begins to move away from 4.This effect is even more pronounced in the second column of ratios, which approach 16(as they should), and then, as the last entry shows, become erratic. Conclusion: Theerror in the entries of the last row of the T table is mainly due to roundoff (rather thandiscretization). Hence0.7468237seems to be the best estimate for I(f) to be gotten with the particular arithmetic used.Sinceandto the number of places shown, we conclude that this estimate is accurate to the numberof places shown. Actually,The discrepancy between this number and our “accurate” estimate is due to the fact thatwe are not dealing with the integrandin our calculations, but rather with a rounded version of f(x), that is, with the functionF(X) = EXP(-X*X)All calculations were carried out in single precision on an IBM 360, which hasparticularly poor rounding characteristics.FORTRAN PROGRAM FOR EXAMPLE 7.9REAL T(l00)EXTERNAL FERRCALL RMBERG( FERR, 0., l., 2, T, 6 )STOPENDSUBROUTINE RMBERG ( F, A, B, MSTART, T, NROW )C CONSTRUCTS AND PRINTS OUT THE FIRST NROW ROWS OF THE ROMBERG TTABLEFOR THE INTEGRAL OF F(X) FROM A TO B , STARTING WITH THECC TRAPEZOIDAL SUM ON MSTART INTERVALS.INTEGER MSTART,NROW, I,K,MREAL A,B,T(NROW,NROW), H,SUMM = MSTARTH = (B-A)/MSUM = (F(A) + F(B))/2.IF (M .GT.
1) THENDO 10 I=l,M-1SUM = SUM + F(A+FLOAT(I)*H)10END IFT(l,l) = SUM*HPRINT 610344DIFFERENTIATION AND INTEGRATION610 FORMAT('l', l0X,'ROMBERG T-TABLE'//)PRINT 611, T(l,l)611 FORMAT(7E15.7)RETURNIF (NROW .LT. 2)C11C1220C6202530630DO 20 K=2,NROWH = H/2.M = M*2SUM = 0.DO 11 I=l,M,2SUM = SUM + F(A+FLOAT(I)*H)T(K, 1) = T(K-1,1)/2. + SUM*HDO 12 J=l,K-1SAVE DIFFERENCES FOR LATER CALC.
OF RATIOST(K-l,J) = T(K(,J) - T(K-1,J)T(K,J+l) = T(K,J) + T(K-l,J)/(4.**J - 1.)PRINT 611, (T(K,J),J=l,K)RETURNIF (NROW .LT. 3)CALCULATE RATIOSPRINT 629FORMAT(///llX,'TABLE OF RATIOS'//)DO 30 K=l,NROW-2DO 25 J=l,KIF (T(K+l,J) .EQ. 0.) THENRATIO= 0.ELSERATIO = T(K,J)/T(K+l,J)END IFT(K,J) = RATIOPRINT 630, (T(K,J),J=l,K)FORMAT(8Fl0.2)RETURNENDREAL FUNCTION FERR(X)REAL XFEAR = EXP(-X*X)RETURNENDEXERCISES7.7-1 Prove that, in Romberg integration,rule; see (7.48).7.7-2 Try to estimate I(f) =of the following cases:(a) f(x) = x 2(b) f(x) = sin 101πx(c) f(x) = 1 + sin 10π xSM, where SM is the composite Simpson’sto within 10-6, using Romberg integration, for each(d) f(x) = |x - 1/3 |aaaa====0, b = 1, M arbitrary0, b = l, M - 10, b = l, M - 10, b = l, M = 1 and M = 3(e) f(x) =a = 0, b = 1, M arbitrary7.7-3 From the data below calculateas accurately as possible using Rombergintegration.
Construct a T table starting with M = 1*7.7x1.01.11.21.31.41.51.61.71.8ROMBERG INTEGRATION345f(x)0.367879440.366158190.361433050.354291330.345235740.334695240.323034430.310561990.297538007.7-4 Obtain Simpson’s rule for Ih(f) =by extrapolating from the midpoint ruleand the trapezoid rule. (Hint: Form the appropriate linear combination of the two equationsI h (f) = T(f) + C T h 2 +Ih(f) = M(f) + CMh2 +to eliminate the h2 terms.
This requires you to find out what the constants CT and CM are.Previous Home NextCHAPTEREIGHTTHE SOLUTION OF DIFFERENTIALEQUATIONSMany problems in engineering and science can be formulated in terms ofdifferential equations. A large part of the motivation for building the earlycomputers came from the need to compute ballistic trajectories accuratelyand quickly. Today computers are used extensively to solve the equationsof ballistic-missile and artificial-satellite theory, as well as those of electrical networks, bending of beams, stability of aircraft, vibration theory, andothers.It is assumed that the student is familiar with the elementary theory ofdifferential equations.
In a first course one learns various techniques forsolving in closed form some selected classes of differential equations. Thevast majority of equations encountered in practice cannot, however, besolved analytically, and recourse must necessarily be made to numericalmethods. Fortunately, there are many good methods available for solvingdifferential equations on computers. In this chapter we shall derive severalclasses of methods, and we shall evaluate them for computationalefficiency.8.1 MATHEMATICAL PRELIMINARIESIt will be useful to review some elementary definitions and concepts fromthe theory of differential equations.
An equation involving a relationbetween the values of an unknown function and one or more of itsderivatives is called a differential equation. We shall always assume that3468.1MATHEMATICAL PRELIMINARIES347the equation can be solved explicitly for the derivative of highest order. Anordinary differential equation of order n will then have the formy (n) (x) = f(x, y(x) ,y´(x), .
. . ,y (n-1) (x))(8.1)By a solution of (8.1) we mean a functionwhich is n times continuously differentiable on a prescribed interval and which satisfies (8.1); thatismust satisfyThe general solution of (8.1) will normally contain n arbitrary constants,and hence there exists an n-parameter family of solutions. If y(x 0 ),y´(x0), . . .
, y(n-1)(x0) are prescribed at one point x = x0, we have aninitial-value problem. We shall always assume that the function f satisfiesconditions sufficient to guarantee a unique solution to this initial-valueproblem. A simple example of a first-order equation is y´ = y. Its generalsolution is y(x) = Cex, where C is an arbitrary constant. If the initialcondition y(x0 ) = y0 is prescribed, the solution can be written y(x) =Differential equations are further classified as linear and nonlinear. Anequation is said to be linear if the function f in (8.1) involves y and itsderivatives linearly.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
