Ïðîñìîòðåòü ôàéë â îòäåëüíîì îêíå: f7cf63efbb16cbaad0b0d3bcd0c9a97d.pdf

3

Roots of Polynomials and

Transcendental Equations

3.1 INTRODUCTION

In 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 have

been 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 find

where the curve C(x) intercepts the x-axis, or, where the surface S(x,y) intercepts

with the x-y plane. Mathematically, these are the problems of finding the roots of

the equations C(x) = 0 and S(x,y) = 0, respectively. The equation to be solved could

be a polynomial of the form P(x) = a1 + a2x + … + aixi–1 + … + aN + 1xN which is

of 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 could

be 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 introduce

computational methods for finding the roots of polynomials and transcendental

equations. Beginning with the very primitive approach of incremental and halfinterval searches, the approximate location of a particular root is to be located. More

refined, systematic methods such as the linear interpolation and Newton-Raphson

methods are then followed to determine the more precise location of the root. A

program called FindRoot incorporating the four methods is to be presented for

interactive solution of a particular root of a given polynomial or transcendental

equation when the upper and lower bounds of the root are provided.

Also discussed is a method called Successive Substitution. A transcendental

equation derived from analysis of a four-bar linkage problem is used to demonstrate

how roots are to be found by application of this method. Another transcendental

equation 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 to

illustrate 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 series

involving the derivatives of F(x), it can be extended to the solution of two-equations

F1(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 has

been developed for this purpose. Also, this generalized Newton-Raphson method

allows 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 LLC

Both QuickBASIC and FORTRAN versions for the above-mentioned programs

are presented. Both the application of the roots m-file of MATLAB in place of the

program Bairstow and direct conversion of the program FindRoot into MATLAB

version are also presented. The Mathematica’s function NSolve is introduced in

place of the program Bairstow if the user prefers. Also the linear interpolation

method used in the program FindRoot has been translated into Mathematica

version. In fact. Mathematica has its own FindRoot based on the Newton-Raphson

method.

3.2 ITERATIVE METHODS AND PROGRAM FINDROOT

Program FindRoot is developed for interactive selection of an iterative method

among the four made available: (1) Incremental Search, (2) Bisection Search, (3)

Linear Interpolation, and (4) Newton-Raphson Iteration. Polynomials are often

encountered 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 of

vibration 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 n

n +1

=

?

a k x k ?1 = 0

(1)

k =1

For n = 1,2,3, there are formulas readily available in standard mathematical

handbooks1 for finding the roots. But for large n values, computer methods are then

necessary to help find the roots of a given polynomial. The methods to be discussed

here are simple and direct and are applicable to not only polynomials but also

transcendental equations such as 5 + 7cosx – cos60° – cos(60° – x) = 0 related to a

linkage design problem2 or x = 40000/{1–0.35sec[40(x/107)0.5]} arisen from buckling study of slender rods.3

INCREMENTAL SEARCH

For 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. This indicates that the P(x) curve must cross the x axis, possibly once or an

odd number of times. Also, the curve may remain above the x-axis or cross it an even

number of times. To further narrowing down the range on the x-axis, in which the root

is located, we can begin to check the sign of P(x) at x = –10 and search toward the

origin using an increment of x equal to 2. That is, we may construct a list such as:

© 2001 by CRC Press LLC

X

P(x)

–10

–

–8

–

–6

–

–4

–

–2

–

0

+

Since P(x) changes sign from x = –2 to x = 0, this incremental search can be

continued using an increment of x equal to 0.2 and the left bound x = –10 by replaced

by x = –2 to obtain:

X

P(x)

–1.8

–

–1.6

–

• • •

–

–0.8

–

–0.6

+

The search continues as follows:

x

P(x)

–0.78

–

–0.76

–

x

P(x)

–0.618

–

x

P(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 root

and it has taken 23 incremental search steps to arrive at this answer. If better accuracy

is 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 of

the options using the above-explained incremental search method, it also has other

methods of finding the roots of polynomials and transcendental equations to be

introduced next.

BISECTION SEARCH

The above example of incremental search shows that if we search from left to

right 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, the

bisectional 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 and

replacing the bounds. The two search methods follow the same procedure. So, the

bisection method would go as follows:

x

P(x)

–10

–

x

P(x)

–0.46875

+

© 2001 by CRC Press LLC

0

+

–5

–

–2.5

–

–0.546875

+

–1.25

–

–0.625

–

–0.5859375

+

–0.3125

+

–0.6054688

+

x

P(x)

–0.6152344

–

–0.6103516

–

–0.6079102

–

x

P(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 INTERPOLATION

Notice that both the incremental and bisection search methods make no use of

the values of the polynomials at the guessed x values. For example, at x = –10 and

x = 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 be

closer to x = 1 than to x = –10. The linear interpolation makes use of the values of

P(x) at the bounds and calculates a new guessing value of the root using the following

formulas derived from the relationship between two similar triangles:

(x ? x L )

[?P(x )] = (x

L

R

? x ) P( x R )

(3)

where xL and xR are the left and right bounds of the root, which in this case are

equal 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 therefore

replace xR = 1 with xR = 0.002688. Linear interpolation involves the continuous use

of Equation 3 and updating of the bounds.

NEWTON-RAPHSON ITERATIVE METHOD

Linear interpolation method uses the value of the function, for which the root

is being sought; Newton-Raphson method goes one step farther by involving with

the 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 guess

the 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) …

2

Since 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 to

be 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 LLC

g

(4)

Equation 4 is to be continuously used to make new guess, (xg)new = xg + x, of

the root, until P(x = (xg)new) is negligibly small.

The major shortcoming of this method is that during the iteration, if the slope

at the guessing point becomes too small, Equation 4 may lead to a very large Dx

so that the xg may fall outside the known bounds of the root. However, this method

has the advantage of extending the iterative procedure to solving multiple equations

of 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 can

select any one of theses methods, edits the equation to be solved, specifies the bounds

of the root, and gives the accuracy tolerance for termination of the root finding. The

programs are listed below along with sample applications.

QUICKBASIC VERSION

Sample Application

All four methods have been applied for searching the roots of the equation

x2–sin(x)–1 = 0 in the intervals (–1,–0.5) and (1,1.5). The negative root equal to

–0.63673 was found after 27, 15, 5, and 3 iterations and the positive root equal to 1.4096

after 29, 15, 4, and 4 iterations by the incremental search, bisection search, linear

interpolation, and Newton-Raphson methods, respectively. An accuracy tolerance of

© 2001 by CRC Press LLC

1.E–5 was used for all cases. For solving this transcendental equation, NewtonRaphson therefore is the best method.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

The interactive question-and-answer process in solving the polynomial 4x3 +

3x + 2x + 1 = 0 using the Newton-Raphson method and the subsequent display on

screen of the iteration goes as follows:

2

Notice that in the FORTRAN program FindRoot, the statement functions F(X)

and FP(X) are defined for calculating the values of the given function and its

derivative at a specified X value. Also, a character variable AS is declared through

a CHARACTER*N with N being equal to 1 in this case when AS can have only

one character as opposed to the general case of having N characters.

© 2001 by CRC Press LLC

MATLAB APPLICATION

A FindRoot.m file can be created and added to MATLAB m files for the purpose

of finding a root of a polynomial or transcendental equation. In this file, the four

methods discussed in the FORTRAN or QuickBASIC versions can all be incorporated. Since some methods require that the left and right bounds, xl and xr, be

provided, the m file listed below includes as arguments these bounds along with the

tolerance and the limited number of iterations:

Notice that the equation for which a root is to be found should be defined in a

mile file called FofX.m, and that if the Newton-Raphson method, i.e., option 4, is

to be used, then the first derivative of this equation should also be defined in a m

file called DFDX.m. We next present four examples demonstrating when all four

methods are employed for solving a root of the polynomial F(x) = 4x3 + 3x2 + 2x

+ 1 = 0 between the bounds x = –1 and x = 0 using a tolerance of 10–5. In addition

to FindRoot.m file, two supporting m files for this case are:

© 2001 by CRC Press LLC

The four sample solutions are (some printout have been shortened for saving

spaces:

© 2001 by CRC Press LLC

Notice that incremental search, half-interval search, interpolation, and NewtonRaphson methods take 28, 17, 16, and 4 iterations to arrive at the root x = –0.6058,

respectively. The last method therefore is the best, but is only for this polynomial

and not necessary so for a general case.

Method of Successive Substitution

As a closing remark, another method called successive substitution is sometimes

a simple way of finding a root of a transcendental equation, such as for solving the

angle in a four-bar linkage problem shown in Figure 1. Knowing the lengths LAB,

LBC amd LCD, and the angle of the driving link AB, the angle of the driven link CD,

can be found by guessing an initial value of ?(0) and then continuously upgraded

using the equation:

[

)]

? 1

?

? ( k +1) = cos?1 ?

L AB cos ? ? L CD + cos ? ? ? ( k ) ?

L

? BC

?

(

(5)

where the superscript k serves as an iteration counter set equal to zero initially. For

? changing from 0 to 360°, it is often required in study of such mechanism to find

the change in ?. This is left as a homework for the reader to exercise.

© 2001 by CRC Press LLC

FIGURE 1. Successive substitution sometimes is a simple way of finding a root of a transcendental equation, such as for solving the angle ? in a four-bar linkage problem.

MATHEMATICA APPLICATIONS

To illustrate how Mathematica can be applied to find a root of F(x) = 1 + 2x

+ 3x2 + 4x3 = 0 in the interval x = [xl,xr] = [–1,0], the linear interpolation is used

below but similar arrangements could be made when the incremental, or, bisection

search, or, Newton-Raphson method is selected instead.

Input[1]: = F[x_]: = 1. + 2*x + 3*x^2 + 4*x^3

Input[2]: = xl = –1; xr = 0; fl = F[xl]; fr = F[xr]; fx = fl;

Input[3]: = Print[“xl = “,xl,” xr = “,xr,” F(xl) = “,fl,” F(xr) = “,fr]

Output[3]: = xl = –1 xr = 0 F(xl) = –2. F(xr) = 1.

Input[4]: = (While[Abs[fx]>0.00001, x = (xr*fl-xl*fr)/(fl-fr);fx = F[x];

Print[“x = “,N[x,5],” F(x) = “,N[fx,5]];

If[fx*fl 2.}

The solution is X = 2. The second example is for finding a root near T = 0.5 for

a transcendental equation described inside the first pair of braces:

In[2]: = FindRoot[{Exp[-.2*T]*Sin[T + 1.37] = = 0.5},{T,0.5}]

Out[2] = {T -> 1.09911}

For solving two simultaneous transcendental equations, two examples are presented below. The first is to find one of the intercepts of two ellipses and the second

is to find one of the intercepts of a circle of radius equal to 2 and a sine curve.

In[3]: = (FindRoot[{(X/3)^2 + (Y/4)^2 = = 1,(X/4)^2 + (Y/3)^2 = = 1},

{X,–2}, {Y,2}]

Out[3] = {X -> –2.4, Y -> 2.4}

In[4]: = FindRoot[{x = = Sqrt[4y^2],y = = Sin[2*x]},{x,1.95},{y,–0.6}]

Out[4] = {x -> 1.90272, y -> –0.616155}

© 2001 by CRC Press LLC

3.4 PROGRAM BAIRSTOW — BAIRSTOW’S METHOD FOR

FINDING POLYNOMIAL ROOTS

Program Bairstow is developed for finding the roots of polynomials based on the

Newton-Raphson’s iterative method for two variables (see program NewRaphG).

Let a nth-order polynomial be denoted as:

P(x) = x N + a1x N ?1 + a 2 x N ?2 + … + a N ?2 x 2 + a N ?1x + a N

(1)

Notice that the highest term xN has a coefficient equal to 1; otherwise the entire

equation must be normalized by dividing by that coefficient. The Bairstow’s method

consists of first selecting a trial divider D(x) = x2 + d1x + d2, and to obtain the

quotient Q(x) = xN–2 + q1xN–3 + q2xN- 4 + ••• + qN–4x2 + qN–3x + qN–2 and a remainder

R(x) = r1x + r2. The objective is to continuously adjust the values of d1 and d2 until

both values of r1 and r2 are sufficiently small. It is apparent that both r1 and r2 are

dependent of d1 and d2. Taylor’s series expansions of r1 and r2 can be written as:

r1 (d1 + ?d1, d 2 + ?d 2 ) = r1 (d1, d 2 ) + r1,d1 (d1, d 2 )?d1

+ r1,d 2 (d1, d 2 )?d 2 + …

(2)

and

r2 (d1 + ?d1, d 2 + ?d 2 ) = r2 (d1, d 2 ) + r2,d1 (d1, d 2 )?d1

+ r2,d 2 (d1, d 2 )?d 2 + …

(3)

where

r1,d1 ? ?r1 ?d1, r2,d 2 ? ?r2 ?d 2

and so on.

The adjustments d1 and d2 are to be calculated so as to make the left-hand

side of Equations 2 and 3 both equal to zero and these adjustments are expected to

be small enough (if the guessed values of d1 and d2 values are sufficiently close to

those which make both r1 and r2 equal to zero) so that the second and higher derivative

terms in Equations 2 and 3 can be dropped. This leads to:

r1,d1 (d1, d 2 )?d1 + r1,d 2 (d1, d 2 )?d 2 = ? r1 (d1, d 2 )

(4)

r2,d1 (d1, d 2 )?d1 + r2,d 2 (d1, d 2 )?d 2 = ? r2 (d1, d 2 )

(5)

and

© 2001 by CRC Press LLC

Cramer’s rule can then be applied to obtain d1 and d2 as:

(

) (r

r

? r1,d 2 r2,d1

)

(6)

(

) (r

r

? r1,d 2 r2,d1

)

(7)

?d1 = ? r1r2,d 2 + r2 r1,d 2

1,d1 2 ,d 2

and

?d 2 = ? r1r2,d1 + r2 r1,d1

1,d1 2 ,d 2

where r1, r2, and their partial derivatives are to be evaluated at (d1,d2). Equations 6

and 7 are to be continuously applied to adjust the guessing values of (d1,d2) until

both r1(d1,d2) and r2(d1,d2) are negligibly small.

To calculate the adjustments d1 and d2 based on Equations 6 and 7, we need

to find the partial derivatives ?r1/?d1, ?r1/?d2, ?r2/?d1, and ?r2/?d2. These derivatives

are, however, depend on the d1 and d2, and the coefficients q’s in the quotient Q(x).

This can be shown by actually carried out the division of P(x) by D(x). The results

are:

q 1 = a 1 ? d1

and

q 2 = a 2 ? q1d1 ? d 2

(8,9)

and

q k = a k ? q k ?1d1 ? q k ?2 d 2 ,

for k = 3, 4,…, N ? 2

(10)

It can also be shown that the coefficients in the remainder R(x) are:

r1 = a N ?1 ? q N ?2 d1 ? q N ?3d 2

and

r2 = a N ? q N ?2 d 2

(11,12)

We notice that Equations 11 and 12 can be included in Equation 10 if k is

extended to N and if the remainder is redefined as:

R(x) = (x + d1 )q N ?1 + q N

(13)

That is, r1 is renamed as qN–1 and r2 is equal to d1qN–1 + qN. As a consequence,

we need to replace r1 and r2 in Equations 6 and 7 by qN–1 and qN. For calculation of

the adjustments d1 and d2, Equation 10 should be used for qN–1 and qN and to

derive their partial derivatives respect to d1 and d2. Since all q’s are functions of d1

and d2, to derive the partial derivatives of the last two q’s we must find the partial

derivatives for all q’s starting with q1. From Equations 8 to 10, we can have:

?q1 ?d1 = ?1,

?q1 ?d 2 = 0,

?q 2 ?d1 = ( ?? q1 ?d1 ) d1 ? q1 = d1 ? q1,

© 2001 by CRC Press LLC

(14,15)

(16)

?q 2 ?d2 = ?(?q 1 ?d2 )q 1 ? 1 = ?

(17)

?q k ?d1 = ?(?q k ?1 ?d1 )d1 ? q k ?1 ? (?q k ?2 ?d1 ) d 2

(18)

?q k +1 ?d 2 = ?(?q k ?d 2 )d1 ? q k ?1 ? (?q k ?1 ?d 2 )d 2

(19)

and for k = 3,4,…,N

It can be concluded from the above results that:

?q k +1 ?d 2 = ?q k ?d1

for k = 1, 2,…, N ? 1

(20)

Now, we can summarize the procedure of Bairstow’s method for factorizing a

quadratic equation from an Nth-order polynomial as follows: (Some changes of

variables are made in the computer programs to be presented next, such as q’s are

changed to b’s, d1 and d2 are changed to u and v, respectively, and c’s are introduced

to represent the derivatives of q’s.)

(1) Specify the values of N, a1 through aN, and a tolerance .

(2) Assume an initial guessing values for d1 and d2 for the divider D(x).

(3) Calculate the coefficients q1 through qN–2 for the quotient Q(x) using

Equations 8 to 10.

(4) Also use Equation 10 to calculate the coefficients qN–1 and qN for the

remainder R(x).

(5) Test the absolute values of qN–1 and qN. If they are both less than , two

root of P(x) are to be calculated by use of the quadratic formulas. The

order of P(x), N, is to be reduced by 2, and q1 through qN–2 are to become

a1 through aN–2, respectively, and return to Step 2. This looping continues

until the quotient Q(x) is of order two or one, for which the root(s) easily

can be calculated.

(6) If the absolute value of either qN–1 or qN is greater than ?, calculate the

partial derivatives of qk with respect to d1, c’s using Equations 14, 16, and

18 for k = 3,4,…,N. The derivatives of qk with respect to d2 are already

available due to Equation 20.

(7) Use Equations 6 and 7 to calculate the adjustments d1 and d2, noticing

that r1 and r2 are to be replaced by qN–1 and qN, respectively. The iteration

is resumed by returning to Step 3.

Both QuickBASIC and FORTRAN versions of the program Bairstow coded

following the steps described above are to be presented next.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

As an example, the polynomial P(x) = x4–5x3 + 13x2–19x + 10 = 0 is solved by

application of the QuickBASIC version of the program Bairstow. The response on

screen is:

The quotient in this case is a quadratic equation:

Q(x) = [x ? (1 + 2i)][x ? (1 ? 2i)] = (x ? 1) ? (2i) = x 2 ? 2 x + 5 = 0.

2

FORTRAN VERSION

© 2001 by CRC Press LLC

2

© 2001 by CRC Press LLC

Sample Application

Consider the polynomial P(x) = x3 + 2x2 + 3x + 4 = 0. When the FORTRAN

version of the program Bairstow is run, the response on screen is:

When the ITERATIONS column indicates 1, it signals that the quotient is of

order one or two. In this case, the quotient is Q(x) = x + 1.65063. In fact, no iteration

has been performed for solving Q(x).

MATLAB APPLICATION

MATLAB has a file called roots.m which can be applied to find the roots of a

polynomial p(x) = 0. To do so, the coefficients of an nth-order p(x) should be ordered

in descending powers of x into a row matrix of order n + 1. For example, to solve

p(x) = x3 + 2x2 + 3x + 4 = 0, we enter:

>> p = [1,2,3,4]; x = roots(p)

and obtain a screen display

As a second example of solving x4–5x3 + 13x2–19x + 10 = 0, MATLAB interactive entries indicated by the leading >> signs and the resulting display are:

© 2001 by CRC Press LLC

Comparing the two examples, we notice that by placing ; after a statement

suppresses the display of the computed value(s). The elements of the first p matrix

(a single row) is not displayed!

It is of interest to introduce the plot capability of MATLAB by use of the results

presented above which involve a polynomial P(x) and its roots. From graphical

viewpoint, the roots are where the polynomial curve crossing the x-axis. MATLAB

has a plot.m file which can be readily applied here. Let us again consider the

polynomial P(x) = x4–5x3 + 13x2–19x + 10 = 0 and plot it for 0?x?3. For adequate

smoothness of the curve, an increment of x equal to 0.1 can be selected for plotting.

The interactive MATLAB commands entered for obtaining Figure 4 are:

>> p=[1,–5,13,–19,10]; x=[0:0.1:3]; y=polyval (p,x);

>> plot(x,y), hold

>> XL=[0 3]; YL=[0 0]; plot(XL,YL)

Notice that another m file polyval of MATLAB has been employed above. The

statement y = polyval(p,x) generates a array of y values using the polynomial defined

by the coefficient vector p and calculated at the values specified in the x array. The

hold statement put the current plot “on hold” so that an additional horizontal line

connecting the two points defined in the XL and YL arrays can be superimposed.

The first plot statement draws the curve and axes and tic marks while the second

plot statement draws the horizontal line.

The horizontal line drawn at y = 0 help to show the intercepts of the polynomial

curve and the x-axis, by observation near x = 1 and x = 2 which confirm the result

found by the MATLAB file roots.m.

MATHEMATICA APPLICATIONS

For finding the polynomial roots, Mathematica’s function NSolve can be

applied readily.

Keyboard input (and then press shift and Enter keys)

NSolve[x^3 + 2x^2 + 3x + 4 = = 0,x]

The Mathematica response is:

Input[1]: =

NSolve[x^3 + 2x^2 + 3x + 4 = = 0,x]

© 2001 by CRC Press LLC

FIGURE 4.

Output[1] =

{{x -> –1.65063}, {x -> –0.174685 — 1.54687 I},

{x -> –0.174685 — 1.54687 I}}

Keyboard input (and then press Shift and Enter keys)

NSolve[x^4–5x^3 + 13x^2–19x + 10 = = 0,x]

The Mathematica response is:

Input[2]: =

NSolve[x^4–5x^3 + 13x^2–19x + 10 = = 0,x]

Output[2] =

{{x -> 1. – 2. I}, {x -> –1 + 2. I}, {x -> 1.}, {x -> 2.}}

To show the locations of the roots of a polynomial, Mathematica’s function

Plot can be applied to draw the polynomial. The following statements (Keyboard

input will hereon be omitted since it is always repeated in the Input response) enable

Figure 5 to be generated:

© 2001 by CRC Press LLC

Input[3]: =

Plot[x^4–5x^3 + 13x^2–19x + 10,{x,0,3},

Frame->True, AspectRatio->1]

Output[3] =

FIGURE 5.

Notice that {x,0,3} specifies the range of x for plotting, Frame->True requests that

the plot be framed, and AspectRatio-> requests that the scales in horizontal and vertical

directions be equal. The graph clearly shows that there are two roots at x = 1 and x = 2.

3.5 PROBLEMS

FINDROOT

1. A root of F(x) = 3x–2e0.5x = 0 is known to exist between x = 1 and x = 2.

Calculate the guessed locations of this root twice by application of the

linear interpolation method.

2. A root is known to exist between x = 0 and x = 1 for the polynomial

P(x) = x3–4.5x2 + 5.75x–1.875 = 0 because P(x = 0) = –1.875 and P(x =

1) = 3.75. What will be the next two guessed root values if linear interpolation method is used? Show details of your calculation.

3. A root is known to exist between x = 1 and x = 2 for the polynomial x3

+ 0.5x2 + 3x–9 = 0.

Based on the linear interpolation method, make two successive guesses

of the location of the root. Show details of the calculations.

© 2001 by CRC Press LLC

4. For finding a root of the polynomial x3–8.9x2–21.94x + 128.576 = 0 within

the bounds x = 0 and x = 4, the linear interpolation method is to be

applied. Show only the details involved in computation of two successive

trial guesses of the root.

5. Use the Newton-Raphson iterative method to find the root of 2X3–5 = 0

between X = 1 and X = 2.

6. Complex roots of a polynomial can be calculated by application of the

program FindRoot simply by treating the variable X in the polynomial

F(X) as a complex variable. Using a complex number which has a real

part and an imaginary part as an initial guess for X to evaluate F(X) and

its derivatives, both values will also be complex. The Newton-Raphson

iterative process is to be continued until both the real and imaginary parts

of F(X) are sufficiently small. According to this outline, modify program

FindRoot to generate a new program NewRaphC for determining a

complex root for the polynomial X4 + 5X2 + 4 = 0.

7. In solving eigenvalue problems (see programs CharacEq and EigenODE), the characteristic equation of an engineering system is in the form

of a polynomial. Physically, the roots of this polynomial may have the

meaning of frequency, or, buckling load, or others. In the program EigenODE, a vibrational problem leads to a characteristic equation of 3–50

2 + 600 – 1000 = 0. Apply the program FindRoot to find a root between

? equal to 1 and 2 accurate to three significant figures. This root represents

the lowest frequency squared.

8. Apply the Newton-Raphson method to find a root of the polynomial f(x) =

3x3 + 2x2–x–30 = 0 by first guessing it to be equal to 3.0. Carry out two

iterative steps by hand calculation to obtain the adjustments that need to

be made in guessing the value of this root.

9. Apply the program FindRoot to solve Problem 8 given above.

10. Apply the linear interpolation method to find a root of the polynomial

f(x) = 3x3 + 2x2–x–30 = 0 between x = 1 and x = 3. Carry out two iterative

steps by hand calculation to obtain the new bounds.

12. The well known secant formula for column bucking3 relating the average

unit load P/A to the eccentricity ratio ec/r2 is:

{ (

) [

P A = ? max 1 + ec r 2 sec ( L r )( P EA) 1 2 2

]}

where ?max is the proportional limit of the column, L/r is the slenderness

ratio, and E is Young’s modulus of elasticity. Solve the above transcendental equation by using ?max = 620 MPa and E = 190 GPa to find P/A

for ec/r2 = 0.1 and L/r = 20.

13. Solve the friction factor f from the Colebrook and White equation6 for

the flow in a pipe (1/f)1/2 = 1.74–0.868{(2K/D) + [18.7/Re(f)1/2]} where

Re is the Reynold’s number and K/D is the relative roughness parameter.

Plot a curve of f vs. Re, and compare the result with the Moody’s diagram.

14. Find the first five positive solution of the equation XJ0(X)–2J1(X) = 0

where J0 and J1 are the Bessel functions of order 0 and 1, respectively.7

© 2001 by CRC Press LLC

15. Write a program SucceSub for implementing the successive substitution

method and apply it to Equation 5 for solving the angle ? for ? changing

from 0° to 360° in equal increment of 15°.

16. Apply the program SucceSub to solve Problem 1.

17. Apply the program SucceSub to solve Problem 12.

18. Rise time is defined as the time required for the response X(t) to increase

its value from 0.1 to 0.9, referring to Equation 8 and Figure 1. For a

second-order system with a1 = 1, a2 = 0.2 sec–1, a3 = 1 sec–1, and a4 = 1.37

radians, use Equation 8 to calculate the rise time of the response X(t) by

applying the computer program FindRoot.

19. Write a m file for MATLAB and name it FindRoot.m and then apply it

for solving Problem 12.

20. Apply FindRoot.m for solving Problem 14.

21. Apply FindRoot.m for solving Problem 18.

22. Apply Mathematica to solve Problem 12.

23. Apply Mathematica to solve Problem 14.

24. Apply Mathematica to solve Problem 18.

NEWRAPHG

1. Shown below are two ellipses which have been drawn using the equations:

2

© 2001 by CRC Press LLC

2

2

? x ? +? y ? =1

? 28 ? ? 24 ?

and

2

? x? ? ? y? ?

+

=1

? 9.5 ? ? 32 ?

where the coordinate axes x and y are the result of rotating the x and y

axes by a counterclockwise rotation of = 30°. The new coordinates can

be expressed in term of x and y coordinates as:

x? = x cos ? + y sin ?

and

y? = ? x sin ? + y cos ?

Use an appropriate pair of values for xP and yP as initial guesses for

iterative solution of the location of the point P which the two ellipses

intercept in the fourth quadrant by application of the program NewRaphG.

2. The circle described by the equation x2 + y2 = 22 and the sinusoidal curve

described by the equation y = sin2x intercepted at two places as shown

below. This drawing obtained using the MATLAB command

axis(‘square’) actually is having a square border when it is shown on

screen but distorted when it is printed because the printer has a different

aspect ratio. Apply the Newton-Raphson iterative method to find the

intercept of these two curves near (x,y) = (2,- 0.5).

© 2001 by CRC Press LLC

3. Apply program NewRaphG for finding a root near X = 0.4 and Y = 0.6

from the equations SinXSinY + 5X–7Y = –0.77015 and e0.1X-X2Y + 3Y =

2.42627. The solutions should be accurate up to 5 significant digits.

4. If one searches near x = 2 and y = 3 for a root of the equations f(x,y) =

(x–1)y + 2x2(y–1)2–35 and g(x,y) = x3–2x2y + 3xy2 + y3–65 what should

be the adjustments of x and y based on the Newton-Raphson method?

5. Write a MATLAB m file and call it NewRaph2.m and then apply it to

solve Problems 1 to 4.

6. Apply Mathematica to solve Problems 1 to 4.

BAIRSTOW

1. Is x2–x + 1 a factor of 4x4–3x3 + 2x2–x + 5? If not, calculate the adjustments for u and v which are equal to –1 and 1, respectively, based on the

Bairstow’s method.

2. Apply plot.m of MATLAB to obtain a plot of P(x) = x6 + x5–8x4 + 14x3

+ 13x2–111x + 90 = 0 vs. x for –4?x?3 as shown below.

© 2001 by CRC Press LLC

3. Find the roots of P(x) = x6 + x5–8x4 + 14x3 + 13x2–111x + 90 = 0 by

application of the program Bairstow. Are the real roots graphically verified in the plot shown in Problem 2?

4. Apply root.m of MATLAB to find the roots of P(x) = x6 + x5–8x4 + 14x3

+ 13x2 111x + 90 = 0.

5. Apply the program Bairstow to solve the characteristic equations derived

by the program CharacEq: (1) 3–15 2–18 = 0, (2) 3–182 + 109

–222 = 0, and (3) 3–122 + 47 –60 = 0.

6. Apply the program Bairstow to solve for the characteristic equations

derived for the program EigenODE.Stb:

?2 ? 36? + 243 = 0,

?3 ? 96?2 ? 2560? ? 16384 = 0,

?4 ? 200?3 + 13125?2 ? 312500? + 1953125 = 0,

and

?5 ? 360?4 + 466563 ? 2612736?2 + 5.878656x10 7 ? ? 3.62791x108 = 0

7. Apply roots.m of MATLAB to find the roots of the polynomials given

in Problem 5 and then graphically verify the locations of these roots by

plotting the polynomial curves using plot.m.

8. Repeat Problem 7 except for those polynomials given in Problem 6.

9. Apply NSolve of Mathematica to solve Problems 2 to 6.

3.6 REFERENCES

1. M. Abramowitz and I. A. Stegun, editors, “Handbook of Mathematical Functions

with Formulas, Graphs, and Mathematical Tables,” Applied Mathematics Series 55,

National Bureau of Standards, Washington, DC, 1964.

2. Y. C. Pao, Elements of Computer-Aided Design and Manufacturing, CAD/CAM, John

Wiley & Sons, New York, 1984.

3. A. Higdon, E. H. Ohlsen, W. B. Stiles, J. A. Weese, and W. F. Riley, Mechanics of

Materials, John Wiley & Sons, New York, Fourth Edition, 1985.

4. J. F. Traub, Iterative Methods for Solution of Equation, Prentice-Hall, Englewood

Cliffs, NJ, 1964.

5. A. E. Taylor, Advanced Calculus, Ginn & Co., Boston, 1955.

6. W. H. Li, Fluid Mechanics in Water Resources Engineering, Allyn & Bacon, Boston,

1983.

7. C. R. Wylie, Jr., Advanced Engineering Mathematics, McGraw-Hill, New York, 1960.

© 2001 by CRC Press LLC

Roots of Polynomials and

Transcendental Equations

3.1 INTRODUCTION

In 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 have

been 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 find

where the curve C(x) intercepts the x-axis, or, where the surface S(x,y) intercepts

with the x-y plane. Mathematically, these are the problems of finding the roots of

the equations C(x) = 0 and S(x,y) = 0, respectively. The equation to be solved could

be a polynomial of the form P(x) = a1 + a2x + … + aixi–1 + … + aN + 1xN which is

of 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 could

be 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 introduce

computational methods for finding the roots of polynomials and transcendental

equations. Beginning with the very primitive approach of incremental and halfinterval searches, the approximate location of a particular root is to be located. More

refined, systematic methods such as the linear interpolation and Newton-Raphson

methods are then followed to determine the more precise location of the root. A

program called FindRoot incorporating the four methods is to be presented for

interactive solution of a particular root of a given polynomial or transcendental

equation when the upper and lower bounds of the root are provided.

Also discussed is a method called Successive Substitution. A transcendental

equation derived from analysis of a four-bar linkage problem is used to demonstrate

how roots are to be found by application of this method. Another transcendental

equation 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 to

illustrate 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 series

involving the derivatives of F(x), it can be extended to the solution of two-equations

F1(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 has

been developed for this purpose. Also, this generalized Newton-Raphson method

allows 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 LLC

Both QuickBASIC and FORTRAN versions for the above-mentioned programs

are presented. Both the application of the roots m-file of MATLAB in place of the

program Bairstow and direct conversion of the program FindRoot into MATLAB

version are also presented. The Mathematica’s function NSolve is introduced in

place of the program Bairstow if the user prefers. Also the linear interpolation

method used in the program FindRoot has been translated into Mathematica

version. In fact. Mathematica has its own FindRoot based on the Newton-Raphson

method.

3.2 ITERATIVE METHODS AND PROGRAM FINDROOT

Program FindRoot is developed for interactive selection of an iterative method

among the four made available: (1) Incremental Search, (2) Bisection Search, (3)

Linear Interpolation, and (4) Newton-Raphson Iteration. Polynomials are often

encountered 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 of

vibration 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 n

n +1

=

?

a k x k ?1 = 0

(1)

k =1

For n = 1,2,3, there are formulas readily available in standard mathematical

handbooks1 for finding the roots. But for large n values, computer methods are then

necessary to help find the roots of a given polynomial. The methods to be discussed

here are simple and direct and are applicable to not only polynomials but also

transcendental equations such as 5 + 7cosx – cos60° – cos(60° – x) = 0 related to a

linkage design problem2 or x = 40000/{1–0.35sec[40(x/107)0.5]} arisen from buckling study of slender rods.3

INCREMENTAL SEARCH

For 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. This indicates that the P(x) curve must cross the x axis, possibly once or an

odd number of times. Also, the curve may remain above the x-axis or cross it an even

number of times. To further narrowing down the range on the x-axis, in which the root

is located, we can begin to check the sign of P(x) at x = –10 and search toward the

origin using an increment of x equal to 2. That is, we may construct a list such as:

© 2001 by CRC Press LLC

X

P(x)

–10

–

–8

–

–6

–

–4

–

–2

–

0

+

Since P(x) changes sign from x = –2 to x = 0, this incremental search can be

continued using an increment of x equal to 0.2 and the left bound x = –10 by replaced

by x = –2 to obtain:

X

P(x)

–1.8

–

–1.6

–

• • •

–

–0.8

–

–0.6

+

The search continues as follows:

x

P(x)

–0.78

–

–0.76

–

x

P(x)

–0.618

–

x

P(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 root

and it has taken 23 incremental search steps to arrive at this answer. If better accuracy

is 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 of

the options using the above-explained incremental search method, it also has other

methods of finding the roots of polynomials and transcendental equations to be

introduced next.

BISECTION SEARCH

The above example of incremental search shows that if we search from left to

right 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, the

bisectional 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 and

replacing the bounds. The two search methods follow the same procedure. So, the

bisection method would go as follows:

x

P(x)

–10

–

x

P(x)

–0.46875

+

© 2001 by CRC Press LLC

0

+

–5

–

–2.5

–

–0.546875

+

–1.25

–

–0.625

–

–0.5859375

+

–0.3125

+

–0.6054688

+

x

P(x)

–0.6152344

–

–0.6103516

–

–0.6079102

–

x

P(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 INTERPOLATION

Notice that both the incremental and bisection search methods make no use of

the values of the polynomials at the guessed x values. For example, at x = –10 and

x = 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 be

closer to x = 1 than to x = –10. The linear interpolation makes use of the values of

P(x) at the bounds and calculates a new guessing value of the root using the following

formulas derived from the relationship between two similar triangles:

(x ? x L )

[?P(x )] = (x

L

R

? x ) P( x R )

(3)

where xL and xR are the left and right bounds of the root, which in this case are

equal 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 therefore

replace xR = 1 with xR = 0.002688. Linear interpolation involves the continuous use

of Equation 3 and updating of the bounds.

NEWTON-RAPHSON ITERATIVE METHOD

Linear interpolation method uses the value of the function, for which the root

is being sought; Newton-Raphson method goes one step farther by involving with

the 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 guess

the 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) …

2

Since 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 to

be 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 LLC

g

(4)

Equation 4 is to be continuously used to make new guess, (xg)new = xg + x, of

the root, until P(x = (xg)new) is negligibly small.

The major shortcoming of this method is that during the iteration, if the slope

at the guessing point becomes too small, Equation 4 may lead to a very large Dx

so that the xg may fall outside the known bounds of the root. However, this method

has the advantage of extending the iterative procedure to solving multiple equations

of 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 can

select any one of theses methods, edits the equation to be solved, specifies the bounds

of the root, and gives the accuracy tolerance for termination of the root finding. The

programs are listed below along with sample applications.

QUICKBASIC VERSION

Sample Application

All four methods have been applied for searching the roots of the equation

x2–sin(x)–1 = 0 in the intervals (–1,–0.5) and (1,1.5). The negative root equal to

–0.63673 was found after 27, 15, 5, and 3 iterations and the positive root equal to 1.4096

after 29, 15, 4, and 4 iterations by the incremental search, bisection search, linear

interpolation, and Newton-Raphson methods, respectively. An accuracy tolerance of

© 2001 by CRC Press LLC

1.E–5 was used for all cases. For solving this transcendental equation, NewtonRaphson therefore is the best method.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

The interactive question-and-answer process in solving the polynomial 4x3 +

3x + 2x + 1 = 0 using the Newton-Raphson method and the subsequent display on

screen of the iteration goes as follows:

2

Notice that in the FORTRAN program FindRoot, the statement functions F(X)

and FP(X) are defined for calculating the values of the given function and its

derivative at a specified X value. Also, a character variable AS is declared through

a CHARACTER*N with N being equal to 1 in this case when AS can have only

one character as opposed to the general case of having N characters.

© 2001 by CRC Press LLC

MATLAB APPLICATION

A FindRoot.m file can be created and added to MATLAB m files for the purpose

of finding a root of a polynomial or transcendental equation. In this file, the four

methods discussed in the FORTRAN or QuickBASIC versions can all be incorporated. Since some methods require that the left and right bounds, xl and xr, be

provided, the m file listed below includes as arguments these bounds along with the

tolerance and the limited number of iterations:

Notice that the equation for which a root is to be found should be defined in a

mile file called FofX.m, and that if the Newton-Raphson method, i.e., option 4, is

to be used, then the first derivative of this equation should also be defined in a m

file called DFDX.m. We next present four examples demonstrating when all four

methods are employed for solving a root of the polynomial F(x) = 4x3 + 3x2 + 2x

+ 1 = 0 between the bounds x = –1 and x = 0 using a tolerance of 10–5. In addition

to FindRoot.m file, two supporting m files for this case are:

© 2001 by CRC Press LLC

The four sample solutions are (some printout have been shortened for saving

spaces:

© 2001 by CRC Press LLC

Notice that incremental search, half-interval search, interpolation, and NewtonRaphson methods take 28, 17, 16, and 4 iterations to arrive at the root x = –0.6058,

respectively. The last method therefore is the best, but is only for this polynomial

and not necessary so for a general case.

Method of Successive Substitution

As a closing remark, another method called successive substitution is sometimes

a simple way of finding a root of a transcendental equation, such as for solving the

angle in a four-bar linkage problem shown in Figure 1. Knowing the lengths LAB,

LBC amd LCD, and the angle of the driving link AB, the angle of the driven link CD,

can be found by guessing an initial value of ?(0) and then continuously upgraded

using the equation:

[

)]

? 1

?

? ( k +1) = cos?1 ?

L AB cos ? ? L CD + cos ? ? ? ( k ) ?

L

? BC

?

(

(5)

where the superscript k serves as an iteration counter set equal to zero initially. For

? changing from 0 to 360°, it is often required in study of such mechanism to find

the change in ?. This is left as a homework for the reader to exercise.

© 2001 by CRC Press LLC

FIGURE 1. Successive substitution sometimes is a simple way of finding a root of a transcendental equation, such as for solving the angle ? in a four-bar linkage problem.

MATHEMATICA APPLICATIONS

To illustrate how Mathematica can be applied to find a root of F(x) = 1 + 2x

+ 3x2 + 4x3 = 0 in the interval x = [xl,xr] = [–1,0], the linear interpolation is used

below but similar arrangements could be made when the incremental, or, bisection

search, or, Newton-Raphson method is selected instead.

Input[1]: = F[x_]: = 1. + 2*x + 3*x^2 + 4*x^3

Input[2]: = xl = –1; xr = 0; fl = F[xl]; fr = F[xr]; fx = fl;

Input[3]: = Print[“xl = “,xl,” xr = “,xr,” F(xl) = “,fl,” F(xr) = “,fr]

Output[3]: = xl = –1 xr = 0 F(xl) = –2. F(xr) = 1.

Input[4]: = (While[Abs[fx]>0.00001, x = (xr*fl-xl*fr)/(fl-fr);fx = F[x];

Print[“x = “,N[x,5],” F(x) = “,N[fx,5]];

If[fx*fl 2.}

The solution is X = 2. The second example is for finding a root near T = 0.5 for

a transcendental equation described inside the first pair of braces:

In[2]: = FindRoot[{Exp[-.2*T]*Sin[T + 1.37] = = 0.5},{T,0.5}]

Out[2] = {T -> 1.09911}

For solving two simultaneous transcendental equations, two examples are presented below. The first is to find one of the intercepts of two ellipses and the second

is to find one of the intercepts of a circle of radius equal to 2 and a sine curve.

In[3]: = (FindRoot[{(X/3)^2 + (Y/4)^2 = = 1,(X/4)^2 + (Y/3)^2 = = 1},

{X,–2}, {Y,2}]

Out[3] = {X -> –2.4, Y -> 2.4}

In[4]: = FindRoot[{x = = Sqrt[4y^2],y = = Sin[2*x]},{x,1.95},{y,–0.6}]

Out[4] = {x -> 1.90272, y -> –0.616155}

© 2001 by CRC Press LLC

3.4 PROGRAM BAIRSTOW — BAIRSTOW’S METHOD FOR

FINDING POLYNOMIAL ROOTS

Program Bairstow is developed for finding the roots of polynomials based on the

Newton-Raphson’s iterative method for two variables (see program NewRaphG).

Let a nth-order polynomial be denoted as:

P(x) = x N + a1x N ?1 + a 2 x N ?2 + … + a N ?2 x 2 + a N ?1x + a N

(1)

Notice that the highest term xN has a coefficient equal to 1; otherwise the entire

equation must be normalized by dividing by that coefficient. The Bairstow’s method

consists of first selecting a trial divider D(x) = x2 + d1x + d2, and to obtain the

quotient Q(x) = xN–2 + q1xN–3 + q2xN- 4 + ••• + qN–4x2 + qN–3x + qN–2 and a remainder

R(x) = r1x + r2. The objective is to continuously adjust the values of d1 and d2 until

both values of r1 and r2 are sufficiently small. It is apparent that both r1 and r2 are

dependent of d1 and d2. Taylor’s series expansions of r1 and r2 can be written as:

r1 (d1 + ?d1, d 2 + ?d 2 ) = r1 (d1, d 2 ) + r1,d1 (d1, d 2 )?d1

+ r1,d 2 (d1, d 2 )?d 2 + …

(2)

and

r2 (d1 + ?d1, d 2 + ?d 2 ) = r2 (d1, d 2 ) + r2,d1 (d1, d 2 )?d1

+ r2,d 2 (d1, d 2 )?d 2 + …

(3)

where

r1,d1 ? ?r1 ?d1, r2,d 2 ? ?r2 ?d 2

and so on.

The adjustments d1 and d2 are to be calculated so as to make the left-hand

side of Equations 2 and 3 both equal to zero and these adjustments are expected to

be small enough (if the guessed values of d1 and d2 values are sufficiently close to

those which make both r1 and r2 equal to zero) so that the second and higher derivative

terms in Equations 2 and 3 can be dropped. This leads to:

r1,d1 (d1, d 2 )?d1 + r1,d 2 (d1, d 2 )?d 2 = ? r1 (d1, d 2 )

(4)

r2,d1 (d1, d 2 )?d1 + r2,d 2 (d1, d 2 )?d 2 = ? r2 (d1, d 2 )

(5)

and

© 2001 by CRC Press LLC

Cramer’s rule can then be applied to obtain d1 and d2 as:

(

) (r

r

? r1,d 2 r2,d1

)

(6)

(

) (r

r

? r1,d 2 r2,d1

)

(7)

?d1 = ? r1r2,d 2 + r2 r1,d 2

1,d1 2 ,d 2

and

?d 2 = ? r1r2,d1 + r2 r1,d1

1,d1 2 ,d 2

where r1, r2, and their partial derivatives are to be evaluated at (d1,d2). Equations 6

and 7 are to be continuously applied to adjust the guessing values of (d1,d2) until

both r1(d1,d2) and r2(d1,d2) are negligibly small.

To calculate the adjustments d1 and d2 based on Equations 6 and 7, we need

to find the partial derivatives ?r1/?d1, ?r1/?d2, ?r2/?d1, and ?r2/?d2. These derivatives

are, however, depend on the d1 and d2, and the coefficients q’s in the quotient Q(x).

This can be shown by actually carried out the division of P(x) by D(x). The results

are:

q 1 = a 1 ? d1

and

q 2 = a 2 ? q1d1 ? d 2

(8,9)

and

q k = a k ? q k ?1d1 ? q k ?2 d 2 ,

for k = 3, 4,…, N ? 2

(10)

It can also be shown that the coefficients in the remainder R(x) are:

r1 = a N ?1 ? q N ?2 d1 ? q N ?3d 2

and

r2 = a N ? q N ?2 d 2

(11,12)

We notice that Equations 11 and 12 can be included in Equation 10 if k is

extended to N and if the remainder is redefined as:

R(x) = (x + d1 )q N ?1 + q N

(13)

That is, r1 is renamed as qN–1 and r2 is equal to d1qN–1 + qN. As a consequence,

we need to replace r1 and r2 in Equations 6 and 7 by qN–1 and qN. For calculation of

the adjustments d1 and d2, Equation 10 should be used for qN–1 and qN and to

derive their partial derivatives respect to d1 and d2. Since all q’s are functions of d1

and d2, to derive the partial derivatives of the last two q’s we must find the partial

derivatives for all q’s starting with q1. From Equations 8 to 10, we can have:

?q1 ?d1 = ?1,

?q1 ?d 2 = 0,

?q 2 ?d1 = ( ?? q1 ?d1 ) d1 ? q1 = d1 ? q1,

© 2001 by CRC Press LLC

(14,15)

(16)

?q 2 ?d2 = ?(?q 1 ?d2 )q 1 ? 1 = ?

(17)

?q k ?d1 = ?(?q k ?1 ?d1 )d1 ? q k ?1 ? (?q k ?2 ?d1 ) d 2

(18)

?q k +1 ?d 2 = ?(?q k ?d 2 )d1 ? q k ?1 ? (?q k ?1 ?d 2 )d 2

(19)

and for k = 3,4,…,N

It can be concluded from the above results that:

?q k +1 ?d 2 = ?q k ?d1

for k = 1, 2,…, N ? 1

(20)

Now, we can summarize the procedure of Bairstow’s method for factorizing a

quadratic equation from an Nth-order polynomial as follows: (Some changes of

variables are made in the computer programs to be presented next, such as q’s are

changed to b’s, d1 and d2 are changed to u and v, respectively, and c’s are introduced

to represent the derivatives of q’s.)

(1) Specify the values of N, a1 through aN, and a tolerance .

(2) Assume an initial guessing values for d1 and d2 for the divider D(x).

(3) Calculate the coefficients q1 through qN–2 for the quotient Q(x) using

Equations 8 to 10.

(4) Also use Equation 10 to calculate the coefficients qN–1 and qN for the

remainder R(x).

(5) Test the absolute values of qN–1 and qN. If they are both less than , two

root of P(x) are to be calculated by use of the quadratic formulas. The

order of P(x), N, is to be reduced by 2, and q1 through qN–2 are to become

a1 through aN–2, respectively, and return to Step 2. This looping continues

until the quotient Q(x) is of order two or one, for which the root(s) easily

can be calculated.

(6) If the absolute value of either qN–1 or qN is greater than ?, calculate the

partial derivatives of qk with respect to d1, c’s using Equations 14, 16, and

18 for k = 3,4,…,N. The derivatives of qk with respect to d2 are already

available due to Equation 20.

(7) Use Equations 6 and 7 to calculate the adjustments d1 and d2, noticing

that r1 and r2 are to be replaced by qN–1 and qN, respectively. The iteration

is resumed by returning to Step 3.

Both QuickBASIC and FORTRAN versions of the program Bairstow coded

following the steps described above are to be presented next.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

As an example, the polynomial P(x) = x4–5x3 + 13x2–19x + 10 = 0 is solved by

application of the QuickBASIC version of the program Bairstow. The response on

screen is:

The quotient in this case is a quadratic equation:

Q(x) = [x ? (1 + 2i)][x ? (1 ? 2i)] = (x ? 1) ? (2i) = x 2 ? 2 x + 5 = 0.

2

FORTRAN VERSION

© 2001 by CRC Press LLC

2

© 2001 by CRC Press LLC

Sample Application

Consider the polynomial P(x) = x3 + 2x2 + 3x + 4 = 0. When the FORTRAN

version of the program Bairstow is run, the response on screen is:

When the ITERATIONS column indicates 1, it signals that the quotient is of

order one or two. In this case, the quotient is Q(x) = x + 1.65063. In fact, no iteration

has been performed for solving Q(x).

MATLAB APPLICATION

MATLAB has a file called roots.m which can be applied to find the roots of a

polynomial p(x) = 0. To do so, the coefficients of an nth-order p(x) should be ordered

in descending powers of x into a row matrix of order n + 1. For example, to solve

p(x) = x3 + 2x2 + 3x + 4 = 0, we enter:

>> p = [1,2,3,4]; x = roots(p)

and obtain a screen display

As a second example of solving x4–5x3 + 13x2–19x + 10 = 0, MATLAB interactive entries indicated by the leading >> signs and the resulting display are:

© 2001 by CRC Press LLC

Comparing the two examples, we notice that by placing ; after a statement

suppresses the display of the computed value(s). The elements of the first p matrix

(a single row) is not displayed!

It is of interest to introduce the plot capability of MATLAB by use of the results

presented above which involve a polynomial P(x) and its roots. From graphical

viewpoint, the roots are where the polynomial curve crossing the x-axis. MATLAB

has a plot.m file which can be readily applied here. Let us again consider the

polynomial P(x) = x4–5x3 + 13x2–19x + 10 = 0 and plot it for 0?x?3. For adequate

smoothness of the curve, an increment of x equal to 0.1 can be selected for plotting.

The interactive MATLAB commands entered for obtaining Figure 4 are:

>> p=[1,–5,13,–19,10]; x=[0:0.1:3]; y=polyval (p,x);

>> plot(x,y), hold

>> XL=[0 3]; YL=[0 0]; plot(XL,YL)

Notice that another m file polyval of MATLAB has been employed above. The

statement y = polyval(p,x) generates a array of y values using the polynomial defined

by the coefficient vector p and calculated at the values specified in the x array. The

hold statement put the current plot “on hold” so that an additional horizontal line

connecting the two points defined in the XL and YL arrays can be superimposed.

The first plot statement draws the curve and axes and tic marks while the second

plot statement draws the horizontal line.

The horizontal line drawn at y = 0 help to show the intercepts of the polynomial

curve and the x-axis, by observation near x = 1 and x = 2 which confirm the result

found by the MATLAB file roots.m.

MATHEMATICA APPLICATIONS

For finding the polynomial roots, Mathematica’s function NSolve can be

applied readily.

Keyboard input (and then press shift and Enter keys)

NSolve[x^3 + 2x^2 + 3x + 4 = = 0,x]

The Mathematica response is:

Input[1]: =

NSolve[x^3 + 2x^2 + 3x + 4 = = 0,x]

© 2001 by CRC Press LLC

FIGURE 4.

Output[1] =

{{x -> –1.65063}, {x -> –0.174685 — 1.54687 I},

{x -> –0.174685 — 1.54687 I}}

Keyboard input (and then press Shift and Enter keys)

NSolve[x^4–5x^3 + 13x^2–19x + 10 = = 0,x]

The Mathematica response is:

Input[2]: =

NSolve[x^4–5x^3 + 13x^2–19x + 10 = = 0,x]

Output[2] =

{{x -> 1. – 2. I}, {x -> –1 + 2. I}, {x -> 1.}, {x -> 2.}}

To show the locations of the roots of a polynomial, Mathematica’s function

Plot can be applied to draw the polynomial. The following statements (Keyboard

input will hereon be omitted since it is always repeated in the Input response) enable

Figure 5 to be generated:

© 2001 by CRC Press LLC

Input[3]: =

Plot[x^4–5x^3 + 13x^2–19x + 10,{x,0,3},

Frame->True, AspectRatio->1]

Output[3] =

FIGURE 5.

Notice that {x,0,3} specifies the range of x for plotting, Frame->True requests that

the plot be framed, and AspectRatio-> requests that the scales in horizontal and vertical

directions be equal. The graph clearly shows that there are two roots at x = 1 and x = 2.

3.5 PROBLEMS

FINDROOT

1. A root of F(x) = 3x–2e0.5x = 0 is known to exist between x = 1 and x = 2.

Calculate the guessed locations of this root twice by application of the

linear interpolation method.

2. A root is known to exist between x = 0 and x = 1 for the polynomial

P(x) = x3–4.5x2 + 5.75x–1.875 = 0 because P(x = 0) = –1.875 and P(x =

1) = 3.75. What will be the next two guessed root values if linear interpolation method is used? Show details of your calculation.

3. A root is known to exist between x = 1 and x = 2 for the polynomial x3

+ 0.5x2 + 3x–9 = 0.

Based on the linear interpolation method, make two successive guesses

of the location of the root. Show details of the calculations.

© 2001 by CRC Press LLC

4. For finding a root of the polynomial x3–8.9x2–21.94x + 128.576 = 0 within

the bounds x = 0 and x = 4, the linear interpolation method is to be

applied. Show only the details involved in computation of two successive

trial guesses of the root.

5. Use the Newton-Raphson iterative method to find the root of 2X3–5 = 0

between X = 1 and X = 2.

6. Complex roots of a polynomial can be calculated by application of the

program FindRoot simply by treating the variable X in the polynomial

F(X) as a complex variable. Using a complex number which has a real

part and an imaginary part as an initial guess for X to evaluate F(X) and

its derivatives, both values will also be complex. The Newton-Raphson

iterative process is to be continued until both the real and imaginary parts

of F(X) are sufficiently small. According to this outline, modify program

FindRoot to generate a new program NewRaphC for determining a

complex root for the polynomial X4 + 5X2 + 4 = 0.

7. In solving eigenvalue problems (see programs CharacEq and EigenODE), the characteristic equation of an engineering system is in the form

of a polynomial. Physically, the roots of this polynomial may have the

meaning of frequency, or, buckling load, or others. In the program EigenODE, a vibrational problem leads to a characteristic equation of 3–50

2 + 600 – 1000 = 0. Apply the program FindRoot to find a root between

? equal to 1 and 2 accurate to three significant figures. This root represents

the lowest frequency squared.

8. Apply the Newton-Raphson method to find a root of the polynomial f(x) =

3x3 + 2x2–x–30 = 0 by first guessing it to be equal to 3.0. Carry out two

iterative steps by hand calculation to obtain the adjustments that need to

be made in guessing the value of this root.

9. Apply the program FindRoot to solve Problem 8 given above.

10. Apply the linear interpolation method to find a root of the polynomial

f(x) = 3x3 + 2x2–x–30 = 0 between x = 1 and x = 3. Carry out two iterative

steps by hand calculation to obtain the new bounds.

12. The well known secant formula for column bucking3 relating the average

unit load P/A to the eccentricity ratio ec/r2 is:

{ (

) [

P A = ? max 1 + ec r 2 sec ( L r )( P EA) 1 2 2

]}

where ?max is the proportional limit of the column, L/r is the slenderness

ratio, and E is Young’s modulus of elasticity. Solve the above transcendental equation by using ?max = 620 MPa and E = 190 GPa to find P/A

for ec/r2 = 0.1 and L/r = 20.

13. Solve the friction factor f from the Colebrook and White equation6 for

the flow in a pipe (1/f)1/2 = 1.74–0.868{(2K/D) + [18.7/Re(f)1/2]} where

Re is the Reynold’s number and K/D is the relative roughness parameter.

Plot a curve of f vs. Re, and compare the result with the Moody’s diagram.

14. Find the first five positive solution of the equation XJ0(X)–2J1(X) = 0

where J0 and J1 are the Bessel functions of order 0 and 1, respectively.7

© 2001 by CRC Press LLC

15. Write a program SucceSub for implementing the successive substitution

method and apply it to Equation 5 for solving the angle ? for ? changing

from 0° to 360° in equal increment of 15°.

16. Apply the program SucceSub to solve Problem 1.

17. Apply the program SucceSub to solve Problem 12.

18. Rise time is defined as the time required for the response X(t) to increase

its value from 0.1 to 0.9, referring to Equation 8 and Figure 1. For a

second-order system with a1 = 1, a2 = 0.2 sec–1, a3 = 1 sec–1, and a4 = 1.37

radians, use Equation 8 to calculate the rise time of the response X(t) by

applying the computer program FindRoot.

19. Write a m file for MATLAB and name it FindRoot.m and then apply it

for solving Problem 12.

20. Apply FindRoot.m for solving Problem 14.

21. Apply FindRoot.m for solving Problem 18.

22. Apply Mathematica to solve Problem 12.

23. Apply Mathematica to solve Problem 14.

24. Apply Mathematica to solve Problem 18.

NEWRAPHG

1. Shown below are two ellipses which have been drawn using the equations:

2

© 2001 by CRC Press LLC

2

2

? x ? +? y ? =1

? 28 ? ? 24 ?

and

2

? x? ? ? y? ?

+

=1

? 9.5 ? ? 32 ?

where the coordinate axes x and y are the result of rotating the x and y

axes by a counterclockwise rotation of = 30°. The new coordinates can

be expressed in term of x and y coordinates as:

x? = x cos ? + y sin ?

and

y? = ? x sin ? + y cos ?

Use an appropriate pair of values for xP and yP as initial guesses for

iterative solution of the location of the point P which the two ellipses

intercept in the fourth quadrant by application of the program NewRaphG.

2. The circle described by the equation x2 + y2 = 22 and the sinusoidal curve

described by the equation y = sin2x intercepted at two places as shown

below. This drawing obtained using the MATLAB command

axis(‘square’) actually is having a square border when it is shown on

screen but distorted when it is printed because the printer has a different

aspect ratio. Apply the Newton-Raphson iterative method to find the

intercept of these two curves near (x,y) = (2,- 0.5).

© 2001 by CRC Press LLC

3. Apply program NewRaphG for finding a root near X = 0.4 and Y = 0.6

from the equations SinXSinY + 5X–7Y = –0.77015 and e0.1X-X2Y + 3Y =

2.42627. The solutions should be accurate up to 5 significant digits.

4. If one searches near x = 2 and y = 3 for a root of the equations f(x,y) =

(x–1)y + 2x2(y–1)2–35 and g(x,y) = x3–2x2y + 3xy2 + y3–65 what should

be the adjustments of x and y based on the Newton-Raphson method?

5. Write a MATLAB m file and call it NewRaph2.m and then apply it to

solve Problems 1 to 4.

6. Apply Mathematica to solve Problems 1 to 4.

BAIRSTOW

1. Is x2–x + 1 a factor of 4x4–3x3 + 2x2–x + 5? If not, calculate the adjustments for u and v which are equal to –1 and 1, respectively, based on the

Bairstow’s method.

2. Apply plot.m of MATLAB to obtain a plot of P(x) = x6 + x5–8x4 + 14x3

+ 13x2–111x + 90 = 0 vs. x for –4?x?3 as shown below.

© 2001 by CRC Press LLC

3. Find the roots of P(x) = x6 + x5–8x4 + 14x3 + 13x2–111x + 90 = 0 by

application of the program Bairstow. Are the real roots graphically verified in the plot shown in Problem 2?

4. Apply root.m of MATLAB to find the roots of P(x) = x6 + x5–8x4 + 14x3

+ 13x2 111x + 90 = 0.

5. Apply the program Bairstow to solve the characteristic equations derived

by the program CharacEq: (1) 3–15 2–18 = 0, (2) 3–182 + 109

–222 = 0, and (3) 3–122 + 47 –60 = 0.

6. Apply the program Bairstow to solve for the characteristic equations

derived for the program EigenODE.Stb:

?2 ? 36? + 243 = 0,

?3 ? 96?2 ? 2560? ? 16384 = 0,

?4 ? 200?3 + 13125?2 ? 312500? + 1953125 = 0,

and

?5 ? 360?4 + 466563 ? 2612736?2 + 5.878656x10 7 ? ? 3.62791x108 = 0

7. Apply roots.m of MATLAB to find the roots of the polynomials given

in Problem 5 and then graphically verify the locations of these roots by

plotting the polynomial curves using plot.m.

8. Repeat Problem 7 except for those polynomials given in Problem 6.

9. Apply NSolve of Mathematica to solve Problems 2 to 6.

3.6 REFERENCES

1. M. Abramowitz and I. A. Stegun, editors, “Handbook of Mathematical Functions

with Formulas, Graphs, and Mathematical Tables,” Applied Mathematics Series 55,

National Bureau of Standards, Washington, DC, 1964.

2. Y. C. Pao, Elements of Computer-Aided Design and Manufacturing, CAD/CAM, John

Wiley & Sons, New York, 1984.

3. A. Higdon, E. H. Ohlsen, W. B. Stiles, J. A. Weese, and W. F. Riley, Mechanics of

Materials, John Wiley & Sons, New York, Fourth Edition, 1985.

4. J. F. Traub, Iterative Methods for Solution of Equation, Prentice-Hall, Englewood

Cliffs, NJ, 1964.

5. A. E. Taylor, Advanced Calculus, Ginn & Co., Boston, 1955.

6. W. H. Li, Fluid Mechanics in Water Resources Engineering, Allyn & Bacon, Boston,

1983.

7. C. R. Wylie, Jr., Advanced Engineering Mathematics, McGraw-Hill, New York, 1960.

© 2001 by CRC Press LLC