Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 62
Текст из файла (страница 62)
1)STOPH = (XEND - XBEGIN)/NSTEPSXN = XBEGINYN = YBEGINDERIV = F(XN,YN)N =0PRINT 601, N, XN, YN, DERIV601 FORMAT(lX,I3,3E21.9)DO 10 N=l,NSTEPSKl = H*F(XN,YN)K2 = H*F(XN+H,YN+Kl)YN = YN + .5*(Kl+K2)XN = XBEGIN + N*HDERIV = F(XN,YN)10PRINT 601, N, XN, YN, DERIVGO TO 1ENDREAL FUNCTION F(X,Y)REAL X,YF = (1./X - Y)/X - Y*YRETURNEND366THE SOLUTION OF DIFFERENTIAL EQUATIONSExample 8.4 Solve the problem of Example 8.2 by the second-order Runge-Kuttamethod (8.36) and by the fourth-order Runge-Kutta method (8.37).In the machine results given in Sec.
8.3, (8.36) is called method 2 and (8.37) method3. We see that the second-order Runge-Kutta method gives results which are entirelycomparable with the Taylor algorithm of order 2 (method 1). The fourth-order RungeKutta method, however, yields remarkably improved results correct to six decimal placesfor h = 1/16 and to seven or eight places for other values of h. The computationalefficiency of methods 2 and 3 may be compared by considering the number of functionevaluations required for each. Method 2 requires two function evaluations per step andfor h = 1/128 requires in all 256 evaluations.
Method 3 requires four function evaluationsand for h = 1/16 a total of only 64 function evaluations and yet produces considerablymore accurate results. The fourth-order Runge-Kutta method is clearly a more efficientmethod to use for this problem, and this is generally true.EXERCISES8.5-1 For the equation y´ = x + y, y(0) = 1, calculate the local error of method (8.36).Compare this with the error of Taylor’s algorithm of order 2. Which would you expect to givebetter results over the interval [0, l]?8.5-2 Carry out a few steps of the integration of y´ = x + y, y (0) = 1, using (8.36) and a stepsize of h = 0.01; then write a program to solve this problem on a computer from x = 0 tox = 1.8.5-3 To Eqs. (8.35) add the additional condition that the coefficients of fxx in (8.34) and(8.34a) must agree.
Solve the resulting system of equations for a, b, α, β. Determine the errorterm of the second-order Runge-Kutta method obtained from this choice of a, b, α, β.8.5-4 It can be shown that the error of the fourth-order Runge-Kutta method satisfies for astep size h a relation of the formy n (h) - y(b) = A(b)h 4 +as h goes to zero, where b = x0 + nh, hence n(h) = (b - x0 ) / h and the constant A(b) doesnot depend on h. Use an extrapolation procedure as in the case of Romberg integration toobtain an approximation to y(b) for which the error is8.6 STEP SIZE CONTROL WITH RUNGE-KUTTA METHODSIn Section 8.5 we considered two Runge-Kutta (RK) methods, one oforder 2 and one of order 4.
Runge-Kutta methods of any order can bederived, although the derivation can become exceedingly complicated. Animportant consideration in using one-step methods of Runge-Kutta type isthat of estimating the local error and of selecting the proper step size toachieve a required accuracy. There is no reason why the step size h needsto be kept fixed over the entire interval as we did in Example 8.4.Estimating the accuracy using different fixed step sizes as we did inExample 8.4 may be very inefficient.
In this section we will examinemethods for estimating the local error and for varying the step sizeaccording to some error criterion.8.6STEP SIZE CONTROL WITH RUNGE-KUTTA METHODS361The first method is based on interval halving. Let us assume that weare using an RK method of order p and that we have arrived at a point xnwith h = xn - xn-1 .
We now integrate from xn to xn+1 = xn + h twice,once using the current step h and again using two steps of length h/2. Wewill thus obtain two estimates yh(xn+1) and yh/2(xn+1) of the value of y(x)at x = xn+1 and a comparison of these two estimates will yield an estimateof the error.
To derive the estimate we first note that a Runge-Kuttamethod of order p has a local asymptotic error expansion of the formy h ( xn + mh) = y(xn + mh) + C(xn + mh)hP +(8.38)Here, yh (xn + mh) denotes the approximation to the solution y(x) at thepoint x = xn + mh obtained after m h-steps of the Runge-Kutta method,starting from the exact value yn = y(xn). Further, the constant C( xn + mh)does not depend on h, though it does depend on f(x,y) and on the pointx = xn + mh. Therefore,y h (x n + 1 ) = y(x n + 1 ) + C(x n + 1 )h P +y h / 2 (x n + 1 ) = y(x n + 1 ) + C(x n + 1 )(h/2) P +(8.39 a)(8.39b)On subtracting (8.39a) from (8.39b) we find that the principal part of theerror in (8.39b) can be estimated asThe quantity(8.40)thus provides us with a computable estimate of the error in the approximation yh/2 (xn+1 ) and it can be used to help us decide whether the step hbeing used is just right, too big, or too small.Suppose now that we are given some local error tolerance ε and thatwe wish to keep the estimated error Dn below the local error tolerance perunit step, i.e., we wantDn < εh(8.41)Assume that we have computed yh(xn+1), yh/2(xn+1), and Dn.
We must nowdecide on whether to accept the value yh/2(xn+1) and on what step h to usefor the next integration. From the given error tolerance ε, we compute alower error bound ε´ < ε in a manner to be described later. We have thefollowing possibilities:In this case we accept the value yh/2(xn+1), and continue the integration from xn+1 using the same step’ size h.368THE SOLUTION OF DIFFERENTIAL EQUATIONS(ii)In this case the error is too large, hence we must reduce h—say to h/2—and integrate again from the point x = xn.(iii)In this case we are getting more accuracy than required.
We acceptthe value y h/2 (x n+1 ) replace h—by say 2h—and integrate from xn+1.If we restrict the interval step size to halving or doubling, then the lowerbound ε´ can be set toε ´ = ε/2P + 1for a pth order method since halving the step size reduces the order byapproximately 1/2P+1 . For the Runge-Kutta method of order 4 we havep = 4, hence ε´ = ε/32.
Actually, it is not advisable to change the step sizetoo often, and to be safe one might use ε´ = ε/50.A more sophisticated form of step size control, which does not restricth to doubling or halving, takes the following form. From (8.40) we have(8.42 a)Our goal is to choose a step size for the next step. Since the principalpart of the error at the next step will bewe must choose sothat the error tolerance (8.41) is satisfied, hence we must have(8.426)Assuming again that C n does not change much, we can eliminate C nbetween (8.42a) and (8.42b) as follows: From (8.42b) we have(8.42c)Thus if we have already successfully integrated with a step h, the nextintegration step size should be h or perhaps, to be safe, a little smaller. Asan example suppose that we have a method with p = 4, that ε = 10 -6 ,h = 0.1 and Dn is computed to be 10-5. Then8.6STEP SIZE CONTROL WITH RUNGE-KUTTA METHODS369These conditions would thus require a much smaller value of h. On theother hand, if again p = 4, h = 0.1, ε = 10-6 and we compute Dn = 10-8,thenso that the step size can be almost doubled.
The use of variable step sizesadds considerably to the complexity of a program and leads to results at aset of nonuniformly spaced points which to a user may be disconcerting.Halving and doubling intervals is generally more acceptable to the user.On the other hand, programs with automatic step size control provide theuser with very good estimates of accuracy, and are overall quite efficient.The major disadvantage of this method of error control is the substantial additional effort required.
In recent years several new variations ofRunge-Kutta methods suitable for step size control have been introduced.Some names associated with these new variations are Merson, Verner, andFehlberg. We describe briefly the method proposed by Fehlberg which wedenote by RKF 45 [28]. This method requires six function evaluations perstep but it provides an automatic error estimate and at the same timeproduces better accuracy than the standard fourth-order method. Fehlbergshowed that four of these function values combined with one set ofcoefficients could be used to produce a fourth-order method while all sixvalues combined with another set of coefficients could be used to producea fifth-order method. Comparison of the values produced by the fourthorder and fifth-order methods then leads to an estimate of the error whichcan be used for step size control.We describe very briefly the approach taken by Fehlberg.
We assumethat we have integrated the equation y´ = f(x,y) up to a point xn with astep size h, and we now wish to find an estimate of y(x) at x = xn+1. Oneestimate will be given by the formula(8.43a)for certain coefficients ci and a second estimate will be given by(8.43b)for another set of coefficients ci *. The error estimate for step size control isthen computed as follows:and it can be used as described earlier to estimate the proper step h for thenext integration.
The functions ki are the same in both formulas and can370THE SOLUTION OF DIFFERENTIAL EQUATIONSbe expressed in the formThere are many possible choices of the coefficients αi and β ij that will leadto Runge-Kutta methods of order 5. Fehlberg proposed one particular setof coefficients which we will not reproduce here. The interested reader isreferred to [28] for further details about this method.Another Runge-Kutta method with step size control, due to Verner, isthe basis of a very successful differential-equation-solving subroutinenamed DVERK which is widely available in subroutine libraries. Verner’smethod, which we denote by RKV 56, requires eight function evaluationsper step, and from these, two estimates of y(x) are obtained, one based ona fifth-order approximation and one based on a sixth-order approximation.A comparison of these two estimates then provides a basis for step sizeselection.
Some of the initial testing of this method was done at theUniversity of Toronto [29]. The method was later incorporated into thesubroutine DVERK and disseminated by IMSL Inc., Houston, Texas.IMSL, which stands for International Mathematical and Statistical Library,is a collection of thoroughly tested subroutines for a wide variety ofmathematical and statistical problems.
The library is available on a subscription basis and is available for almost all medium- and large-scalecomputers, including those of IBM, CDC, UNIVAC, Burroughs, andHoneywell. Since most computing installations now subscribe to the IMSLcollection, we shall not reproduce the code for DVERK here. Since we willuse this subroutine to solve several problems in this chapter, we willdescribe briefly the parameters in the call statement and the variousavailable options.In normal usage under default options and after initialization, theheart of the program to solve a first-order differential equation y´ = f(x,y)from x = XBEGIN to x = XM consists of a DO loop of the form:X = XBEGINY = YBEGINDO 10 K=l ,MXEND = XBEGIN + FLOAT(K)*(XM - XBEGIN)/FLOAT(M)CALL DVERK ( N, FCN1, Y, Y, XEND, TOL, IND, C, NW, W, IER )PRINT 600, XEND, Y(l), C(24)FORMAT(F19.6,E2l.8,Fl6.0)60010 CONTINUEThe parameters in the subroutine have the following meanings:N = the number of equations to be solved (here N = 1)FCN1 = the name of the subroutine for f(x,y); to be supplied by theuser as an external subprogramX = the initial value of the independent variable8.6 STEP SIZE CONTROL WITH RUNGE-KUTTA METHODS 371Y = the initial value of the dependent variableXEND = the value of x at which the solution is to be outputTOL = tolerance for error control; while different types of errortolerance specifications are possible, the default option triesto keep the relative global error less than TOLIND = 1 causes all default options to be used= 2 allows options to be selectedC = communications vector of length 24; some of these can be setby the user if IND was set to 2; these choices allow differenttypes of error control, minimum or maximum step sizes,limits on the number of function evaluations, etc.NW = the first dimension of the workspace matrix W, must be atleast as large as NW = workspace matrix whose first dimension is NW and whosesecond dimension must be greater than or equal to 9IER = an error flag, used to denote various types of errors encounteredIn the DO loop above, the points XEND are those values of x atwhich the solution is outputted.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
