Nonlinear Equations

Chapter 4 Nonlinear equations 4.1 Root finding Consider the problem of solving any nonlinear relation g(x) = h(x) in the real variable x. We rephrase this problem as one of finding the zero (root) of a function, here f (x) = g(x) − h(x). The minimal assumption we need on f, g, h is that they’re continuous. We have at our disposal any number of evaluations of f and its derivative ' f. 1. Method 1: bisection. The bisection methods starts from two points a0 and b0 such that f (a0 ) > 0, and f (b0 ) < 0. Because f is continous, there must exist a root in the interval [a0 , b0 ]. At stage k, assume that we have obtained an interval [ak , bk ] such that the same sign properties hold: f (ak ) > 0 and f (bk ) < 0. The bisection method consists in subdividing the interval [ak , bk ] in two and discard the half in which there may not be a root. Let mk = (ak + bk )/2. • If f (mk ) < 0, then it is the interval [ak , mk ] which is of interest. We put ak+1 = ak and bk+1 = mk . • If f (mk ) > 0, then it is the interval [mk , bk ] which is of interest. We put ak+1 = mk and bk+1 = bk . • If f (mk ) = 0, then mk is a root and we stop the algorithm. 1 CHAPTER 4. NONLINEAR EQUATIONS In practice, this iteration is stopped once f (mk ) gets small enough. Let x∗ be the unknown root. The error obeys |x∗ − mk | ≤ |bk − ak | = 2−k |b0 − a0 |. Every step of the bisection discovers a new correct digit in the binary expansion of x∗ . The advantage of the bisection method is that it is guaranteed to con­ verge to a root, by construction. On the other hand, convergence is rather show compared to the next 2 methods we now present. If there are several roots, the bisection method will converge toward one of them (we may not have no control over which root the method chooses.) 2. Method 2: Newton-Raphson. This method is very important: it is the basis of most optimization solvers in science and engineering. Let us first present the Newton-Raphson method for solving a single scalar equation f (x) = 0. Newton’s method fits a tangent line to the point (xn , f (xn )) on the graph of f , and defines xn+1 at the intersection of this tangent line with the x axis. We have 0 = f (xn ) + (xn+1 − xn )f ' (xn ), from which we isolate xn+1 = xn − f (xn ) . f ' (xn ) For instance, we can find the decimal expansion of positive root of f (x) = x2 − 2. The iteration reads xn+1 = xn − √ 2 by finding the x2n − 2 xn 1 = + . 2 xn 2xn Starting with x0 = 1, we get x1 = 32 = 1.5, x2 = 17 = 1.4167..., 12 √ 577 x3 = 408 = 1.4142157... The true value of 2 is 1.4142135... Convergence is very fast, when it occurs. Assume that f '' is continuous, and that f ' (x) = 0 in some neighborhood of the root x∗ (large enough so that all our iterates stay in this neighborhood.) Put En = xn − x∗ . 2 4.1. ROOT FINDING Then we can perform a Taylor expansion of f around xn , and evaluate it at x = x∗ : 1 0 = f (x∗ ) = f (xn ) + (x∗ − xn )f ' (xn ) + (x∗ − xn )2 f '' (ξ), 2 for some ξ ∈ int(xn , x∗ ) (the notation int refers to the interval generated by xn and x∗ , i.e., either [xn , x∗ ] or [x∗ , xn ].) We also have the equation defining xn+1 : 0 = f (xn ) + (xn+1 − xn )f ' (xn ). Subtracting those 2 equations, we get 1 0 = −En+1 f ' (xn ) + E2n f '' (ξ), 2 En+1 = 1 f '' (ξ) 2 E . 2 f ' (xn ) n rr Our assumptions ensure that the ratio ffr (x(ξ) exists and converges to n) ' ∗ '' ∗ some limit (f (x )/f (x )) as n → ∞. Hence the sequence is bounded uniformly in n, and we can write |En+1 | ≤ CE2n , where C > 0 is some number (which depends on f but not on n.) It follows that 1 k |En | ≤ (CE0 )2 . C We say the method “converges quadratically” because the exponent of En is 2. The number of correct digits is squared at each iteration. In contrast, the bisection method only converges linearly. We also some­ times refer to “linear convergence” as first-order convergence, although the meaning of the expression is completely different from what is was in the previous chapters. Convergence is ensured as soon as the starting point x0 is close enough k to the (unknown) root x∗ , in the sense that |CE0 | < 1, so that (CE0 )2 → 0 as k → ∞. If the condition |CE0 | < 1 is not satisfied, Newton’s method may very well diverge. For instance, we expect problems when the derivative is very small: following the tangent can take us to a region very far away from the root. An example of a function f (x) for 3 CHAPTER 4. NONLINEAR EQUATIONS which Newton’s method diverges is atan(x), when x0 is chosen to be too far from the origin. On the plus side, Newton’s method is fast. On the minus side, Newton’s method only converges to a root only when you’re already quite close to it. 3. Method 3: the secant method. If we do not know the derivative, we cannot set up Newton’s method, but we can approximate it by replacing the derivative by (let fn = f (xn )) fn − fn−1 f [xn−1 , xn ] = . xn − xn−1 Hence we define xn+1 by xn+1 = xn − fn f [xn−1 , xn ] . The geometrical idea is to replace the tangent line at xn by the secant line supported by xn−1 and xn . The secant method requires two points x0 and x1 as starting guesses. Notice that at each step, only one evaluation of f is necessary, because f (xn−1 ) is already known from the previous iteration. If we were to form a finite difference approximation of the derivative with a very small grid step h, we may be more accurate but that requires two evaluations of f rather than one. Let us check the convergence properties of the secant method. The line joining the two points (xn−1 , fn−1 ) and (xn , fn ) is the degree-1 interpolant in the interval [xn−1 , xn ]: p(x) = fn + f [xn−1 , xn ](x − xn ). Outside of this interval, it is an extrapolant. Regardless of whether x ∈ [xn−1 , xn ] or not, the difference between p and f is known from a theorem we saw in the previous chapter: 1 f (x) − p(x) = f '' (ξ)(x − xn )(x − xn−1 ), 2 4 4.1. ROOT FINDING where ξ is in the smallest interval containing x, xn−1 , and xn . Evaluat­ ing this relation at the root x = x∗ , we get 1 0 = fn + f [xn−1 , xn ](x∗ − xn ) + f '' (ξ)(x∗ − xn )(x∗ − xn−1 ). 2 On the other hand the definition of xn+1 gives 0 = fn + f [xn−1 , xn ](xn+1 − xn ). Subtracting these two equations we get En+1 = 1 f '' (ξ) En En−1 . 2 f [xn−1 , xn ] Again, thanks to the same assumptions on f as in Newton’s method, f rr (ξ) the ratio f [xn−1 ,xn ] has a finite limit as n → ∞, hence is bounded by some number C > 0. We get |En+1 | ≤ C|En ||En−1 |. The decay of En is somewhere between first (linear) and second (quadratic) order. To obtain a more precise rate of decay, we guess that the in­ equality above should be reducible to the form |En | ≤ C|En−1 |p for some p. Using this equation and |En+1 | ≤ C|En |p above, we get 2 |En−1 |p ≤ C|En−1 |p |En−1 |. The exponents match on the left and the right provided p2 = p + 1, which has for positive solution √ 1+ 5 p= . 2 (a number sometimes called the golden ratio). We check that p = 1.618..., a number between 1 and 2, Hence the secant method is faster than bisection, but slower than Newton’s method. The secant method inherits the problem of Newton’s method: it only converges when the starting guesses x0 and x1 are sufficiently close to the root. 5 CHAPTER 4. NONLINEAR EQUATIONS We can also set up Newton’s method in several dimensions. A system of nonlinear equations is of the form fi (x1 , . . . , xn ) = 0, i = 1, . . . , n. We take the same number of equations and unknowns, so that we may be in a situation where there is one solution (rather than a continuum of solutions or no solution at all.) Whether the system has zero, one or several solutions is still a question that needs to be addressed separately. The shorthand notation for the system is f (x) = 0 By analogy with the 1D case we perform a Taylor expansion about xn : 0 = f (x∗ ) = f (xn ) + Vf (xn )(x∗ − xn ) + O(�x∗ − xn �2 ). With indices this equation is written as n : : ∂fi 0 = fi (x ) = fi (xn ) + (xn )(xj,n − x∗j ) + O( (xj,n − xj∗ )2 ). ∂xj j=1 j ∗ (Watch the subscript n which indicates the n-th iterate while the subscript j indicates the j-th component.) The next iterate xn+1 is defined by neglecting the quadratic error and isolating x∗ . A linear system of equations needs to be solved: the Jacobian matrix Vf (xn ) is inverted and we get xn+1 = xn − [Vf (xn )]−1 f (xn ). The geometrical interpretation of this equation is that we can fit the tangent plane to each of the surfaces y = fi (x1 , . . . , xn ) in Rn+1 , find the line at the intersection of all these planes, and check where this line intersects the (hyper)plane y = 0. Newton’s method is still quadratically convergent in multiple dimensions, and special care must still be taken to make sure that we start close enough to a root. Example 10. x21 + x22 = 1, x2 = sin(x1 ). Write this as a root-finding problem: f1 = f2 = 0 with f1 (x1 , x2 ) = x21 + x22 − 1, 6 f2 (x1 , x2 ) = x2 − sin(x1 ). 4.2. OPTIMIZATION PROBLEMS The Jacobian matrix is � J = Vf (x) = ∂f1 ∂x1 ∂f2 ∂x1 ∂f1 ∂x2 ∂f2 ∂x2 � 2x1 2x2 . − cos(x1 ) 1 = Use the formula for the inverse of a 2-by-2 matrix: a b c d −1 = 1 ad − bc d −b , −c a to obtain J −1 = 1 2x1 + 2x2 cos(x1 ) 1 −2x2 . cos(x1 ) 2x1 The Newton iteration is therefore x1,n+1 x2,n+1 4.2 = x21,n + x22,n − 1 x1,n − J −1 . x2,n x2,n − sin x1,n Optimization problems Another important recurring problem in science and engineering is that of finding a minimum or a maximum of a function F (x). A point x∗ is a local minimum when F (y) ≥ F (x∗ ) for all y in a neighborhood of x∗ . It is a global minimum when F (y) ≥ F (x∗ ) for all y. We write min F (x) x for the minimum value F (x∗ ). We then call x∗ the argument of the minimum. Maximizing F (x) instead is the same as minimizing −F (x), so it suffices to talk about minimization. When F (x) is smooth, and x is allowed to run over all real numbers (not restricted to an interval or other set), then it suffices to solve F ' (x) = 0 (and check that F '' (x) > 0) in order to find a local minimum. Hence it suffices to apply Newton’s method or any other root-finding method to the function f (x) = F ' (x). We obtain xn+1 = xn − F ' (xn ) . F '' (xn ) 7 CHAPTER 4. NONLINEAR EQUATIONS In multiple dimensions, we minimize a scalar function F (x1 , . . . , xn ). The optimality condition, obeyed at the minimum x∗1 , . . . x∗n , is that all the partial derivatives of F vanish, i.e., VF (x∗1 , . . . x∗n ) = 0. Newton’s method, also called Newton descent, follows from considering these ∂F . We get equations as a nonlinear system fi (x1 , . . . , xn ) = 0 with fi = ∂x i xn+1 = xn − [VVF (xn )]−1 VF (xn ). The matrix VVF of second partial derivatives of F is called the Hessian. In index notation, ∂ 2F (VVF )ij = . ∂xi ∂xj Compare Newton’s method with simple gradient descent: xn+1 = xn − αVF (xn ), for some sufficiently small scalar α. Gradient descent is slower but typically converges from a larger set of initial guesses than Newton’s method. Example 11. Consider F (x1 , x2 ) = x21 + (log x2 )2 . This function has a unique minimum for x1 ∈ R and x2 > 0. We compute   2x1 ∇F (x1 , x2 ) = 2 log x2 x2 and ∇∇F (x1 , x2 ) = Newton’s iteration is therefore     x1,n+1 x1,n = − x2,n+1 x2,n 1 2 0 2 0 ! 0 2−2 log x2 x22 0 x22,n 2−2 log x2,n ! . 2x1,n x2,n 2 logx2,n ! . Notice that x1 goes in one step to zero, because a quadratic function is exactly minimized by Newton’s method. 8