Algebraic Equations: Solving with Bisection and Newton-Raphson Methods, Study notes of Physics

Download Algebraic Equations: Solving with Bisection and Newton-Raphson Methods and more Study notes Physics in PDF only on Docsity! Chapter 4 Roots of Algebraic Equations In this chapter we will • learn how to solve algebraic equations using the – Bisection method – Newton-Raphson method – Secant method • introduce the concept of convergence of iterative methods One of the most obvious uses of computer algorithms in the physical and mathe- matical sciences is in solving algebraic equations. A computational method often provides the only way of obtaining the roots of an equation, since analytic solutions exist only for a handful of cases. In other situations, a computational method may provide a faster route in obtaining the roots of an equation that has known analytic solutions. In this chapter, we will consider general equations of the form f(x) = 0 , (4.1) where f(x) is any function of the unknown variable x. Every equations can be trivially put in this form by subtracting its right-hand side from its left-hand side. The function f(x) may be an analytic expression in closed form, such as f(x) = xex − 1 . (4.2) Alternatively, it may be a complicated computational procedure that returns a real number for each value of the argument x. In this chapter, we will consider only equations that involve real quantities and search for roots that are real numbers. A root x0 of an equation f(x) = 0 is said to have multiplicity k if there is a Root Multiplicityfunction g(x) such that f(x) = (x − x0) kg(x) . (4.3) Alternatively, a root x0 of an equation f(x) = 0 is said to have a multiplicity k if the k−th order derivative of the function f(x) evaluated at x0 is zero, i.e., if dkf(x) dxk ∣ ∣ ∣ ∣ x0 = 0 . (4.4) 4-1 4-2 CHAPTER 4. ROOTS OF ALGEBRAIC EQUATIONS Figure 4.1: An example of functions that have a root of multiplicity (left) one, (center) two, and (right) three. Figure 4.1 shows a few illustrative examples of functions with roots of multiplicity one, two, and three. In general, a double root corresponds to a function that is tangent to the x−axis, whereas a root of multiplicity three or higher corresponds to a function that has an inflection point on the x−axis. 4.1 The Bisection Method This is the simplest and most robust method for finding the root of an equation, provided that we can guess an initial interval in which one and only one root exists. The algorithm successively divides the interval in half (bisects) keeping the solution within its limits, until it reaches a desired level of accuracy. The main drawback of the method is its slow rate of convergence. The bisection algorithm is based on the intermediate value theorem: Consider a function f(x) that is continuous in the interval [a, b] with f(a)·f(b) < 0.Intermediate Value Theorem Then, there exists a value x0 in the interval [a, b] such that f(x0) = 0. In other words, this theorem states that, if a continuous function f(x) changes sign in the interval [a, b], then there is at least one root of the equation f(x) = 0 within this interval. It is worth paying close attention to the requirement that the function f(x) is continuous in the interval [a, b]. If it is not, then the theorem does not guarantee the existence of a solution. Consider, for example, the function f(x) = 1 x − 2 . (4.5) Clearly, the function changes sign in the interval [0, 4]. However, we cannot apply the intermediate value theorem for this function, because it is discontinuous at x = 2. Indeed, there exists no value of x for which f(x) = 0. When using the bisection algorithm to find iteratively the solution of an equation of the form f(x) = 0 we take the following steps: 1. We start with an initial interval [a, b] in which we know a priori that the solution exists 4.1. THE BISECTION METHOD 4-5 Figure 4.2: The function f(x) = xex − 1 plotted against the variable x. A quick glance reveals that the equation f(x) = 0 has a unique solution in the interval [−1, 1]. Finding the initial interval in which a unique solution to an equation lies is more Guessing the Initial Intervalof a guess than an algorithm. In most cases, calculating the limiting values of the function at ±∞ or examining the graph of the function gives the best handle on bracketing its solutions. These two techniques are illustrated with the following example. Example: Solving the equation xex−1 = 0 with the bisection method In order to study whether this equation has roots, we will first check the limits of the function f(x) = xex − 1 as x → ±∞. Calculating one of the limits is simple, because lim x→∞ (xex − 1) = +∞ . (4.12) Calculating the limiting value of the function as x → −∞, however, requires the use of L’Hospital’s rule lim x→−∞ (xex − 1) = lim x→−∞ (xex) − 1 = lim x→−∞ ( x e−x ) − 1 = −1 . (4.13) Because the two limits have opposite signs and the function is continuous, the equation has at least one solution for some real value of the variable x. If the function is monotonic and has two limits of opposing signs, then it has at most one solution. We study the monotonicity of the function by evaluating its derivative df dx = (x + 1)ex . (4.14) This is positive when x > −1 and negative when x < −1. At x = −1, the function takes the value f(−1) ≃ −1.368, which is negative. As a result, between x = −∞ and x = −1 the function is negative and decreasing, it has a minimum at x = −1, and between x = −1 and x = +∞ the function is increasing and changes sign. Therefore, the equation has only one solution in the interval [−1, +∞). In order to obtain a finite upper bound on the interval in which the solution exists, we recognize the fact that if x > 1, f(x) = xex − 1 > 1 · e1 − 1 > 1 . (4.15) 4-6 CHAPTER 4. ROOTS OF ALGEBRAIC EQUATIONS The function, therefore, changes sign in the interval [−1, 1] in which the single root of the equation lies. Alternatively, we can guess an initial interval for the bisection algorithm by plotting the function and studying graphically its behavior. A quick glance at the figure shows that the equation has only one solution, which lies in the interval [−1, 1], as we inferred earlier using analytical methods. Inserting this interval as an initial guess in a bisection function allows us to find the solution to the equation as x0 = 0.567143 with an accuracy of 1 × 10 −6. The bisection method is a very robust and simple to implement algorithm forAssessment of the Bisection Method finding the root of an equation. It guarantees finding the solution at the required accuracy, within a finite number of steps. The method, however, has a number of disadvantages: (i) It cannot find the solutions to all equations. For example, a polynomial that has a root of multiplicity 2 does not change sign across the solution. More specifically, the equation (x − 2)2 = 0 (4.16) cannot be solved with the bisection method because the function f(x) = (x − 2)2 is positive or zero for all real values of x. The same is also true for many equations that are not polynomials of x. For example, the equation ex − x − 1 = 0 (4.17) has an obvious root at x = 0, even though the function f(x) = ex − x − 1 is never negative. (ii) It requires an a priori knowledge of an interval in which a unique solution exists. This not necessarily a disadvantage of the method, because it forces us to study thoroughly the equation before we solve it and, hence, protects us from making mistakes. It is worth asking, however, which root does the bisection method find when there are more than one roots in the initial interval. The answer, of course, depends on the particular equation and the choice of initial interval. One can answer this question probabilistically, however, considering a large number of application of the method to various equations with various initial intervals[1]. The result is fascinating! Consider an equation f(x) = 0 that has N simple roots in the interval [a, b], which we will denote by x1 < x2 < ... < xN . The bisection algorithm will find the odd numbered roots (x1, x3, etc.) with equal probability and the even numbered roots (x2, x4, etc.) with zero probability. It is fair to say that these are not good chances to play with! (iii) It converges slowly. Although the uncertainty in the solution is reduced by a factor of two after each iteration, the bisection method is the slowest among the root finding algorithms that we discuss in this chapter. Before embarking into a more detailed study of this issue, however, we need to discuss in general the order of convergence of an iterative method. 4.2 Order of Convergence of an Iterative Method When we employ an iterative method to solve a problem, we start with an initial guess of the solution and improve it systematically with each cycle of the algorithm. 4.3. THE NEWTON-RAPHSON METHOD 4-7 Our aims is to reach a given accuracy for the solution with the smallest number of cycles. We often measure the efficiency of an iterative method to achieve this by calculating its order of convergence. For concreteness, we will consider an iterative method that aims to find the root of an algebraic equation x0 and denote by xi, the approximate solution after i iterations. We will also use the symbol ǫi to denote the fractional error between the approximate and the correct solution after i iterations. If a constant number λ exists such that lim i→∞ |xi+1 − x0| |xi − x0|N = λ (4.18) then we say that our iterative method converged with an order N and an asymptotic error constant λ. If N = 1, we call the convergence linear. If N = 2, we call the convergence quadratic. The higher the order of convergence of an iterative method, the faster the so- lution converges to a required accuracy. Indeed, after the first few iterations, the definition (4.18) implies that the error in the solution achieved by a method with order of convergence N is ǫi+1 ≃ λǫ N i . (4.19) We can use this definition to calculate the order of convergence of the bisection method. We will denote by ∆ ≡ b − a the size of the initial interval in which the single root exists. This is also an estimate of the uncertainty ǫ0 with which we know the solution a priori . After one iteration, the interval and hence the uncertainty is reduced to half, i.e., ǫ1 = ∆/2. After two iterations, it is reduced by one more factor of two, i.e., ǫ2 = (∆/2)/2 = ∆/2 2, and so on. In general, ǫi+1 = 1 2 ǫi (4.20) and, therefore, the bisection method converges linearly with an asymptotic error constant of 0.5. We can use the linear order of converge of the bisection method to calculate how many iterations we will need to complete in order to reach an uncertainty δ in the final estimate of the solution. It is easy to show by induction that after i iterations the uncertainty is ǫi = ∆/2 i. Setting this to δ and solving for i we obtain i = log10 ∆ − log10 δ log10 2 . (4.21) We can, therefore, reduce the uncertainty with which we know the solution to an algebraic equation by six orders magnitude by using the bisection method for 20 iterations. 4.3 The Newton-Raphson Method The Newton-Raphson method is an iterative algorithm for solving algebraic equa- tions of the form f(x) = 0. It overcomes two disadvantages of the bisection method. It does not require knowledge of an interval in which a unique solution lies and it converges quadratically. However, it does not guarantee convergence to a solution and requires an analytic knowledge of the first derivative of the function f(x). In order to introduce this method, we will consider a function f(x) that is continuous and such that the equation f(x) = 0 has a root at x = x0. We will also assume that we have made an initial guess for the value of the root, which we will 4-10 CHAPTER 4. ROOTS OF ALGEBRAIC EQUATIONS double foo(double x); // prototypes for functions double fooprime(double x); double newtonraphson(double root) /* Uses the Newton-Raphson method to find the solution to an algebraic equation of the form f(x)=0. The function f(x) and its derivative are provided by the functions foo(x) and fooprime(x) and are supplied by the user. The parameter ACCURACY is used in determining whether the method has converged to a solution. Upon completion, it returns the root of the equation. */ #define ACCURACY 1.0e-6 { double delta; // Newton-Raphson correction double fx,fprime; // f(x) and its derivative fx=foo(root); // f(x) at x=root do { fprime=fooprime(root); // f’(x) at x=root if (fabs(fprime)>ACCURACY) // if not extremum/inflection delta=-fx/fprime; // calculate correction else // otherwise { // print error message and exit printf("Newton-Raphson cannot converge\n"); return 0.0; } root+=delta; // update root fx=foo(root); // f(x) at root } // until converged while (fabs(delta)>ACCURACY || fabs(fx)>ACCURACY); return root; } method converges quadratically. In order to show this, we will consider two suc- cessive approximations xi and xi+1 evaluated after the i-th iteration and the one following it. These two trials are related by equation (4.23) or xi+1 − xi = − f(xi) f ′(xi) . (4.24) If we denote by ǫi ≡ xi − x0 and ǫi+1 ≡ xi+1 − x0 the deviations of the two approximations from the true solution x0, then ǫi+1 = ǫi − f(xi) f ′(xi) . (4.25) We now Taylor expand the function f(x) around the true solution x0 to obtain f(x) = f(x0) + f ′(x0)(x − x0) + 1 2 f ′′(x0)(x − x0) 2 + ... = f ′(x0)(x − x0) + 1 2 f ′′(x0)(x − x0) 2 + ... , (4.26) where we have used the fact that f(x0) = 0. We then take the derivative of the expanded function with respect to the variable x and evaluate it at xi to obtain f ′(xi) = f ′(x0) + f ′′(x0)(xi − x0) + ... (4.27) 4.3. THE NEWTON-RAPHSON METHOD 4-11 Inserting equations (4.26) and (4.27) into equation (4.25) we obtain ǫi+1 = ǫi − f ′(x0)ǫi + 1 2f ′′(x0)ǫ 2 i f ′(x0)[1 + f ′′(x0)ǫi f ′(x0) ] . (4.28) Note that this last equation is value only if f ′(x0) 6= 0, which is the reason why we required the multiplicity of the solution to be one. We now expand the denominator in the fraction using the approximate relation (1 + x)a ≃ 1 + ax + ... (4.29) that is valid when x ≪ 1, rearrangement some terms keeping only those that are up to second order in the small parameter ǫi and arrive at the final result ǫi+1 = ǫ 2 i [ f ′′(x0) 2f ′(x0) ] . (4.30) This last relation proves that the Newton-Raphson method converges quadratically to a solution of multiplicity one of an algebraic equation. This is, of course, true as long as the initial guess is not significantly different than the true solution. If the multiplicity of the root is equal to m > 1, then it can be shown that the Newton-Raphson method converges only linearly with an asymptotic error constant of l = 1 1 − m . (4.31) In general, the Newton-Raphson method is an algorithm that converges fast Assessment of the Newton-Raphson method to the solution of an equation starting from a simple initial guess. It also does not require an a priori knowledge of an interval in which the solution exists. This method, however, also has a number of disadvantages. (i) It requires an analytic knowledge of the derivative f ′(x). In many situations this may not be possible, if the function f(x) itself is the result of a numerical calculation. In other cases, calculating the derivative may be computationally so expensive that it counters the gain in speed achieved by the quadratic convergence of the method. (ii) It does not converge globally. Indeed, our proof of quadratic convergence relies on two assumptions. That the current guess to the solution is not very differ- ence from the true value and that the derivative of the function f ′(x) at the guess is not zero. If either of these conditions is not satisfied, the method may not converge at all to the true solution. (iii) In very symmetric situations, it enters a cyclic behavior without converging to a solution. Further Reading 1. Which root does the bisection algorithm find? Corliss, G., SIAM Review 19, 2, 325 (1977) 2. Historical Development of the Newton-Raphson Method , Ypma, T. J., SIAM Re- view 37, 4, 531 (1995)