Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 66
Текст из файла (страница 66)
Moreover, as the step size was made smaller, these errors for a fixed value of xactually became larger rather than smaller. To illustrate this behavior, letus consider the method derived in Sec. 8.7,(8.66)y n + 1 = y n - l + 2hfnfor which the discretization error isWe would expect thismethod to give more accurate results for h fixed than Euler’s method,whose error isConsider, however, the following simple problem,y ´ = - 2y + 1(8.67)y(0) = 1whose exact solution is y = 1/2 e-2x + 1/2.The results given in Table 8.1 were obtained by the computer, using astep size of h = 1/32. The first column gives selected values of x at which thesolution is printed, Y(N) denotes the exact solution, Y1(N) denotes thesolution obtained by Euler’s method, Y2(N) the solution obtained by(8.66), and E1(N), E2(N) their respective errors.
Method (8.66) requiresTable 8.1X(N)Y(N)Y1(N)0.0.03125000.5000000l.00000001.50000002.00000002.25000002.50000003.00000003.50000003.75000003.78125003.81250003.84375003.87500003.90625003.93750003.96875004.00000001.00000000.96970650.68393970.56766760.52489350.50915780.50555450.50336900.50123940.50045590.50027650.50025980.50024400.50022930.50021540.50020230.50019010.50017850.50016771.00000000.96875000.67803700.56339430.52257300.50803760.50479620.50286200.50101900.50036280.50021650.50020290.50019030.50017840.50016720.50015680.50014700.50013780.5001292Y2(N)1.00000000.96970650.68408170.56782470.52513280.50970070.50642640.50479040.50507590.51086690.51743370.48199950.51968370.47953910.52224130.47675890.52514650.47361560.5284445E1(N)E2(N)0.-0.0009565-0.0059027-0.0042733-0.0023205-0.0011202-0.0007583-0.0005070-0.0002203-0.0000931-0.0000601-0.000568-0.0000538-0.0000509-0.0000482-0.0000456-0.0000431-0.0000408-0.00003860.-0.00000000.00014200.00015710.00023920.00054290.00087190.00142140.00383650.01041100.0171571-0.01826030.0194397-0.02069020.0220260-0.02344340.0249564-0.02656300.0282768390THE SOLUTION OF DIFFERENTIAL EQUATIONStwo starting values y0 and y1.
For y1 we take the exact value as computedfrom the exact solution.The error columns show that E2(N) is considerably smaller thanE1(N) for the first few steps but grows rapidly, so that at x = 2.25, E2(N)is greater than E1(N). Asthe solution approaches the steady-statevalue y = 1/2. Euler’s method actually approaches this steady-state solutionwith monotonically decreasing error, whereas for method (8.66) the error isgrowing exponentially. Moreover, as the last few steps (where the resultsare printed at every integration step) show, the errors E2(N) oscillate insign.
Beyond x = 4, Y2(N) would have no significant digits of accuracy.The phenomenon exhibited in this example is known as numerical instability.To help us understand this behavior, let us examine the differenceequation (8.66) more closely. For the example being considered, f n =-2yn + 1, and hence (8.66) becomes(8.68)y0 = ly n + 1 + 4h y n - y n - 1 = 2hWe can solve this difference equation explicitly, using the methods of Sec.8.2. The general solution of (8.68) is(8.69)where β 1, β 2 are the roots of the characteristic equationβ 2 + 4hβ - 1 = 0These roots areIf we expandin a Taylor’s series through linear terms, theseroots can be expressed in the formSubstituting into (8.69), we have(8.70)In the calculus it is shown thatUsing this limit and the fact that n = xn/h, it follows for xn fixed that*8.10STABILITY OF NUMERICAL METHODS391and similarly thatHence, asthe solution (8.70) approaches(8.71)Thus the first term tends to the true solution of the differential equation.The second term is extraneous and arises only because we have replaced afirst-order differential equation by a second-order difference equation.Imposing the initial conditions will, if all arithmetic operations are exact,result in choosing C2 = 0 so that the correct solution will be selected from(8.71).
In practice, however, some errors will be introduced, primarily dueto roundoff or to inexact starting values, and hence C2 will not be exactlyzero. A small error will therefore be introduced at each step of theintegration, and this error will subsequently be magnified because it isbeing multiplied by the exponentially increasing factorBecausethe major part of the true solution is exponentially decreasing, the errorintroduced from the extraneous solution will eventually dominate the truesolution and lead to completely incorrect results.Loosely speaking, we can say that a method is unstable if errorsintroduced into the calculations grow at an exponential rate as the computation proceeds.One-step methods like those of the Runge-Kutta type do not exhibitany numerical instability for h sufficiently small.
Multistep methods may,in some cases, be unstable for all values of h, and m other cases for a rangeof values of h. To determine whether a given multistep method is stablewe can proceed as follows: If the multistep method leads to a differenceequation of order k, find the roots of the characteristic equation corresponding to the homogeneous difference equation. Call these β i (i =1, . . . , k). The general solution of the homogeneous difference equation isthen(8.72)One of these solutions, say β 1n, will tend to the exact solution of thedifferential equation asAll the other solutions are extraneous.
Amultistep method is defined to be strongly stable if the extraneous rootssatisfy asthe conditioni = 2, 3, . . . , k|β i | < 1Under these conditions any errors introduced into the computation willdecay as n increases, whereas if any of the extraneous β i are greater thanone in magnitude, the errors will grow exponentially.For the general differential equation y´ = f(x,y), it will be impossibleto obtain the roots β i of the characteristic equation. A consideration of the392THE SOLUTION OF DIFFERENTIAL EQUATIONSspecial equation y´ = λy, λ constant, is usually considered sufficient, however, to give an indication of the stability of a method.We consider first the Adams-Bashforth fourth-order method. If in(8.47) we set f(x,y) = λy we obtain(8.73)The characteristic equation for this difference equation isThe roots of this equation are of course functions of hλ. It is customary towrite the characteristic equation in the form(8.74)where ρ(β) and σ(β) are polynomials defined byρ ( β ) = β4 - β 3We see that as(8.74) reduces to ρ(β) = 0, whose roots are β 1 = 1,β 2 = β 3 = β 4 = 0.
Forthe general solution of (8.73) will have theformwhere the β i are solutions of (8.74). It can be shown that approaches thedesired solution of y´ = λy aswhile the other roots correspond toextraneous solutions. Since the roots of (8.74) are continuous functions ofh, it follows that for h small enough, |β i | < 1 for i = 2, 3, 4, and hencefrom the definition of stability that the Adams-Bashforth method isstrongly stable. All multistep methods lead to a characteristic equation inthe form (8.74) whose left-hand side is sometimes called the stabilitypolynomial.
The definition of stability can be recast in terms of the stabilitypolynomial. A method is strongly stable if all the roots of ρ(β) = 0 havemagnitude less than one except for the simple root β = 1.We investigate next the stability properties of Milne’s method (8.64b )given by(8.75)Again setting f(x,y) = λy we obtainand its characteristic equation becomesρ ( β ) + hλσ(β) = 0(8.76)*8.10STABILITY OF NUMERICAL METHODS393ρ(β) = β 2 - 1withσ(β) = β 2 + 4β + 1This time ρ(β) = 0 has the roots β 1 = 1, β 2 = -1, and hence by thedefinition above, Milne’s method is not strongly stable.
To see the implications of this we compute the roots of the stability polynomial (8.76). For hsmall we have(8.77)Hence the general solution of (8.75) isIf we set n = xn/h and letthis solution approaches(8.78)In this case stability depends upon the sign of λ. If λ > 0 so that thedesired solution is exponentially increasing, it is clear that the extraneoussolution will be exponentially decreasing so that Milne’s method will bestable. On the other hand if λ < 0, then Milne’s method will be unstablesince the extraneous solution will be exponentially increasing and willeventually swamp the desired solution. Methods of this type whose stability depends upon the sign of λ for the test equation y´ = λy are said to beweakly stable.
For the more general equation y´ = f(x,y) we can expectweak stability from Milne’s method wheneveron the interval ofintegration.In practice all multistep methods will exhibit some instability for somerange of values of the step h. Consider, for example, the Adams-Bashforthmethod of order 2 defined byIf we apply this method to the test equation y´ = λy, we will obtain thedifference equationand from this the stability polynomialor the equation394THE SOLUTION OF DIFFERENTIAL EQUATIONSIf λ < 0, the roots of this quadratic equation are both less than one inmagnitude provided that -1 < hλ < 0.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
