Thompson - Computing for Scientists and Engineers (523188), страница 52
Текст из файла (страница 52)
Later, Newton about 1680, Huygens about 1691, and Taylor (of Taylor’s theorem) in 17 15also provided proofs.After the Great Fire of London in 1666, Hooke and Christopher Wren were architectural consultants for the rebuilding of St. Paul’s cathedral, and Hooke probably told Wren about the stability of the catenary shape for arches and domes. For asolid arch, in the approximation that it is continuous and of uniform composition,the tension forces will all be tangential to the surface of the arch if it is of the regularcatenary shape. This property will improve the stability of the arch.
It is unlikely,however, that Wren used his colleague’s advice when designing the cathedral dome.8.3 NUMERICAL METHODS FOR SECOND-ORDERDIFFERENTIAL EQUATIONSIn the next four sections of this chapter we emphasize numerical methods for solvingsecond-order differential equations, with several worked examples and applicationsof the methods. Beginning in this section, we fist extend the Euler-type algorithmsfrom the first-order equations developed in Sections 7.4 and 7.5 to second-orderequations.
Then we set the scene for the Noumerov method by showing how firstderivatives can be removed from second-order linear equations, since the Noumerovalgorithm assumes that first derivatives are absent. This algorithm is then derived.Project 8A (Section 8.4) is on programming and testing the Euler-type algorithms, while Project 8B (Section 8.5) emphasizes programming and testing theNoumerov algorithm, then applying it to solve numerically the differential equationfor the quantum harmonic oscillator. The problem of the numerics of stiff differential equations in Section 8.6 completes our study of numerical methods for secondorder equations.280SECOND-ORDER DIFFERENTIAL EQUATIONSEuler-type algorithms for second-order equationsWe now describe three Euler-type algorithms that are variants of the forward-difference predictor formula (7.38) used for first-order equations.
As we see in Section 8.4, these algorithms are not very accurate for a given stepsize h. Their utilitylies in their general applicability, because they make very few assumptions about thetype of second-order differential equation that is to be solved numerically.The Euler methods for second-order equations are based on the fact that the second-order differential equation(8.54)is equivalent to a pair of first-order equations, namely(8.55)and, as an identity,(8.56)Therefore, in principle, any second-order equation can be solved as a pair of firstorder equations.
In the practice of numerical methods for estimating the solutions ofsuch equations, careful analysis and programming, as well as considerable skill andexperience, are necessary if accurate results are to be obtained efficiently. The difficulties arise because with a finite stepsize h the numerical methods are necessarilyapproximate, as (7.38) indicates.We now describe three variants of the forward Euler predictor that may be usedfor estimating the solutions of the differential-equation pair (8.55) and (8.56).These methods are generally of increasing accuracy, as indicated by the dependenceof their error estimates on the stepsize h.Method 1. Here we use the forward Euler predictors directly for advancing thesolutions y and y' according to (7.38) namely(8.57)and(8.58)The first equation neglects the change in F between steps k and k + 1, while thesecond neglects the change of y' in this interval, so that it neglects a term of orderFh.2 The errors in both the function and its derivative are therefore of order h2.8.3 NUMERICAL METHODS FOR SECOND-ORDER DIFFERENTIAL EQUATIONS281Method 2.
In this variant we advance the derivative according to (8.57), but wethen average the derivative just obtained with the derivative at the start of the intervalin order to predict the new y value. That is,(8.59)The averaging of the two derivatives allows for what physicists would call an acceleration, since if x is time and y is displacement, then y' is speed and F is acceleration. The order of the error in y is h3 in this method.Method 3. Here the function is advanced by using terms up through its secondderivative, that is, through F, so its error is of order h3. Then the derivative is advanced using the average of the F values at the two ends, which gives an error in thederivative of order h3.
Thus(8.60)then the derivative is advanced by(8.61)Notice that the order of evaluation is significant here if F is a function of y. Further,in (8.57) — (8.61) mathematically there should besigns or computationally theyshould be assignment statements (“:=“ in Pascal).The three Euler-type methods are summarized in Table 8.3. The programmingof these algorithms is presented in Project 8A in Section 8.4.TABLE 8.3 Euler-type algorithms for numerical solution of second-order differential equaions.In method 3 the formulas must be evaluated in the order shown. The dependence of the errors onstepsize, h, is also indicated.282SECOND-ORDER DIFFERENTIAL EQUATIONSExercise 8.18Make Taylor expansions of y and of y' about xk in order to derive (8.57)through (8.61).
Include enough terms in each expansion that the lowest neglected power of h is indicated. Use this to show that the error estimate for each formula has the approximate dependence on h that is claimed. nA simple example will clarify these Euler-type algorithms for second-order equations. Suppose, for purposes of comparison, that we know the exact solution(8.62)Let us compare the exact solutions with those obtained from the three methods inTable 8.3 when they are used to advance the solutions from x1 = 1 to x2 = 2 in asingle step with h = 1, assuming that all values are exact at x = 1. Working backwards, we have by differentiating (8.62) twice with respect to x that(8.63)Starting with y1 = 1.6667, y'1 = 2.5000, and F1 = 2.0000, the three methodsproduce the values given in Table 8.4 and displayed in Figures 8.13 and 8.14.TABLE 8.4 Example of the Euler-type methods for advancing one step of the differentialequation solution.
The second derivative is (8.63) and the exact solution is (8.62). The value ofx = 2.In order to check that you understand how each of the Euler-type methodsworks, try the following exercise.Exercise 8.19Use each of the algorithms in Table 8.3 in turn to predict the values of y2 andUsing that h = 1 and thaty2, given the exact values in the text for y1 andthe neglected terms can be calculated from the derivatives F' = 1 and F” = 0,are the discrepancies between estimated and exact values correct? nThe behavior of the algorithms and their results are also shown in Figures 8.13and 8.14 for y and its derivative y' respectively. Notice that the constancy of the8.3 NUMERICAL METHODS FOR SECOND-ORDER DIFFERENTIAL EQUATIONS283derivative of F gives rise to estimates of y2 that are the same (but not exact) formethods 2 and 3.
For similar reasons, methods 1 and 2 agree (inexactly) on thevalue of y', while method 3 predicts this derivative exactly because only the secondand higher derivatives of F (which are zero) have been neglected.FIGURE 8.13 Numerical solutions by the three Euler methods for the exponential function,which is shown as the solid curve, y (x).This completes our derivation and preliminary discussion of the Euler-type algorithms for second-order differential equations. Numerical analyses that are muchmore extensive are part of Project 8A in Section 8.4.FIGURE 8.14 Numerical solutions by the three Euler methods for the derivative of the exponential function, whose derivative is shown as the solid curve, y'(x).284SECOND-ORDER DIFFERENTIAL EQUATIONSRemoving first derivatives from second-order linear equationsThe Noumerov algorithm for second-order linear differential equations that we derive in the next subsection requires that the first derivative, y', be absent.
Thiswould seem to be a severe limitation of the algorithm, because it might not be able tohandle such interesting differential equations as that for a damped harmonic oscillator (Section 8.1) where the damping is proportional to y'. The purpose of thissubsection is to show how the first derivative in a second-order linear differentialequation can always be transformed away.Consider the general second-order linear differential equation(8.64)The method of eliminating the first derivative is to find a transformation that makesthe first two terms part of a single second derivative. Consider transforming y(x)by multiplication by some, as yet unknown, function g(x), to produce z(x)(8.65)(Since unknowns are piling up at a great rate, this is clearly mathematics.) Also consider multiplication of the original equation (8.64) throughout by g(x). Is there achoice of g(x) that will remove the first derivative?Exercise 8.20(a) Show that(8.66)and therefore that if we force the first- and second-derivative terms in this equation and g times the original equation to coincide, we must have(8.67)(b) Show that the solution of this equation is(8.68)in which any constant of integration may be ignored because it just affects theoverall scale of the differential equation solution.(c) Thence show that the transformed variable z satisfies the following lineardifferential equation, in which its second derivative does not occur(8.69)8.3 NUMERICAL METHODS FOR SECOND-ORDER DIFFERENTIAL EQUATIONS285(d) Express the second derivative g” in this equation by using (8.67) in order toderive the final form of the transformed differential equation(8.70)in which all variables are implicitly functions of x.
nThus, if the original equation (8.64) has a first-derivative term, then the function gcalculated from (8.68) produces a differential equation without a first derivative,namely (8.70) for z. From the solution of this equation the original function, y, canbe recovered by dividing z by g because of (8.65).An example will make this formal analysis much clearer. Consider (8.64) forthe damped harmonic oscillator, with the coefficient of the second derivative dividedout for simplicity:(8.7 1)In the notation of (8.64) we have a(x) = b, a constant, so that a' = 0.b(x) = k, again a constant. From (8.68) the function g is therefore justAlso,(8.72)The modified differential equation (8.70) becomes(8.73)where K = k - b2 /4.
By inspection, the solution for z is just a complex exponential (assuming that K is positive). Division by g produces the solution of the originaldamped-oscillator equation(8.74)which is in agreement with (8.13). You will notice that elimination of the first derivative makes the differential equation much easier to solve analytically. As the Noumerov algorithm that we now derive shows, a similar result obtains for numericalsolution of these differential equations.Deriving the Noumerov algorithm for second-order equationsBoris V. Noumerov (1891 - 1943), a prominent Russian astronomer and geophysicist, derived the following algorithm as part of his numerical studies of the perturbations of planetary motions.
The algorithm was first published in English in the reference given at the end of this chapter. It is accurate and efficient for the solution of286SECOND-ORDER DIFFERENTIAL EQUATIONSsecond-order linear equations with no first-derivative term. As we showed in thepreceding subsection, the latter restriction is not really an impediment.We seek a recurrence relation relating y and its second derivative at a given xk interms of its values at the two preceding points, namely(8.75)where A, B, C, D, and E are to be chosen to minimize the remainder R.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
