Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 45

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 45)

Then for any value of the parameter α,a2 = αa1 = 1 - αgives a formula with agreement in all terms involving h0. If we further require thatthe choicegives a formula with agreement in all terms involving h1. Thenyields a family of one-step methods of order 2 whenand f is sufficiently smooth.Some of the members of this family of formulas have names.

Euler’s methodhas α = 0 and the order p = 1. Heun’s method (the improved Euler method) is thecase α = l/2, and the midpoint or modified Euler method is the case α = 1. The broadapplicability of these formulas is seen when we ask what is needed for the convergencetheorem to be valid. The continuity of R obviously follows from that of f. It is apleasant fact that the Lipschitz condition on R also follows from that on f:for all 0 < h < h0, and we may take the Lipschitz constant for R to be226CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSTherefore, if the differential equation satisfies the conditions of Theorem 6.1, and if thefunction f has two continuous derivatives [which implies the solutionany member of the family withis convergent of order 2.Higher order procedures involving more substitutions can be derived in the sameway, although naturally the expansions become (very) tedious. As it happens, p evaluations off per step lead to procedures of order p for p = 1, 2, 3, 4 but not for 5.

For thisreason, fourth order formulas were often preferred for constant step size integrations.As in the second order case, there is a family of fourth order procedures depending onseveral parameters. The classical choice of parameters leads to the algorithmy0 = A,and for n = 0, l, . . . ,This is formulated for the first order system of equationsY(a) = AY´ = F(x,Y)in a natural way:Y0 = A,and for n = 0, l, . . . ,the classical Runge-Kutta algorithm is6.3 ONE-STEP METHODS227Another quite similar fourth order procedure is(6.21)There is little reason to prefer one of these procedures over the other as far as a singlestep is concerned.

In the next section we shall learn how R. England exploited (6.21)to achieve an estimate of the error made in the step.EXERCISES6.10 A definite integralcan be evaluated by solvwhat is the largest value of h for which the absoluteing an initial value problem for an ordinary differentialerror will be less than 10-6 in magnitude?equation. Let y(x) be the solution of6.13 Consider the linear equationy´ = f(x), a < x < by(a) = 0.Theny´ = P1(x)y + Q1(x).Show that the derivatives needed in Taylor seriesmethods can be obtained fromy(r) = Pr(x)y + Qr(x),whereRunge-Kutta methods integrate the more generalequation y´ = f(x,y). In this special case, they assumePr(x) = P´r-1(x) + P1(x)Pr-1(x)forms familiar in quadrature.

Referring to Chapter 5,Qr(x) = Q´r-1(x) + Q1(x)Pr-1(x), r = 2, 3, . . . .identify the familiar procedure to which both the clasUse this to develop a fifth order formula for computsical fourth order formula and England’s formula deing Dawson’s integral.generate.6.11 Implement Heun’s method to estimate the solution 6.14 An interesting fact about Runge-Kutta methods is thatof the initial value problem in Exercise 6.8b. Usethe error depends on the form of the equation as wellh = 1/40 and h = 1/80.

Compute the errors at x = 0.5as on the solution itself. To see an example of this,and x = 1.0 to see if they are roughly quartered as h isshow that y(x) = (x + 1)2 is the solution of each of thehalved. How small an h would you estimate is neededtwo problemsin order for the absolute error to be less than 10-6 iny´ = 2(x + 1),y(0) = 1magnitude?y(0) = 1.y´ = 2y/(x + 1),6.12 Implement the classical Runge-Kutta algorithm to esThen show that Heun’s method is exact for the firsttimate the solution of the initial value problem in Exequation. Prove that the method is not exact when apercise 6.8b.

Use h = 1/40 and h = 1/80. Computeplied to the second equation, although it has the samethe errors at x = 0.5 and x = 1.0. By what factor aresolution.the errors reduced as h is halved? Roughly speaking,228CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONS6.4 ERRORS-LOCAL AND GLOBALModern codes for the initial value problem do not use a fixed step size h. The errormade at each step is estimated and h is adjusted both to obtain an approximation that issufficiently accurate and to carry out the integration efficiently.

There is an unfortunateconfusion on the part of many users of codes with error estimates as to what is beingmeasured and what its relation to the true error is.The function y(x) denotes the unique solution of the problemy´ = f(x,y)y(a) = A.The true or global error at xn+l isy(x n+l ) - y n+1 .Unfortunately, it is relatively difficult and expensive to estimate this quantity. This isnot surprising, since in the step to xn+1 the numerical procedure is only supplied withxn,yn and the ability to evaluate f. The local solution at xn is the solution u(x) ofu´ = f(x,u)u(xn) = yn.The local error isThis is the errorinating at (xn,yn)reasonable to askhas on the globalmade approximating the solution of the differential equation origin a single step.

These errors are illustrated in Figure 6.1. It isthat a numerical procedure keep this error small. What effect thiserror depends on the differential equation itself. After all,(6.22)The quantityy(xn+1) - u(xn+1)is a measure of the stability of the differential equation since it is a consequence (atx n + l ) of the initial difference y(xn) - yn at xn. If this quantity increases greatly, theproblem is poorly posed or ill-conditioned or unstable.Example 6.7. Considery´ = αyfor a constant a. A little calculation shows that6.4 ERRORS-LOCAL AND GLOBAL229Figure 6.1 Local and global errors.furthermore,If α > 0, the solution curves spread out (Figure 6.2a), the more so as a is large.

Fromthe expression (6.23) it is clear that a small local error at every step does not imply asmall global error. On the other hand, if α < 0, the curves come together (Figure 6.2b)and (6.23) shows that controlling the local error will control the global error.

Forgeneral functions f(x,y) the Lipschitz condition alone cannot predict this behavior,nsince for this example the Lipschitz constant is |α| in either case.The local error is related to the local truncation error. Indeed, it is just h times thelocal truncation error, µ, for the local solution u(x):local error = u(xn+1) - yn+lFor example, when y(x) is a solution of y´ = f(x,y), we have seen that Euler’s methodhasApplying this to u(x), we havelocal error =230CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSFigure 6.2 Solution curves for (a) y´ = 2y and (b) y´ = -2y.Similarly, for the Runge-Kutta formulas of orderwe havethen the numerical approximations satisfyThis leads tolocal error =A little reflection about these expressions suggests a way to estimate the localerror.

Suppose we compute yn+1 with Euler’s method and we also compute an approximate solutionwith one of the second order Runge-Kutta formulas. Theexpressions above show thatThat is, the discrepancy between the two values estimates the error in the lower orderformula. This is the same principle used in Chapter 5 to estimate quadrature errors. Ingeneral, suppose that in addition to the valuewith truncation error µ n = O(hP), we compute another approximation6.4 ERRORS-LOCAL AND GLOBALwith truncation error231of higher order, q > p.

Then by definitionand, similarly,which, on subtracting, shows thatBecausegoes to zero faster thanwe can estimate the local error bylocal errorWe would like to approximate the local solution u(xn+1). In view of the fact thatwe have a good estimate of the error in yn+l, why not try to improve it by takingout the error? This process, called local extrapolation, is here formally equivalent tobecauseadvancing the integration with the higher order approximationThis tells us that local extrapolation will raise the effective order of the pair from p to q.Thus we can think of what is happening in two ways.

One is that a formula of order p isbeing used with its result improved by local extrapolation. The other is that a formulaof order q is being used with the step size selected conservatively by pretending thatthe step is being taken with a formula of a lower order p. Because local extrapolationincreases the accuracy at no increase in cost, all the current production-grade codesbased on explicit Runge-Kutta methods use it.A Runge-Kutta formula of order 4 requires (at least) four evaluations of F per stepand a companion formula of order 5 requires at least six. Just as with Gauss-Kronrodquadrature, the trick to being efficient is to derive the formulas as a pair in whichfunction evaluations are used in both formulas.

R. England published such a pair offormulas in [5]. To advance from xn to xn + h, he takes a step of length h/ 2 with (6.2 1)to get the fourth order resultand then another step of length h/2to get the fourth order resultBy working with two half-steps, he hasenough function evaluations available that with only one extra evaluation, he is ableto form a fifth order approximationIn this way, one extra functionevaluation is made every two half-steps to get an error estimate. An error estimate isused to control the local error and so give some credibility to the computed solution.It also allows the code to select the largest step size that will result in the error testbeing passed. Except in unusual cases, adaptation of the step size to the solution inthis way greatly increases the efficiency of integration.

Характеристики

Список файлов книги