Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 61
Текст из файла (страница 61)
. . Based on this, write a FORTRAN program which findsthe value y(3) of the solution y(x) of the problem in Example 8.1 to within 10-6, usingAlgorithm 8.1.8.4 ERROR ESTIMATES AND CONVERGENCE OF EULER’SMETHODTo solve the differential equation y´ = f(x,y), y(x o ) = y 0 by Euler’smethod, we choose a constant step size h, and we apply the formulayn+1 = yn + hf(xn,yn)n = 0, l, . . .(8.23)where xn = x0 + nh. We denote the true solution of the differential equation at x = xn by y(xn ), and the approximate solution obtained by applying (8.23) as yn. We wish to estimate the magnitude of the discretizationerror en, defined by(8.24)e n = y(x n ) - y nWe note that, if y0 is exact, as we shall assume, then e0 = 0. Assuming thatthe appropriate derivatives exist, we can expand y(xn+1 ) about x = xn ,using Taylor’s theorem with remainder:xn < ξn < xn+1 (8.25)The quantity (h 2 /2)y´´( ξ n ) is called the local discretization error, i.e., theerror committed in the single step from xn, to xn+1, assuming that y and y´were known exactly at x = xn.
On a computer there will also be an error incomputing yn+1 , using (8.23), due to roundoff. Round-off errors will beneglected in this section.360THE SOLUTION OF DIFFERENTIAL EQUATIONSOn subtracting (8.23) from (8.25) and using (8.24), we obtain(8.26)By the mean-value theorem of differential calculus, we havewhereis between yn and y(xn). Hence (8.26) becomes(8.27)We now assume that over the interval of interest,|f y (x,y)| < L|y´´(x)| < Ywhere L and Y are fixed positive constants. On taking absolute values in(8.27), we obtain(8.28)We will now show by induction that the solution of the difference equation(8.29)with ξ0 = 0 dominates the solution of (8.27); i.e., we will show thatξn > |en|n = 0, l, .
. .(8.30)Since e0 = ξ0 = 0, (8.30) is certainly true for n = 0. Assuming the truth of(8.30) for an integer n, it then follows from (8.29), since ξn > |en| and(1 + hL) > 1, thatξn+1 > |en+1|completing the induction.The solution ξn of the nonhomogeneous difference equation (8.29)therefore provides an upper bound for the discretization error en. From thetheory of difference equations given in Sec. 8.2, the solution of (8.29) isξ n = c(1 + hL) n - B(8.31)where c is an arbitrary constant, andTo satisfy the condition ξ0 = 0, we see that we must choose c = + B, sothat (8.31) becomesξ n = B(1 + hL) n - B8.4ERROR ESTIMATES AND CONVERGENCE OF EULER’S METHOD361We infer from Sec. 1.7 that ex = 1 + x + eξ x 2/2; hence ex > 1 + x, forall x. It follows that 1 + hL < ehL and therefore also that (1 + hL) n <enhL.
Using this in (8.31), we can therefore assert thatξ n < B( enhL - 1)where we have used the fact that nh = xn - x0. Since |en| < ξ n, we haveproved the following theorem.Theorem 8.2 Let yn be the approximate solution of (8.19) generated byEuler’s method (8.23). If the exact solution y(x) of (8.19) has acontinuous second derivative on the interval [x0, b], and if on thisinterval the inequalities|f y (x,y)| < L|y´´(x)| < Yare satisfied for fixed positive constants L and Y, the error en = y(xn)-yn of Euler’s method at a point xn = x0 + nh is bounded as follows:(8.32)This theorem shows that the error isthat is, the error tends tozero aslike ch for some constant c if x = xn is kept fixed.
It must beemphasized that the estimate (8.32) provides an upper bound rather than arealistic bound. Its primary importance is to establish convergence of themethod rather than to provide us with a realistic a priori error estimate.Example 8.3 Determine an upper bound for the discretization error of Euler’s method insolving the equation y´ = y, y(0) = 1 from x = 0 to x = 1.hence we can take L = 1. Also since y = ex,SOLUTION Here f(x,y) = y,xthen y´´ = e and |y´´(x)| < e for 0 < x < 1. To find a bound for the error at x = 1, wehave xn - x0 = 1, y = e1, and from (8.32)< 2.4hThus the error e(1) at x = 1 is bounded by 2.4h.
To see how realistic this bound is, weshall obtain the exact solution of Euler’s method for this problem. Thusy n + 1 = y n + hf(x n ,y n )= (1 + h) ynThe solution of this difference equation satisfying y(0) = 1 isyn = (1 + h) n362THE SOLUTION OF DIFFERENTIAL EQUATIONSNow if h = 0.1, n = 10, we find on expanding (1.1)10 that Euler’s method givesy 1 0 y (1) = 2.5937. On subtracting this from the exact solution y (1) = e = 2.71828,we find the error to be 0.1246, compared with the bound of 0.24 obtained by using(8.32).EXERCISES8.4-l For the equation y´ = - y,2 0 < x < 1, y(0) = 1:(a) Find an upper bound on the error at x = 1 in terms of the step size h, using (8.32).(b) Solve the difference equation which results from Euler’s method.(c) Compare the bound obtained from (a) with the actual error as obtained from (b) atx = 1 for h - 0.1, h = 0.01.(d) How small a step size h would have to be taken to produce six significant figures ofaccuracy at x = 1, using Euler’s method (assuming no round-off error)?8.4-2 The error en of an integration method is known to satisfy a difference inequality|en+2| < a1|en+1| + a2|en| + Awhere a1, a2, A are positive constants with e1 = e0 = 0.
Letequationbe a solution of the differencewith ξ 1 = ξ0 = 0. Show by induction that|en| < ξnfor all n8.5 RUNGE-KUTTA METHODSAs mentioned previously, Euler’s method is not very useful in practicalproblems because it requires a very small step size for reasonable accuracy.Taylor’s algorithm of higher order is unacceptable as a general-purposeprocedure because of the need to obtain higher total derivatives of y(x).The Runge-Kutta methods attempt to obtain greater accuracy, and at thesame time avoid the need for higher derivatives, by evaluating the functionf(x,y) at selected points on each subinterval.
We shall derive here thesimplest of the Runge-Kutta methods. A formula of the following form issought:yn+l = yn + ak1 + bk2where(8.33)k 1 = hf(x n ,y n )k 2 = h f(xn + αh,yn + βk 1)and a, b, α, β are constants to be determined so that (8.33) will agree withthe Taylor algorithm of as high an order as possible. On expanding y(xn+1)8.5RUNGE-KUTTA METHODS363in a Taylor series through terms of order h3, we obtain(8.34)where we have used the expansions (8.21), and the subscript n means thatall functions involved are to be evaluated at {xn,yn}.On the other hand, using Taylor’s expansion for functions of twovariables (see Sec. 1.7), we find thatwhere all derivatives are evaluated at {xn,yn}.If we now substitute this expression for k2 into (8.33) and note thatk1 = hf(xn,yn), we find upon rearrangement in powers of h thaty n + l = y n + (a + b)hf + bh 2 (αfx + βf f y )(8.34 a)On comparing this with (8.34) we see that to make the correspondingpowers of h and h2 agree we must havea + b = lb α = bβ = ½(8.35)Although we have four unknowns, we have only three equations, andhence we still have one degree of freedom in the solution of (8.35).
Wemight hope to use this additional degree of freedom to obtain agreement ofthe coefficients in the h 3 terms. It is obvious, however, that this isimpossible for all functions f(x,y).There are many solutions to (8.35), the simplest perhaps beinga=b=½α = β =1Algorithm 8.2: Runge-Kutta method of order 2 For the equationy´ = f(x,y)y(x 0 ) = y 0generate approximations y n to y(x 0 + nh), for h fixed and n =364THE SOLUTION OF DIFFERENTIAL EQUATIONS0, l, . . . , using the recursion formulayn+l = yn + 1/2(k1 + k2 ) with k1 = hf(xn,yn)(8.36)k 2 = hf(xn + h,yn + k1 )Algorithm 8.2 may be pictured geometrically as in Fig.
8.1. Euler’smethod yields an increment P1 P0 = hf(xn ,yn ) to yn ; P2 P0 = hf ( xn + h,yn+ hf(xn ,yn )) is another increment based on the slope obtained at xn+1.Taking the average of these increments leads to formula (8.36).The local error of (8.36) is of the formThe complexity of the coefficient in this error term is characteristic of allRunge-Kutta methods and constitutes one of the least desirable features ofsuch methods since local error estimates are very difficult to obtain. Thelocal error of (8.36), is, however, of order h 3 , whereas that of Euler’smethod is h2. We can therefore expect to be able to use a larger step sizewith (8.36).
The price we pay for this is that we must evaluate the functionf(x,y) twice for each step of the integration. Formulas of the Runge-Kuttatype for any order can be derived by the method used above. However, thederivations become exceedingly complicated. The most popular and mostcommonly used formula of this type is contained in Algorithm 8.3.Algorithm 8.3: Runge-Kutta method of order 4 For the equation y´ =f(x,y), y(x0 ) = y0 , generate approximations yn to y( x0 + nh) for hfixed and for n = 0, 1, 2, .
. . , using the recursion formulayn+1 = yn + 1/6(k1 + 2k2 + 2k3 + k4 )Figure 8.1(8.37)8.5RUNGE-KUTTA METHODS366whereThe local discretization error of Algorithm 8.3 isAgain theprice we pay for the favorable discretization error is that four functionevaluations are required per step. This price may be considerable incomputer time for those problems in which the function f(x,y) is complicated. The Runge-Kutta methods have additional disadvantages, whichwill be discussed later.
Formula (8.37) is widely used in practice withconsiderable success. It has the important advantage that it is self-starting:i.e., it requires only the value of y at a point x = xn to find y and y´ atx = xn+l.A general-purpose FORTRAN program based on Algorithm 8.2 for asingle differential equation is given below. To use this program, the usermust include a subroutine for evaluating the function f(x,y), and mustspecify the initial value y(x0) = y0, the final point x NSTEP, and the totalnumber of steps NSTEPS.FORTRAN PROGRAM FOR ALGORITHM 8.2CCCCCCFORTRAN PROGRAM TO SOLVE THE FIRST ORDER DIFFERENTIAL EQUATIONY´ (X) = F(X,Y)WITH INITIAL CONDITION OFY(XBEGIN) = YBEGINTO THE POINT XEND , USING THE SECOND ORDER RUNGE-KUTTA METHOD.A FUNCTION SUBPROGRAM CALLED 'F' MUST BE SUPPLIED.INTEGER I,N,NSTEPSREAL DERIV, H, Kl, K2, XBEGIN, XN, XEND, YBEGIN, YN1 READ 501, XBEGIN, YBEGIN, XEND, NSTEPS501 FORMAT(3Fl0.5,13)IF (NSTEPS .LT.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
