Download Newton's Method: Real and Complex Cases and more Study notes Introduction to Business Management in PDF only on Docsity!
Theory of Equations
Lesson 6
by
Barry H. Dayton
Northeastern Illinois University
Chicago, IL 60625, USA
www.neiu.edu/˜bhdayton/theq/
These notes are copyrighted by Barry Dayton, 2002. The PDF files are freely available on the web and may be copied into your hard drive or other suitable electronic storage devices. These files may be shared or distributed. Single copies may be printed for personal use but must list the website www.neiu.edu/˜bhdayton/theq/ on the title page along with this paragraph.
“Maple” and “MAPLE” represent registered trademarks of Waterloo Maple Inc.
2.11. NEWTON’S METHOD, REAL CASE 71
2.11 Newton’s Method, Real Case
Today, the method of choice for finding roots of polynomials is Newton’s method. This method is sometimes known under other names such as the “Newton-Raphson” method, but I believe it is proper to call this Newton’s method as it first appears in full generality in Newton’s “Method of Fluxions.” In the case of quadratics of the form x^2 − a this method (whose result is
a) was known to the Babylonians. Newton’s method is well adapted to the calculator or computer, is fast, will gener- ally be as accurate as the internal accuracy of the calculator or computer and, although this method does not work as well on roots of multiplicity greater than 1, it is the only algorithm that we will study extensively that works at all for multiple roots.
The algorithm goes as follows. Let f (x) be your polynomial, choose an approxi- mate root x 0. Then (hopefully) a better approximation is x 1 = x 0 − f (x 0 )/f ′(x 0 ). A still better approximation is x 2 = x 1 − f (x 1 )/f ′(x 1 ). In general we get a sequence of approximations defined by
xn+1 = xn −
f (xn) f ′(xn)
After a relatively small number of iterations, the sequence {xn} will usually converge to a root.
Consider the example we used to illustrate Horner’s method. Here f (x) = x^4 − 2 x^3 + x^2 − 3 , f ′(x) = 4x^3 − 6 x^2 + 2x. Choose as a starting point x 0 = 2. Then
x 1 = 2 −. 083333333 =1. 916666667 x 2 =1. 916666667 − .008723748=1. 907942919 x 3 =1. 907942919 − .000089648=1. 907853271 x 4 =1. 907853271 − .000000009=1. 907853262
At this point our computer tells us f (x 4 ) = 0 so the terms xn have stabilized at x 4 and so we may stop.
One way we may justify this method is graphically. Given a value x 0 relatively close to a root with y 0 = f (x 0 ) we construct the tangent line to the graph at (x 0 , y 0 ). The slope is f ′(x 0 ) and using the point-slope form of the equation of a line the equation of the tangent line is
(y − y 0 ) = f ′(x 0 )(x − x 0 ).
If x 0 is close to a root the tangent line crosses the x-axis even closer to that root. To find this point, set y = 0 in the equation above and solve for x. We get f ′(x 0 )x =
2.11. NEWTON’S METHOD, REAL CASE 73
instance, after a change of variables we can assume c = 0 and f (x) = xm^ + · · ·. For x close to 0 we can ignore the higher order terms so we can take f (x) = xm. A calculation shows that the iteration formula is xn+1 = ((m − 1)/m) ∗ xn or xn = ((m − 1)/m)n^ ∗ x 0. As (m − 1)/m < 1 we know that ((m − 1)/m)n^ goes to 0 as n gets large and thus Newton’s method will converge to the root. For a root of multiplicity two the error is halved each time so Newton’s method is about the same speed as the bisection method (as applied to simple roots, recall that bisection will not work for roots of even multiplicity), for higher multiplicity convergence is even slower. There are variations of Newton’s method which will speed convergence in this case, or better yet one can use the Euclidean algorithm to reduce your polynomial to one where the multiplicities of the roots are smaller. As a practical matter, however, polynomials with multiple roots are rare in nature and thus this problem will not occur often.
Maple Implementation
MAPLE does not have a built in Newton’s method, using instead more sophisticated solving methods. It is easy enough to write your own Newton’s method procedure to watch the method converge.
Newt := proc(f,z0,n) local g,z,j; g := x - f/diff(f,x); z := evalf(z0); for j to n do z := evalf(subs(x=z,g)); print(z); end do; end;
Here f is the polynomial in indeterminate x whose root is to be found, z0 is the initial point and n is the number of iterations (10 is a good number to start with). As I mentioned earlier in regards to iteration, you are warned to not use an indeterminate or algebraic expression as your initial value!
Exercise 2.11.1 [20 points] Find all 5 real roots of f (x) = x^5 − 8 x^3 + 10x + 1 to 6 significant digits using Newton’s method. For each root give the initial point used and the number of iterations needed.
2.12 Newton’ Method, Complex Case
Perhaps the best thing about Newton’s method is that not only does the algorithm work in the complex case but it works much better in the complex plane, converging promptly to a root from almost any initial point. The method is the same:
- Choose an initial point z 0.
- Iterate using zn+1 = zn − f (zn)/f ′(zn).
Here z is the variable we often use to denote a complex number and the arithmetic here is that of complex numbers. Starting with an imaginary initial point the sequence {zn} almost always will converge quickly to a root (real or imaginary) of f (z). In fact, if f (z) has real coefficients Newton’s method will often converge more quickly to a real root starting from an imaginary initial point. The reason is that the real Newton’s method can get hung up on a critical point (i.e. a zero of f ′(z)) whereas the complex Newton’s method has room to go around it. It is more difficult to estimate the location of imaginary roots than real roots. For- tunately with the complex Newton’s method this is not so vital. Using Newton’s upper bound one can find upper bounds on the modulus of complex roots. Often convergence is more direct if one starts with an initial estimate slightly larger in modulus than the desired root. It should be remembered that for polynomials with real coefficients, the imaginary roots come in conjugate pairs. The complex Newton’s method is easiest to implement using a computer language such as MAPLE or FORTRAN which supports complex arithmetic, however with care the algorithm can be implemented using other languages. Generally implementations of the complex Newton’s method use Horner’s process to calculate f (z) and f ′(z) efficiently. This method works the same in the complex case as in the real case. Sometimes because Horner’s process is used, algorithms to do the complex Newton’s method are, incorrectly, called Horner’s algorithm. When Newton’s method finds a root c of f (z) then f (z) factors as f (z) = (z − c)q(z). q(z) is the first quotient calculated by Horner’s process when calculating f (c). The method of ”deflation”, once a root is found, replaces f (z) with q(z). Then after finding a root of q(z) deflation is used again and so on, each time we get a polynomial of lower degree to solve. By the time a linear or quadratic polynomial is reached then it can be solved directly. In this way one gets quickly and in an orderly fashion all the complex roots of a polynomial. The quotients q(z) obtained by deflation are only approximate quotients since the roots found are only approximate roots, thus one may wish to go
Exercise 2.12.2 [10 points] Use the complex Newton’s method with deflation (if you use Maple you can deflate manually) to calculate the 7 complex roots of p(x) = 2 − 3 x + x^4 + 8x^7. In a previous Exercise we showed all roots α satisfy |α| < 1. For each root α of p(x) find an initial point z 0 with |z 0 | > 1 for which Newton’s method applied to p(x), not a deflation, converges to α and give the number of iterations required.
2.13 The Newton–Barstow algorithm
Because the complex Newton’s method requires complex arithmetic it is a little bit awkward to implement in a computer language, such as versions of BASIC or C, which does not support complex arithmetic. Since the coefficients of the original polynomial are usually real the Newton-Barstow algorithm which can find imaginary roots of real polynomials using only real arithmetic is sometimes an appealing alternative. The idea behind the Newton-Barstow method is interesting, although a bit com- plicated. We only sketch it here. We start with a real polynomial p(x) so imaginary roots occur in conjugate pairs. In particular each imaginary root corresponds to a real quadratic factor of p(x), thus the goal of this algorithm is to find such factors of p(x). Thus we start with a guess x^2 + ax + b of a quadratic factor. Dividing p(x) by this guess we get a quotient and a linear remainder cx + d. We want this remainder to be 0 so the procedure is to refine our initial guess in order to get a smaller remainder. We can regard the coefficients of the remainder c, d as functions of the two variables a, b. By a variation on Horner’s procedure we can calculate the partial derivatives of c, d in terms of a, b. We can then apply a two dimensional real version of Newton’s method to produce new values of a, b which will tend to give smaller values of c, d. Iterating this procedure c, d become small enough to be considered 0 and so x^2 + ax + b is a quadratic factor of p(x). The roots of x^2 + ax + b can then be found by the quadratic formula. Although the idea is complicated, the implementation of this algorithm is simple. We give a MAPLE procedure.
Maple Implementation
NBarstow := proc(f,a1,a0) local i,k, u,n,s1,s2, jac,b0,b1,b2,c0,c1,c2,d1,d2; n := degree(f,’x’); if n < 3 then RETURN(f); fi; s1 := -a1; s2 :=-a0; for k to 80 do
2.13. THE NEWTON–BARSTOW ALGORITHM 77
b1 := coeff(f,’x’,n); c1 := 0; b0 := coeff(f,’x’,n-1) + s1b1; c0 := b1; for i from n-2 by -1 to 0 do b2 := b1; c2 := c1; b1 := b0; c1 := c0; b0 := coeff(f,’x’,i) + s2b2 + s1b1; c0 := b1 + s2c2 + s1c1; od; jac := c1ˆ2 - c0c2; d1 := (b1c1 -b0c2)/jac; d2 := (b0c1 -b1c0)/jac; s1 := evalf(s1-d1); s2 := evalf(s2-d2); if evalf(d1ˆ2 + d2ˆ2) < 10ˆ(-12)(s1ˆ2+s2ˆ2) then RETURN( (’x’)ˆ2-s1’x’ - s2); fi; od; RETURN(FAIL); end; For this procedure, f is given as an expression in the in- determanate x and the initial guess of a quadratic factor is xˆ2 + a1 x + a0. If successful, the procedure returns an ac- tual (numerically approximate) quadratic factor. The Newton-Barstow algorithm has several drawbacks. It does not work at all well in the case of multiple roots and if a polynomial of odd degree has only one real root, this real root is not a root of any real quadratic factor and hence is invisible to the method, although it can be found by deflation. And unlike Newton’s method, this method cannot be adapted to find zeros of complex functions other than polynomials. However for polynomials the drawbacks can be compensated for, for instance one can use the real Newton’s method to find invisible real roots. And if one suspects multiple roots then these multiple roots are also roots of p′(x) so one strategy is to apply the Newton-Barstow method to the lower degree p′(x), and if that still doesn’t work to p”(x) etc. The Newton-Barstow algorithm is a useful alternative to Newton’s Method in the D’Alembert’s Theorem situation, i.e., one wants real quadratic factors rather than roots. See, for example, section 5.5.
Exercise 2.13.1 [10 points] Factor x^5 − 6 x + 3 as a product of real linear and quadratic factors using the Newton-Barstow algorithm.
2.15. POLYNOMIAL INTERPOLATION 79
In some Linear Algebra books it is suggested that to find the eigenvalues of a square matrix one should calculate the characteristic polynomial and find its roots. For those who still insist on this method we will discuss an efficient way to calculate the charac- teristic polynomial in section 4.7. In practice, however, the modern mathematician will use iterative methods to directly calculate the eigenvalues of a matrix. Since a computer is essential in such a calculation, and software for this problem is widely available, the modern mathematician who wishes to find roots often finds that the quickest way (in human time, not computer time) to find all roots of a polynomial is to construct a matrix whose characteristic polynomial is the desired polynomial and then to pop this matrix into the computer and find the eigenvalues. The first part is easy, for given a monic polynomial p(x) = a 0 +a 1 x+a 2 x^2 +· · ·+xn the companion matrix is the n × n matrix
Cp =
−a 0 −a 1 −a 2 −a 3 · · · −an− 1
and det(Cp) = (−1)np(x) (see, for instance, Birkhoff and MacLane for more de- tails.) In either case the roots are the same as those of p(x) so the eigenvalues of the companion matrix are the roots of the polynomial. Having discussed a variety of methods, some of which show how roots should not be found, some of which show how roots should be found, and this last which shows how roots usually are found we have completed our discussion of root finding algorithms.
2.15 Polynomial Interpolation
No discussion of Numerical Analysis and polynomials would be complete without mentioning polynomial interpolation. Since polynomials are the simplest type of func- tions (other than linear functions) it is desirable to approximate functions by polynomi- als. For any suitably differentiable function a good way to get a polynomial aproxima- tion near a point c is to take the first n terms of the Taylor series. However this is only accurate near c, usually, and we would like an approximation that is close in an entire interval. The method used is interpolation. Given a function f (x) and an interval [a, b] one chooses n + 1 points usually evenly spaced in an interval a = x 0 < x 1 < x 2 <
· · · < xn = b and then constructs a polynomial p(x) of degree n so that p(xj ) = f (xj ) for j = 0, 1 , ..., n. One nice feature of this is that it is not necessary to know much about the function f (x) except its values yj = f (xj ) at these n + 1 points. This is particularly useful when f (x) is known only by experimental data, for we can then find a polynomial which agrees with our data and if necessary we can calculate values of the function for values of x inside our interval other than our experimental values of x. Extrapolation is the word used to describe the method of approximating the value of f (x) for x outside the interval by calculating p(x). However closely the interpolating polynomial approximates the value of f (x) inside the interval [a, b] it usually diverges quickly from f (x) outside the interval. Thus while interpolation is extremely accurate, those who extrapolate usually make fools of themselves.
In light of the above discussion our problem is: given n+1 values of x, x 0 , x 1 , x 2 ,... , xn and n + 1 values of y, y 0 , y 1 ,... yn to find a polynomial p(x) of degree n so that p(xj ) = yj for j = 0, 1 ,... , n. The following theorem of Lagrange shows that this is always possible.
Theorem 2.15.1 Given x 0 , x 1 ,... , xn distinct real or complex numbers and y 0 , y 1 ,... , yn arbitrary real or complex numbers (several may be the same) , there is a unique poly- nomial p(x) of degree n or less satisfying p(xj ) = yj for all j_._
Proof: We define the kth Lagrange interpolating polynomial by
Ln,k(x) =
(x − x 0 ) · · · (x − xk− 1 )(x − xk+1) · · · (x − xn) (xk − x 1 ) · · · (xk − xk− 1 )(xk − xk+1) · · · (xk − xn)
noting that these depend only on the xj ’s, not the yj ’s. Ln,k(x) is a polynomial that takes the value 0 at x 0 ,... xk− 1 , xk+1,... , xn but the value 1 at xk. The desired interpolating polynomial is then
p(x) = y 0 Ln, 0 (x) + y 1 Ln, 1 (x) + · · · + ynLn,n(x)
Uniqueness follows from Theorem 1.7.2 as in the proof of 1.7.3.
From a practical standpoint although one can use the Lagrange interpolating poly- nomials in the proof above to calculate the interpolating polynomial however there are better methods (see a Numerical Analysis text for some.) The best way, however, is to treat this as a linear algebra problem. Letting p(x) = a 0 + a 1 x + a 2 x^2 + · · · + anxn^ we
The Maple procedure interp solves the interpolation prob- lem. For example, to find a polynomial p(x) of degree 3 with p(0) = 18, p(1) = 15, p(3) = − 2 and p(5) = 17 one would enter
p:= interp([0,1,3,5],[18,15,-2,17],x);
The least squares problem can be solved using Maple by either using the regression procedure in the stats package or by set- ting up an under-determined system of linear equations and solving by using the leastsqrs procedure in the linalg package.
One can also generalize by setting conditions on the derivatives, second derivatives etc. to find polynomials displaying desired features, such as roots of various multi- plicity. The important thing is that to generate a polynomial of degree n one should have n + 1 conditions. But only n of these can involve derivatives, n − 1 second or higher derivatives, and so on. Also at least one condition must not involve a value being zero, otherwise the zero polynomial will be your solution. One gets a system of linear equations which can be solved by one of Maple’s solving procedures. You should be warned, however, that the coeficients and/or values on the graph of the polynomial may be large or weird. To try to get a nice example with integer coefficients and reasonable size values the Maple lattice procedure may be helpful as in the following example. Given a real matrix with independent rows the lattice procedure attempts to find a basis for the rowspace consisting of “small” vectors which are integer linear combina- tions of the original rows, we will see more applications of this procedure in Chapter
Maple Implementation
Here we wish to find a polynomial p(x) of degree 4 with root of multiplicity 2 at x = 2, a turning point at x = 0 and a turning point at x = − 1 with as small a y value as possible. Thus we want the value of the polynomial at 2 to be zero, and the derivative to be zero at x = 2, 0 and − 1. To avoid the zero polynomial as our solution we ask that the value p(−1) not be zero, but be as small as possible. We actually also want our second derivatives at x = − 1 , 0 and 2 to not be zero, but we will trust to luck that this will happen. The condition that the derivative be zero at x = 0 tells us the linear coefficient must be zero so our polynomial can be written
p(x) = ax^4 + bx^3 + cx^2 + d
2.15. POLYNOMIAL INTERPOLATION 83
Our conditions then give the following equations:
p(2) = 0 16 a + 8b + 4c + d = 0 p′(2) = 0 32 a + 12b + 4c = 0 p(−1) = s a − b + c + d = s p(−1) = 0 − 4 a + 3b − 2 c = 0
Where s is to be a small, in absolute value, non-zero number. We set up a matrix as follows. Note the first 4 columns refer to the variables a, b, c, d and the last 4 are the equations. We multiply those equations where we want to get zero on the right by a large number (here 10) and those that we want small on the right by 1 or a not quite so large number. We enter the followng 4 x 8 matrix in Maple.
A:=[[1,0,0,0,160,320,1,-40], [0,1,0,0,80,120,-1,30], [0,0,1,0,40,40,1,-20], [0,0,0,1,10,0,1,0]]:
Now we apply the procedure lattice which will provide us with ”small” vectors in the rowspace. The output is shown:
lattice(A);
[[0,0,0,1,10,0,1,0], [1,-1,-5,12,0,0,9,30],
[1,-1,-4,8,0,40,6,10], [-2,3,7,-20,0,0,-18,30]]
Unfortunately we don’t get a solution, i.e., a vector with 0 in the 5, and 8th column. So we try again, multipy the 5,6 and 8th equations by a larger number, now 100, we are trying to force these entries, which would now have to be multiples of 100, to actually be zero. Apply lattice again:
A:=[[1,0,0,0,1600,3200,1,-400], [0,1,0,0,800,1200,-1,300], [0,0,1,0,400,400,1,-200], [0,0,0,1,100,0,1,0]]:
lattice(A);
