CH-03 (Pao - Engineering Analysis)

3Roots of Polynomials andTranscendental Equations3.1 INTRODUCTIONIn the preceding chapter, we derive equations which fit a given of data either exactly,or, by using a criterion such as the least-squares method. Once such equations havebeen obtained in the form of y = C(x) when the data are two-dimensional, or, z =S(x,y) when the data are three-dimensional. It is next of common interest to findwhere the curve C(x) intercepts the x-axis, or, where the surface S(x,y) interceptswith the x-y plane. Mathematically, these are the problems of finding the roots ofthe equations C(x) = 0 and S(x,y) = 0, respectively.

The equation to be solved couldbe a polynomial of the form P(x) = a1 + a2x + … + aixi–1 + … + aN + 1xN which isof Nth order, or, a transcendental equation such as C(x) = a1sinx + a2sin2x + a3sin3x.As it is well known, a polynomial of Nth order should have N roots which couldbe real, or, complex conjugate pair if the coefficients of the polynomial are all real.Geometrically speaking, only those real roots really pass the x-axis. For a transcendental equation, there may be infinite many roots. In this chapter, we shall introducecomputational methods for finding the roots of polynomials and transcendentalequations. Beginning with the very primitive approach of incremental and halfinterval searches, the approximate location of a particular root is to be located. Morerefined, systematic methods such as the linear interpolation and Newton-Raphsonmethods are then followed to determine the more precise location of the root.

Aprogram called FindRoot incorporating the four methods is to be presented forinteractive solution of a particular root of a given polynomial or transcendentalequation when the upper and lower bounds of the root are provided.Also discussed is a method called Successive Substitution. A transcendentalequation derived from analysis of a four-bar linkage problem is used to demonstratehow roots are to be found by application of this method. Another transcendentalequation has been derived for the unit-step response analysis of a mechanical vibration system and its roots solved by application of the Newton-Raphson method toillustrate how the design specifications are checked in the time domain.Since the Newton-Raphson method for solving F(x) = 0 which can be a polynomial, or, transcendental equation of one variable is based on the Taylor’s seriesinvolving the derivatives of F(x), it can be extended to the solution of two-equationsF1(x,y) = 0 and F2(x,y) = 0 by application of Taylor’s series involving partial derivatives of both F1 and F2 with respect to x and y.

A program called NewRaphG hasbeen developed for this purpose. Also, this generalized Newton-Raphson methodallows the quadratic factors of a higher order polynomial to be iteratively and continuously extracted and their quadratic roots solved by the so-called Bairstow method.For that, a program called Bairstow is made available for interactive application.© 2001 by CRC Press LLCBoth QuickBASIC and FORTRAN versions for the above-mentioned programsare presented. Both the application of the roots m-file of MATLAB in place of theprogram Bairstow and direct conversion of the program FindRoot into MATLABversion are also presented. The Mathematica’s function NSolve is introduced inplace of the program Bairstow if the user prefers.

Also the linear interpolationmethod used in the program FindRoot has been translated into Mathematicaversion. In fact. Mathematica has its own FindRoot based on the Newton-Raphsonmethod.3.2 ITERATIVE METHODS AND PROGRAM FINDROOTProgram FindRoot is developed for interactive selection of an iterative methodamong the four made available: (1) Incremental Search, (2) Bisection Search, (3)Linear Interpolation, and (4) Newton-Raphson Iteration. Polynomials are oftenencountered in engineering analyses such as the characteristic equations in vibrational and buckling problems. The roots of a polynomial are related to some important physical properties of the systems being analyzed, such as the frequencies ofvibration or buckling loads. A nth degree polynomial can be expressed as:P(x) = a1 + a 2 x + a 3x 2 + … + a n x n −1 + a n +1x nn +1=∑a k x k −1 = 0(1)k =1For n = 1,2,3, there are formulas readily available in standard mathematicalhandbooks1 for finding the roots.

But for large n values, computer methods are thennecessary to help find the roots of a given polynomial. The methods to be discussedhere are simple and direct and are applicable to not only polynomials but alsotranscendental equations such as 5 + 7cosx – cos60° – cos(60° – x) = 0 related to alinkage design problem2 or x = 40000/{1–0.35sec[40(x/107)0.5]} arisen from buckling study of slender rods.3INCREMENTAL SEARCHFor convenience of discussion, let us consider a cubic equation:P(x) = 1 + 2 x + 3x 2 + 4 x3 = 0(2)To find a root of P(x), we first observe that P(x = –∞)<0, P(x = 0) = 1, and P(x =–∞)>0.

This indicates that the P(x) curve must cross the x axis, possibly once or anodd number of times. Also, the curve may remain above the x-axis or cross it an evennumber of times. To further narrowing down the range on the x-axis, in which the rootis located, we can begin to check the sign of P(x) at x = –10 and search toward theorigin using an increment of x equal to 2. That is, we may construct a list such as:© 2001 by CRC Press LLCXP(x)–10––8––6––4––2–0+Since P(x) changes sign from x = –2 to x = 0, this incremental search can becontinued using an increment of x equal to 0.2 and the left bound x = –10 by replacedby x = –2 to obtain:XP(x)–1.8––1.6–• • •––0.8––0.6+The search continues as follows:xP(x)–0.78––0.76–xP(x)–0.618–xP(x)–0.6058+–0.616–• • •–• • •––0.62––0.60+–0.606––0.604+If only three significant figures accuracy is required, then x = 0.606 is the rootand it has taken 23 incremental search steps to arrive at this answer.

If better accuracyis required, the root should then be sought between x = –0.6060 and x = –0.6058.Program FindRoot prepared both is QuickBASIC and FORTRAN has one ofthe options using the above-explained incremental search method, it also has othermethods of finding the roots of polynomials and transcendental equations to beintroduced next.BISECTION SEARCHThe above example of incremental search shows that if we search from left toright of the x-axis for the root of 4x3 + 3x2 + 2x + 1 = 0 between x = –2 and x = 0,it would be longer than if we search from right to left because the root is near x =0.

Rather than using a fixed incremental in the incremental search method, thebisectional method uses the mid-point of the two bounds of x in search of the root.It involves the testing of the signs of the polynomial at the bounds of the root andreplacing the bounds. The two search methods follow the same procedure. So, thebisection method would go as follows:xP(x)–10–xP(x)–0.46875+© 2001 by CRC Press LLC0+–5––2.5––0.546875+–1.25––0.625––0.5859375+–0.3125+–0.6054688+xP(x)–0.6152344––0.6103516––0.6079102–xP(x)–0.6060792––0.6057741+–0.6059266-2.68817E-04–0.6066895–If we require only three significant figures accuracy, then –0.606 can be considered as the root after having taken 18 bisection search steps.LINEAR INTERPOLATIONNotice that both the incremental and bisection search methods make no use ofthe values of the polynomials at the guessed x values. For example, at x = –10 andx = 1, the polynomial P(x) has values equal to –3719 and 1, respectively.

Since P(x =1) has a smaller value than P(x = –10), we would certainly expect the root to becloser to x = 1 than to x = –10. The linear interpolation makes use of the values ofP(x) at the bounds and calculates a new guessing value of the root using the followingformulas derived from the relationship between two similar triangles:(x − x L )[−P(x )] = (xLR− x ) P( x R )(3)where xL and xR are the left and right bounds of the root, which in this case areequal to –10 and 1, respectively. Based on Equation 3 and P(xL) = –3719 and P(xR) =1, we can have x = –0.002688 and P(x) = 0.9946.

Since P(x)>0, we can thereforereplace xR = 1 with xR = 0.002688. Linear interpolation involves the continuous useof Equation 3 and updating of the bounds.NEWTON-RAPHSON ITERATIVE METHODLinear interpolation method uses the value of the function, for which the rootis being sought; Newton-Raphson method goes one step farther by involving withthe derivative of the function as well.

For example, the polynomial P(x) = 4x3 + 3x2+ 2x + 1 = 0 has its first-derivative expression P'(x) = 12x2 + 6x + 2. If we guessthe root of P(x) to be x = xg and P(xg) is not equal to zero, the adjustment of xg,calling ∆x, can be obtained by application of the Taylor’s series:() ( ) [ ( ) ][ ( ) ]P x g + ∆x = P x g + P ′ x g 1! ∆x + P ′′ x g 2! ( ∆x) …2Since the intention is to find an adjustment x which should make P(xg + x)equal to zero and ∆x itself should be small enough to allow higher order of x tobe dropped from the above expression. As a consequence, we can have 0 = P(xg) +P'(xg)x, or( ) P ′( x )∆x = − P x g© 2001 by CRC Press LLCg(4)Equation 4 is to be continuously used to make new guess, (xg)new = xg + x, ofthe root, until P(x = (xg)new) is negligibly small.The major shortcoming of this method is that during the iteration, if the slopeat the guessing point becomes too small, Equation 4 may lead to a very large Dxso that the xg may fall outside the known bounds of the root.

However, this methodhas the advantage of extending the iterative procedure to solving multiple equationsof multiple variables (see program NewRaphG).An interactive program called FindRoot has been developed in both QuickBASIC and FORTRAN languages with all four methods discussed above. User canselect any one of theses methods, edits the equation to be solved, specifies the boundsof the root, and gives the accuracy tolerance for termination of the root finding. Theprograms are listed below along with sample applications.QUICKBASIC VERSIONSample ApplicationAll four methods have been applied for searching the roots of the equationx2–sin(x)–1 = 0 in the intervals (–1,–0.5) and (1,1.5).