Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 44
Текст из файла (страница 44)
The true values of the integraly(xn) are taken from [7].To study the convergence of Euler’s method, we relate the error at xn+1 to the errorat xn. Subtracting (6.9) from (6.8) givesDenoting the error at xn by En = y(xn) - yn, the Lipschitz condition on f and thisequation imply thatWith the definitionwe obtain(6.10)Here the term h2M2/2 bounds the error made in the current step and the other termbounds the error propagated from preceding steps.To prove convergence, we bound the worst error that can arise as we step fromx0 = a to xN = b and then show that it tends to zero as h does.
The first order ofbusiness is to see how rapidly the inequality (6.10) allows the error to grow. To do thiswe establish a general result for later use. Suppose there are numbers δ > 0 and M > 0such that the sequence d0, d1 , . . . satisfies6.2 A SIMPLE NUMERICAL SCHEME219The case n = 0,d1 < (1 + δ)d0 + M,can be combined with the case n = 1 to obtaind2 < (l + δ)dl + M < (1 + δ) 2 d0 + M[1 + (1 + δ)].Similarly,d 3 < ( l + δ )d2 + M < (1 + d)3d0 + M[1 + (1 + δ) + (1 + δ) 2 ].At this point we might guess thatdn < (1 + δ) n d0 + M[1 + (1 + δ) + (1 + δ)2+ ··· +(1 + δ) n-1 ].(6.11)To prove this, we use induction.
The inequality (6.11) certainly holds for n = 1, 2, 3.Suppose the inequality is true for the case n = k. Thendk + 1 < (1 + δ)dk + M< (1 + δ)k+1 d0 + M[1 + (1 + δ) + ··· + (1 + δ) k ],which establishes the result for n = k + 1 and completes the induction argument.Lemma 6.1. Suppose there are numbers δ > 0 and M > 0 such that the sequenced0, d1 , . . . satisfiesdk+l < (l + δ) dk + M, k = 0, l, . . .
.Then for any n > 0,(6.12)Proof. Using the identitywith x = 1 + δ, we see that the right-hand side of (6.11) can be rewritten in the form(6.13)Expansion of the exponential function about zero gives for δ > 0,It then follows thatand220CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSThis implies that (6.13) is bounded bynwhich establishes (6.12).Returning now to Euler’s method, we apply Lemma 6.1 to (6.10) and arrive atHowever, nh = xn - a and E0 = y0 - A = 0, so(6.14)Using xn - a < b - a, this implies that(6.15)It is seen that the error of Euler’s method is bounded by a constant times h.
When thevalue of the constant is immaterial, such expressions are written as 0(h).Generally speaking, we try to ignore the effects of finite precision arithmetic.When the tolerance corresponds to a relative accuracy comparable to the word lengthof the computer, this is not possible. Also, if the solution is very hard to approximateaccurately at xn, the step size necessary may be so small that the precision must beconsidered. To gain some insight, note that we do not obtain f(xn,yn) from a subroutine, but rather f(xn,yn) + εn.
Similarly, in computing yn+1 = yn + h[ f(xn,yn) + εn] madditional error ρn is made. The sequence generated computationally is thenyn+l = yn + hf(xn,yn) + hεn + ρn .Let us suppose that |ρn| < ρ and |ρn| < ε for all h < h0. Then the analysis can bemodified to yieldAccording to this bound, roundoff effects get worse as the step size is reduced inan attempt to get a more accurate solution. Clearly there is a maximum accuracypossible that depends on the problem, the numerical method, and the arithmetic of thecomputer used. The effects are more complex than this bound shows, but the boundis qualitatively correct. It is easy to show experimentally that as h is decreased, thenumerical solution is at first more accurate, reaches a best value, and subsequently isincreasingly less accurate.The convergence analysis just presented is the traditional one.
The trouble is thatthis is not the way modern codes work. Rather than take a step of specified lengthh, they select a step size automatically that will produce a solution with a specifiedaccuracy. A reasonable model of the step sizes selected by such codes is that at xnthe code selects a step hnwhereis a piecewise-continuous function6.3 ONE-STEP METHODS221satisfying 0 < q << 1 for a < x < b. With this model it is easy enough to modifythe convergence proof just given to account for variable step size.
The result is that asthe maximum step size H tends to zero,It is not hardto see how the model comes about. The user specifies a tolerance τ. In a step of lengthh from xn, Euler’s method makes an error of approximately h2|y´´(xn)|/2. The largeststep size hn that can be used and still keep the error less than τ is then aboutSpecial rules come into play in the codes when y´´(xn,) is nearly zero, so suppose thaty´´(x) does not vanish in [a,b]. IfandthendefinesNotice that H = O(τ 1/2) so that max |y(xn) - yn| is O(τ 1/2) for Euler’smethod with automatic selection of the step size.EXERCISES6.8 Use Euler’s method on the following problems usinga fixed step size h = 1.0, and then h = 0.5. In eachcase calculate the errors at x = 1.0.(a) y´ = -y/( x + 1) with y(0) = 1, so y(x) = l/(x +1).(b) y´ = -y 3/2 with y(0) = 1, so y(x)6.9 Implement Euler’s method to estimate solutions of theinitial value problem in Exercise 6.8b.
Use h = l/40and h = l/80. Compute the errors at x = 0.5 andx = 1.0 to see if they are roughly halved as h is. Howsmall an h would you estimate is needed in order forthe absolute error to be less than 10-6 in magnitude?6.3 ONE-STEP METHODSLet us now consider general one-step methods and base our assumptions on the successful treatment of Euler’s method. The recipe is to be of the formy0 = A(6.16)222CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSThe method has no memory, so Φ depends only on the arguments listed. Usually fand h are omitted in the notation.
It is assumed that Φ is continuous in x and y. Thetreatment of Euler’s method had Φ(x,y) = f(x,y) and a Lipschitz condition was usedin an important way. So, for the general procedure we assume that(6.17)for a < x < b, all 0 < h < h0 for some h0, any continuous function f satisfying aLipschitz condition, and all u, v.In discussing Euler’s method we used as a starting point the fact that the solutiony(x) almost satisfies the recipe (6.9) used to define the numerical approximation.
Theanalog here is(6.18)with µn “small.” More precisely, if for all xn in [a,b] and all h < h0, there are constantsC and p such that(6.19)then we say that the method is of order p for the equation (6.1). The quantity µ n iscalled the local truncation error.Theorem 6.2.Suppose the initial value problemy´ = f(x,y)y(a) = Aon the finite interval [a,b] is solved by the one-step method (6.16) and suppose that thehypotheses of Theorem 6.1 are satisfied. If Φ(x,y) satisfies (6.17) and if the method isof order p > 1 for y(x), then for any xn = a + nh[a,b],ProofAs before, let En = y(xn) - yn and subtract (6.16) from (6.18) to obtainUsing the Lipschitz condition (6.17) and the fact that the method is of order p, we seethatThe theorem now follows from Lemma 6.1 and the fact that E0 = 0.nAs with our discussion of Euler’s method, the result of this theorem gives convergence of O(hP).
This explains our calling the method of order p for y(x). The term “amethod of order p´´ is used to describe a method that is of order p if f is sufficientlysmooth. The order of convergence is lower when f is not so smooth.As explained in connection with Euler’s method, codes select the step size automatically so as to keep the error smaller than a tolerance τ at each step. At the same6.3 ONE-STEP METHODS223time they try to use an efficiently large step. A reasonable model of such a step sizealgorithm leads to a step size hn at xn given byfor a piecewise-continuous function Θ(x) with 0 < θ < Θ(x) < 1 on [a,b] With stepsizes specified in this way, the convergence proof can be altered easily to conclude thatthe error is O(HP) = 0(τ1 / P ) .The most important task now left is to find functions Φ that are inexpensive toevaluate and of order p for smooth f.
We need, then,with µ n = O(hp). A Taylor series expansion of y(x) shows thatif y(x)case thatSo we find that if the method is of order p, then it must be thewith z(x) = O(hP). Because y(x) is a solution of the differential equation y´(x) =f(x,y(x)), the derivatives of y can be expressed in terms of total derivatives of f.Using the notation f(m) (x,y(x)) for the mth total derivative of f and subscripts forpartial derivatives, the relation isy(m)(x)=f(m-1)(x,y(x)),whereThe expression for Φ(x,y)becomes(6.20)An obvious choice for Φ is the function T(x,y),which yields a family of one-step methods called the Taylor series methods. Euler’smethod is the case p = 1.
When it is possible to evaluate derivatives efficiently, thesemethods are very effective.224CHAPTER 6ORDINARY DIFFERENTIAL EQUATIONSFor a simple equation like that satisfied by Dawson’s integral, and especially whenhigh accuracy is desired, a Taylor series method may be the best way to proceed. Thisequation hasf(x,y) = 1 - 2xy(1)f (x,y) = -2.x + (4x2 - 2)y,and so forth.
Exercise 6.13 develops a simple recursion that makes it easy to use aTaylor series method of very high order for this equation.Runge-Kutta methods use several evaluations of f(x,y) in a linear combination toapproximate y(x). The simplest case is Euler’s method that uses one evaluation. Letus now derive a procedure using the two evaluations f(xn,yn) and f(xn + p1h,yn +p2hf(xn,yn)), where p1 and p2 are parameters. Then for Φ we use the linear combination R(x,y):R(xn ,yn ) = a1 f(xn,yn) + a2 f(xn + p1 h,yn + p 2 hf(xn ,yn )).In this expression we are free to choose any useful values for p1, p2, a1, and a2.
Theaim is to choose the parameters so that the representation (6.19) holds for as large avalue of p as possible. To carry this out we expand all the quantities in Taylor seriesin h and equate coefficients of like powers. To simplify the notation, arguments areshown only if they are different from (xn,yn). The reader familiar with Taylor seriesexpansions of a function of two variables may skip to the result. Otherwise, we canproceed by a succession of familiar one-variable expansions as follows:Now we want to choose the parameters so thator, writing this out,6.3 ONE-STEP METHODS225Equating terms involving the same powers of h, it is found that it is possible to obtainagreement only for h0 and h1 :h0 : a1 + a2 = 1h1 : a2p2 = ½a2p1 = ½.Let a2 = α.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
