Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 32
Текст из файла (страница 32)
Write aroutine Root that first finds a change of sign and thencalls Zero to compute a root. The parameter list ofRoot should be the same as that of Zero except thatB and C are replaced by arguments Z and SCALE.Here the Z input is a guess for a root. If all goes well,on output Z is to be the answer B obtained by Zero.The search algorithm you are to implement is essentially that of the fzero.m program of MATLAB.
However, fzero.m begins searching with an incrementZ/20 if Z 0 and 1/20 otherwise. In Root the initialsearch increment is to be an input variable SCALE.Initialize DZ = SCALE, B = Z - DZ, C = Z + DZ.If f(B)f(C) < 0 (be sure to code this properly), callZero to compute a root. Otherwise, double the increment, DZ = 2 × DZ, expand the search to the left byB = B - DZ, and test again. If this does not result ina change of sign, expand the search to the right byC = C + DZ and test again. If this does not resultin a change of sign, double the increment and repeat.Limit the number of tries. Test your code using one ofthe examples in the text.
After you are sure it works,you might want to use it in other exercises.CONDITION, LIMITING PRECISION, AND MULTIPLE ROOTSIt is important to ask what limitations on accuracy are imposed by finite precisionarithmetic. Since we seek a machine representable number that makesas nearlyzero as possible, the details of the computation off, the machine word length, and theroundoff characteristics play important roles. We have remarked that the computedfunction values may vary erratically in an interval about the zero and we have seen anexample in Figure 1.2.
Let us look at another example in more detail.Example 4.6. Consider the polynomial (x - 1)3 that we evaluate in the form ((x 3 ) x + 3)x - 1 in three-digit decimal chopped floating point arithmetic. For x = 1.00, 1.01,. . ., 1.17 the computed function values are exactly zero with the exception of the value0.0100 at x = 1.01, 1.11, 1.15. For x = 1.18, 1.19,. . . ,1.24 all function values are0.200 except for a value of exactly zero at x = 1.20 and a value of 0.0200 at x = 1.23.The reader might enjoy evaluating the function for x values less than 1 to explore thisphenomenon. It is clear that these erratic values might cause the secant rule to be unstable. Evaluating the derivative shows that Newton’s method can also be unstable atna multiple root like this one.What effect on the accuracy of a root does inaccuracy in the function values have?To get some feeling for this, suppose that the routine for f(x) actually returns a value158CHAPTER 4ROOTS OF NONLINEAR EQUATIONSand for x a machine number near a root a, the best we can say is thatfor a suitable e.
Suppose that z is a machine number andcan z be? If a is of multiplicity m, then= 0. How much in errorSince it is possible for f(z) to be as large as e, we could haveso it is possible that(4.14)For small ε and m > 1 the term ε1 / m is much larger than ε and there is a serious loss ofaccuracy. The other factor plays a role, but generally we must consider multiple rootsto be ill-conditioned (sensitive to errors in the evaluation off).
The ill conditioning ofthe root of multiplicity 3 in (x - 1)3 = 0 is evident in Example 4.6. We saw there thatx = 1.20 led to a function value of exactly zero, and this is certainly a poor approximation to the root at 1.00. The essence of the matter is that at a multiple root, f(x)is almost tangent to the horizontal axis so that shifting the curve vertically by a smallamount shifts its intersection with the axis by a considerable amount. Exercise 4.24 isan example of this effect.Even when m = 1, the root can be poorly determined.
As we have already seen,clusters of roots “look” like multiple roots from a distance, but we are now consideringwhat happens close to a root. Even when well separated and simple, a root is poorlydetermined if f(x) passes through zero with a small slope. More formally, the quantitycan be large when the slope |f´ (α)| is small. A famous example from Wilkinson [11,p. 43] illustrates this dramatically.Example 4.7.Consider the polynomial equation(x - 1)(x - 2)···(x - 19)(x - 20) = 0,which has the roots 1, 2,.
. . , 19, 20. These roots are obviously simple and well separated. The coefficient of xl9 is -210. If this coefficient is changed by 2-23 to become-210.000000119, the roots become those shown in Table 4.2.4.4 CONDITION, LIMITING PRECISION, AND MULTIPLE ROOTS159Notice that five pairs of roots have become complex with imaginary parts of substantial magnitude! There is really no remedy for this ill conditioning except to usenmore digits in the computations.Multiple roots are awkward not only because of their ill conditioning but for otherreasons, too. Bisection cannot compute roots of even multiplicity because there is nosign change. Its rate of convergence to roots of odd multiplicity is not affected bythe multiplicity, but the other methods we have presented slow down drastically whencomputing multiple roots.
If the derivative f´(x) is available, something can be doneabout both these difficulties. Near α,andwhereandThis says that zeros of f(x) of even multiplicity are zeros of f´(x) of odd multiplicity,so they could be computed with a bisection code or Zero. Also, notice thatwhereH(x) = g(x)/G(x), H(α) = l/m0,so that u(x) has only simple zeros. Because of this, solving u(x) = 0 with Zero isfaster than solving f(x) = 0 and allows the code to compute zeros of f(x) of even160CHAPTER 4ROOTS OF NONLINEAR EQUATIONSmultiplicity.
However, it must be appreciated that u(x) has a pole wherever f(x) = 0and f(x) 0.EXERCISES4.22 What is the value of the right-hand side of (4.14) forthe root in Exercise 4.10a?for the condition of this root. Perturb f(x) by 10-7using the form in Exercise 4.10e; then solve4.23 What is the value of the right-hand side of (4.14) forthe root of f(x) = (x- 10)(3x- 1)2 in [0,l]? Assumethat ε = 10-12.f(x) + 10-7 = 04.24 The problem f(x) = (x + 1)(x - 0.8)7 = 0 has 0.8 as aroot of multiplicity 7.
Evaluate the expression (4.14)4.5accurately with Zero (use ABSERR = RELERR =10-10). How much was the root 0.8 perturbed? Compare this to the result of (4.14) with ε = 10-7. Repeatusing the form in Exercise 4.10d.NONLINEAR SYSTEMS OF EQUATIONSA problem occurring quite frequently in computational mathematics is to find some orall of the solutions of a system of n simultaneous nonlinear equations in n unknowns.Such problems are generally much more difficult than a single equation. An obviousstarting point is to generalize the methods we have studied for the case n = 1.
Unfortunately, the bracketing property of the method of bisection does not hold for n > 1.There is a generalization of the secant method, but it is not at all obvious because of themore complicated geometry in higher dimensions. Newton’s method, however, doesextend nicely. Only the case n = 2 is examined because the notation is simpler and thebasic ideas are the same for the general case.Consider the system of equationsf(x,y) = 0g(x,y) = 0,(4.15)which we occasionally write in vector form asTo solve (4.15), Newton’s method for finding the root of a single nonlinear equation is generalized to two dimensions. The functions f(x,y) and g(x,y) are expandedabout a point (x0, y0) using Taylor’s theorem for functions of two variables (see the appendix for a statement of the theorem).
Carrying only terms of degrees 0 and 1 givesthe approximating system(4.16)4.5 NONLINEAR SYSTEMS OF EQUATIONS161which is linear in the variables x - x0 and y - y0. The next approximation (x1 ,y1) tothe solution of (4.16) is found by solving the linear equationsfor ∆x0 = x1 - x0 and ∆y0 = y1 - y0 and formingx1 = ∆x0 + x0 and y1 = ∆y0 + y0.In general (xk+1,yk+l) is obtained from (x k , y k) by adding a correctionobtained by solving a linear system. To summarize, we have derivedNewton’s method for two equations and two unknowns:or in matrix form,The matrix J is called the Jacobian matrix of the system of equations composed of fand g.Example 4.8.Set up Newton’s method for obtaining solutions to the equationsf(x,y) = x2 + xy3 - 9 = 0g(x,y) = 3x2y - y3 - 4 = 0.Sincethe system to be solved at each iteration isThe following table gives some numerical results for different starting points (x0, y0) .In all cases the iterations were terminated when the quantitywas less than 10-6.162CHAPTER 4ROOTS OF NONLINEAR EQUATIONSThese computations show that this system of equations has at least three solutions.Which solution is found depends on the starting guess (x0, y0).nAs with Newton’s method for a function of one variable, it can be shown that if his twice continuously differentiable near a root a of h(w) = 0, if the Jacobian matrix ata, J(a), is not singular, and if w0 is sufficiently close to a, then Newton’s method willconverge to a and the convergence will be quadratic.Advanced references like [7] develop Newton’s method much further.
A seriouspractical difficulty is to find an initial approximation sufficiently close to the desiredroot a that the method will converge. A knowledge of the problem and continuationcan be very helpful in this. A general approach is to connect finding a root w ofh(w) = 0 with minimizing the residual R(w) = f2 (w) + g2(w). Clearly this functionhas a minimum of zero at any root of h(w) = 0. The idea is to regard the change ∆w kcomputed from Newton’s method as giving a direction in which we search for a valueλ such that the iteratewk+l = wk + λ∆w kresults in a smaller value of the residual:This is always possible because until we reach a root,There are many practical details to be worked out.
For example, it is not necessary, oreven desirable, to find the value of λ that minimizes the residual. Methods of this kindare called damped Newton methods. A careful implementation will often convergewhen Newton’s method on its own will not and will have nearly the same efficiencywhen both converge.EXERCISESy2 - 14xz = 0.4.25 Use Newton iteration to find a solution (good to at-4least an absolute error of 10 in magnitude) near4.26 We seek the three parameters α, β, and γ in the model(0.5, l.0, 0.0) of the nonlinear system2x2- x + y2- z = 032x 2 - y 2 - 20z = 0by interpolating the three data points (1, 10), (2, 12),4.6 CASE STUDY 4and (3,18). Use Newton iteration to solve for the pa-163rameters to three significant figures.4.6 CASE STUDY 4The Lotka-Volterra equationsdescribing the populations of predator and prey are studied at length in most modernbooks on ordinary differential equations (see, e.g., Boyce and DiPrima [2]).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
