VBA homework

2 Solution of Nonlinear Equations

2.1 Introduction

Many engineering analyses require the determination of the value(s) of the vari- able x that satisfies a nonlinear equation:

f .x/ D 0: (2.1) The values of x that satisfy Eq. (2.1) are known as the roots of Eq. (2.1); the num- ber of roots may be finite or infinite and may be real or complex depending on the nature of the equation and the physical problem. The function f (x) may or may not be available in explicit form. If it is available in explicit form, it may be in the form of a general nonlinear equation, or a polynomial or a transcendental equation:

x4 � 80x C 120 D 0:.polynomial equation/ (2.2)

and

tan x � tanh x D 0:.transcendental equation/: (2.3)

C H A P T E R C O N T E N T S

2.1 Introduction 2.2 Engineering Applications 2.3 Incremental Search Method 2.4 Bisection Method 2.5 Newton–Raphson Method 2.6 Secant Method 2.7 Regula Falsi Method 2.8 Fixed Point Iteration or Successive

Substitution Method 2.9 Determination of Multiple Roots 2.10 Bairstow’s Method 2.11 Muller’s Method 2.12 Newton–Raphson Method for

Simultaneous Nonlinear Equations

2.13 Unconstrained Minimization 2.14 Convergence of Methods 2.15 Choice of the Method 2.16 Use of Software Packages 2.17 Computer Programs

48

Section 2.1 Introduction 49

The roots of Eq. (2.1) are also known as the zeros of the function f (x). Thus, a function is said to have zeros, while an equation is said to have roots. A transcendental function is one whose value cannot be determined for any specified value of the argument by a finite number of additions, subtractions, multiplications, or divisions. Exponential, logarithmic, trigonometric, and hyperbolic functions are examples of transcendental functions. Any equation containing transcendental functions is called a transcendental equation. If the function is known only implicitly, we do not have an explicit mathematical expression available for f (x), but a rule will be known to evaluate the value of f (x) for any specified value of the argument. The methods of finding the approximate values of the roots of Eq. (2.1) are considered in this chapter.

2.1.1 Polynomials in x

As a special case of the equation f .x/ D 0, we consider a polynomial equation f .x/ D an xn C an�1xn�1 C Ð Ð Ð C a2x2 C a1x C a0 D 0; (2.4)

where n denotes the degree of the polynomial and a0; a1; a2; : : : ; an are the coeffi- cients of the polynomial that are real (complex, in some cases). Equation (2.4), in general, will have n roots, some or all of which may be complex. The roots may be independent, or some of them may be repeated. If x1; x2; : : : ; xn denote the roots of Eq. (2.4), they are related to the coefficients of the polynomial as [2.8, 2.16]:

nX iD1

xi D �an�1 an

;

nX iD1

nX jD1; j 6Di

xi xj D an�2 an

;

nX iD1

nX jD1; j 6Di

nX kD1;k 6D j

xi xj xk D �an�3 an

;

:::;

x1x2 Ð Ð Ð xn�1xn D .�1/n a0 an

: (2.5)

Equations (2.5) are known as Newton’s relations.

2.1.2 Graphical Interpretation

The real roots of an equation can be interpreted graphically as shown in Fig. 2.1. (Graphical interpretation is not possible if the roots are complex.) The roots occur where the graph of f (x) crosses or touches the x axis. If the graph of f (x) touches or coincides with the x axis at x1, it indicates that the equation f .x/ D 0 has a multiple root at x1. (See Fig. 2.2.)

In some cases, the equation f .x/ D 0 can be split into two separate parts. In such a case, the two parts can be represented graphically, and the point of

50 Chapter 2 Solution of Nonlinear Equations

x1 x2 x3 x4 x

f(x)

f(x)

Figure 2.1 x1; x2; x3; x4 - roots of the equation.

f(x)

f(x) � (1�x)2

f(x) � (1�x)3

(a)

(b)

�1

�1

f(x)

x

x

Figure 2.2 Multiple roots of f .x/ D 0.

intersection of the two parts denotes the root of the equation. As an example, consider the equation

f .x/ D x3 � 2x C 1 D 0; which can be rewritten as

x3 D 2x � 1: (2.6) When the graphs of the functions f1 D x3 and f2 D 2x � 1 are drawn as shown in Fig. 2.3, the points of intersection, x1 and x2, indicate the roots of the equation, f .x/ D 0.

Section 2.1 Introduction 51

x1 x2 x

f1 � x 3

f2 � 2x�1

Figure 2.3 Point of intersection of two graphs.

Uniform beam

M

l Figure 2.4 Bar carrying an end mass.

�

Example 2.1

The natural frequencies of axial vibration of a bar, fixed at one end and carrying a mass M at the other end (see Fig. 2.4), are given by the frequency equation [2.1, 2.2]

cot �

!l

c

� D M

m

!l

c ; (E1)

where ! is the natural frequency, l is the length, c D q

E ²

, E is Young’s modulus, ² is the density, m D ² Al is the mass of the bar, and A is the cross sectional area of the bar. Find the first three natural frequencies of the bar using a graphical procedure for Mm D 0:1, 0.5, and 1.0.

Solution

The frequency equation, Eq. (E1), has the parameters Mm and !l c . The left-hand side

of Eq. (E1), cot �

!l c

Ð , can be plotted as a function of

� !l c

Ð , and the relationship is

shown by the solid lines in Fig. 2.5. The right-hand side of Eq. (E1) represents a straight line in

� !l c

Ð with the slope equal to Mm . It is shown by dotted lines for three

different values of Mm . The intersection points of the solid and dotted lines represent

52 Chapter 2 Solution of Nonlinear Equations

�2

�1

0

1

1.429 4.306

1.0 76

0. 86

0

2

3

4

5

6

7

5 � 2�3 �2��2� �3

3.642

7.229

3.427

6.437

� 6.579c l�

c l�

� 0.1M m

� 0.

5

M m

� 1

.0

M m

� �2

Figure 2.5 Solution of frequency equation.

the roots of the frequency equation. The following results can be obtained from Fig. 2.5.

Frequency Value of M

m

parameter 0.1 0.5 1.0

!1l

c 1.429 1.076 0.860

!2l

c 4.306 3.642 3.427

!3l

c 7.229 6.579 6.437

Note that although the graphical method can be used to find the roots of any equation f .x/ D 0, it has the following limitations:

1. It requires extensive computations to plot the graph of f (x).

2. The roots found from the graphical method will be less accurate.

Section 2.2 Engineering Applications 53

However, in some cases, the roots found from the graphical method can be used as good initial approximations for some of the iterative procedures presented in subsequent sections. �

2.2 Engineering Applications

The solutions of many engineering analysis problems require the determination of the roots of nonlinear equations. The following examples illustrate representative applications.

�

Example 2.2

Water is discharged from a reservoir through a long pipe as shown in Fig. 2.6. By neglecting the change in the level of the reservoir, the transient velocity of the water flowing from the pipe, v(t), can be expressed as [2.4]:

v.t/p 2gh

D tanh �

t

2L

p 2gh

� ; (E1)

where h is the height of the fluid in the reservoir, L is the length of the pipe, g is the acceleration due to gravity, and t is the time elapsed from the beginning of the flow. Find the value of h necessary for achieving a velocity of v D 5 m/sec at time t D 3 sec when L D 5 m and g D 9:81 m/sec2.

Solution

Equation (E1) is rewritten as f1 D f2 with f1 D v.t/p2gh D 5p

2.9:81/h D 1:1288p

h and

f2 D tanh �

3 2.5/

p 2.9:81/h

D tanh.1:3288ph/. The functions f1 and f2 are plotted

by varying h as shown in Fig. 2.7. It can be seen that the two graphs intersect at h D 1:45 m, thereby indicating that Eq. (E1) is satisfied when h D 1:45 m. �

Reservoir

Pipe

v(t)

waterh

L Figure 2.6 Discharge of water from reservoir.

54 Chapter 2 Solution of Nonlinear Equations

10 21.45 3 4 5 h

v

2gh f1 �

1

3

2

4

f2 � tanh t

2L 2gh

Figure 2.7 Solution of Eq. (E1).

�

Example 2.3

The length of a belt in an open-belt drive, L , is given by [2.5]:

L D q

4c2 � .D � d/2 C 1 2

fD�D C d�dg; (E1) where

�D D ³ C 2 sin�1 �

D � d 2c

� (E2)

and

�d D ³ � 2 sin�1 �

D � d 2c

� : (E3)

c is the center distance, D is the diameter of the larger pulley, d is the diameter of the smaller pulley, �D is the angle of contact of the belt with the larger pulley, and �d is the angle of contact of the belt with the smaller pulley. (See Fig. 2.8.) If a belt

c

d DD�d�

Figure 2.8 Open belt drive.

Section 2.2 Engineering Applications 55

having a length 11 m is used to connect two pulleys with diameters 0.4 m and 0.2 m, determine the center distance between the pulleys.

Solution

For the given data, the angles of contact of the belt with the two pulleys can be expressed, from Eqs. (E2) and (E3) as

�D D ³ C 2 sin�1 �

0:4 � 0:2 2c

� D 3:1416 C 2 sin�1

� 0:1 c

� (E4)

and

�d D ³ � 2 sin�1 �

0:4 � 0:2 2c

� D 3:1416 � 2 sin�1

� 0:1 c

� : (E5)

Equation (E1) can be rewritten as

f .c/ D q

4c2 � .D � d/2 C 1 2

fD�D C d�dg � L D 0: (E6) Using the given data and Eqs. (E4) and (E5), Eq. (E6) can be expressed as

f .c/ D q

4c2 � .0:4 � 0:2/2 C 0:5

0:4 ²

3:1416 C 2 sin�1 �

0:1 c

�¦

C0:2 ²

3:1416 � 2 sin�1 �

0:1 c

�¦½ � 11 D 0

or

f .c/ D p

4c2 � 0:04 C 0:24 sin�1 �

0:1 c

� � 10:055752 D 0 (E7)

The function f (c) is shown plotted in Fig. 2.9. It can be seen that f .c/ D 0 at the center distance, c D 5:05 m. �

�

Example 2.4

An electrical circuit consists of a resistor with resistance 14.85 ohms, an inductor with inductance 7 henries, and a capacitor with capacitance 142 farad in parallel as shown in Fig. 2.10(a). If the initial current through the inductance, i(0), and the initial voltage across the capacitor, e(0), are specified as 10 amperes and 0 volts, respectively, the response of the circuit is given by [2.7]:

e.t/ D 210e�1:414t sin 2t: (E1) Determine the time at which the response of the circuit will be zero.

Solution

The response of the circuit, e(t), can be plotted as a function of time, t , as shown in Fig. 2.10(b). It can be seen that the response will be zero at t1 D 1:6 sec, t2 D 3:1 sec, etc. �

56 Chapter 2 Solution of Nonlinear Equations

�10

�9

�8

�7

�6

�5

�4

�3

�2

�1

0 21 3 4 5

5.05

6 7 c

f(c)

1 2 3 4 5 6 7 8

Figure 2.9 Graph of Eq. (E7).

t, sec

e(t) � 210 e�1.414 t sin 2t

e(t), volts

1 �10

0

10

20

30

40

50

60

70

80

2 3

(b)

(a)

4

R � 14.85 ohms

L � 7 henries C �

1 42

farad

i(t)

e(t)

Figure 2.10 Electrical circuit and its response.

Section 2.2 Engineering Applications 57

�

Example 2.5

The shear stress induced along the z axis when two spheres are in contact with each other, while carrying a load F , is given by

h.½/ D 0:75 1 C ½2 C 0:65½ tan

�1 �

1 ½

� � 0:65; (E1)

where h.½/ D −zx

pmax and ½ D z

a ;

in which −zx is the shear force,

pmax D 3F2³a2 (E2) is the maximum pressure developed at the center of the contact area, and the radius of the contact area (see Fig. 2.11) is

a D

8>>< >>:0:34125F

� 1

E1 C 1

E2

� �

1 d1

C 1 d2

� 9>>= >>;

1 3

; (E3)

where E1 and E2 are Young’s moduli of the two spheres, and d1 and d2 are diameters of the two spheres. Poisson’s ratios of the spheres was assumed to be 0.3 in deriving Eqs. (E1) and (E3). In many practical applications such as ball bearings, when the contact load F is large, a crack originates at the point of maximum shear stress and propagates to the surface leading to a fatigue failure. Thus it becomes necessary to

F

F

2a

E2, v2

d2

d1

z

a a

pmax

Contact area

pmax

y

x

E1, v1

Contact area

Figure 2.11 Contact stress between spheres.

58 Chapter 2 Solution of Nonlinear Equations

find the point at which the shear stress attains a maximum value in order to locate the origin of the crack. Determine the value of ½ at which the shear stress, given by Eq. (E1), attains its maximum value.

Solution

To find the value of ½ at which the shear stress attains its maximum value, the condition dhd½ D 0 is used. This yields the equation

dh.½/

d½ D f .½/ D 1:5½

.1 C ½2/2 � 0:65 tan �1 �

1 ½

� C 0:65½

.1 C ½2/ D 0: (E4) This equation is plotted in Fig. 2.12. It can be observed that the root of Eq. (E4) lies at ½ ³ 0:48, which implies that the shear stress will attain its maximum value at z ³ 0:48 a. �

2.3 Incremental Search Method

In the incremental search method, the value of x is incremented from an initial value x1, successively until a change in the sign of the function f (x) is observed. The idea

0.5

0.48

1.0 1.5 2.0 2.5 3.0

0.50

0.25

0

�0.25

�0.50

�0.75

�1.00

�1.25

f( )

Figure 2.12 Graph of Eq. (E4).

Section 2.3 Incremental Search Method 59

is that f (x) changes sign between xi and xiC1, if it has a root in the interval [xi ; xiC1] as shown in Fig. 2.13. This implies that

f .xi / f .xiC1/ < 0; (2.7)

whenever a root is crossed. Here, xiC1 D xi C 1x , where 1x 1x.1/ denotes the initial step size. Once the location of the root between xi and xiC1 is observed, the process can be repeated with a smaller step size 1x .2/ from the new starting point, xi . By successive application of this procedure, we can make the interval 1x .i/, in which the root lies, as small as possible to satisfy the specified convergence requirement.

Notes

1. The incremental search method can be used to find only the real roots of a function f (x). The starting value for the incremental search method, x1, is usually taken as the smallest value in the range of interest. After finding one root as x D xŁ1 , additional roots can be found by using xŁ1 as the new starting point, x1, and the original 1x .1/ as the new step length. The procedure is applied until another change in sign is observed that indicates the crossing of the second root, x D xŁ2 .

2. If the graph of f (x) touches the x axis tangentially, the function f (x) does not undergo a sign change. Hence, the incremental search method cannot find the root.

3. The method may not be able to distinguish between a root and a singular point. (See Fig. 2.14.)

4. An accelerated step size can be used to bracket the root quickly. In this case, the step size is increased (for example, doubled) at each stage.

�

Example 2.6

Find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0; (E1)

using the incremental search method with x1 D 0:0 and 1x .1/ D 0:1.

Solution

Since the initial value is x1 D 0:0 and the initial step size is 1x .1/ D 0:1, the function is evaluated at 0.0, 0.1, 0.2, : : : to obtain the following results:

x1 D 0:0; f1 D f .x1/ D �1:021017;

60 Chapter 2 Solution of Nonlinear Equations

x

Root

x

f(x)

xi xi�1

Figure 2.13 Change in sign of f .x/ at a root.

Root Root

Root

xi xxi �1

��� �2

Singular point

f(x) � cot x

0

Figure 2.14 Roots and singular point.

x2 D 0:1; f2 D f .x2/ D �0:744832; x3 D 0:2; f3 D f .x3/ D �0:490343; x4 D 0:3; f4 D f .x4/ D �0:273916; x5 D 0:4; f5 D f .x5/ D �0:103652; x6 D 0:5; f5 D f .x5/ D C0:020354:

Section 2.4 Bisection Method 61

Since the sign of the function f changes between x5 and x6. f5 f6 < 0/, the step size is reduced to 1x .2/ D 0:01, and the procedure is repeated from x5 D 0:40. This yields the following results:

x5 D 0:40; f5 D f .x5/ D �0:103652; x6 D 0:41; f6 D f .x6/ D �0:089227;

::: ;

x13 D 0:48; f13 D f .x13/ D �0:000956; x14 D 0:49; f14 D f .x14/ D C0:009907:

Since a change in sign of f is observed between x13 and x14. f13 f14 < 0/, the procedure can be continued from x13 with the reduced step size, 1x .3/ D 1x.2/10 D 0:001. However, the procedure is terminated at this stage, and hence, the root of Eq. (E1) can be taken as x13 D 0:48 or x14 D 0:49, or the midpoint of x13 and x14 (0.485). �

2.4 Bisection Method

In order to find the roots of the equation f .x/ D 0 using the bisection method, the function f (x) is first evaluated at equally spaced intervals of x until two successive function values are found with opposite signs. Let a D xk and b D xkC1 be the values of x at which the function values f (a) and f (b) have opposite signs. This implies that the function f (x) has a root between a D xk and b D xkC1 (see Fig. 2.15). The interval .xk; xkC1/, in which the root is expected to lie, is called the interval of uncertainty. The midpoint of the current interval of uncertainty (a; b) is computed as

xmid D a C b2 ; (2.8)

x a � xk

0 b � xk�1

f(xk)

f(x)

f(xk�1)

Figure 2.15 Opposite signs of function at a and b.

62 Chapter 2 Solution of Nonlinear Equations

I (1)

I (3)

I (4)

a (1) a (2) � a (3)

b(3)� b(4) b(1)� b(2)

I (2)

0

f(x)

a (4) x

Figure 2.16 Bisection method.

and the function value f .xmid/ is determined. If f .xmid/ D 0; xmid will be a root of f .x/ D 0. If f .xmid/ 6D 0, then the sign of f .xmid/ will coincide with that of f (a) or f (b). If the signs of f .xmid/ and f (a) coincide, then a is replaced by xmid. Otherwise (that is, if the signs of f .xmid/ and f (b) coincide), b is replaced by xmid. Thus the interval of uncertainty is reduced to half of its original value. Again the midpoint of the current interval of uncertainty is computed using Eq. (2.8), and the procedure is repeated until a specified convergence criterion is satisfied. The reduction of the intervals of uncertainty (i.e., the progress of the iterative process) is shown in Fig. 2.16. The following convergence criterion can be used to stop the iterative procedure:

j f .xmid/j ": (2.9) Here " is a specified small number.

Assuming that the values of a and b, at which f (a) and f (b) have opposite signs, are known, the iterative procedure used to find the roots of f .x/ D 0 can be summarized as follows:

1. Set a.1/ D a; b.1/ D b, and i D 0. 2. Set iteration number i D i C 1. 3. Find xmid D a.i/Cb.i/2 . 4. If xmid satisfies the convergence criterion

j f .xmid/j "; take the desired root as xroot D xmid and stop the procedure. Otherwise, go to step 5.

Section 2.4 Bisection Method 63

5. If f .xmid/ Ð f .a.i// > 0, both f .xmid/ and f .a.i// will have the same sign, hence, set a.iC1/ D xmid and b.iC1/ D b.i/, and go to step 2.

6. If f .xmid/ Ð f .a.i// < 0; f .xmid/ and f .a.i// will have opposite signs, hence, set b.iC1/ D xmid and a.iC1/ D a.i/, and go to step 2.

Notes

1. The bisection method works when the initial interval of uncertainty (a, b) contains an odd number of roots. The method will not work if the interval (a, b) contains a double root (as in Fig. 2.2(a)), since f (a) and f (b) will have the same sign.

2. The method will not be able to distinguish between a singularity and a root. (See Fig. 2.14.)

�

Example 2.7

Find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0 (E1)

using the bisection method with a D 0:0; b D 2:0, and " D 0:005.

Solution

The function values at a and b can be found as f .a/ D f .0:0/ D �1:021017 and f .b/ D f .2:0/ D C0:078629. Using a.1/ D a D 0:0 and b.1/ D b D 2:0, we set the iteration number as i D 1 and find that

xmid D a .1/ C b.1/

2 D 1:0 with f .xmid/ D f .1:0/ D C0:189491:

Since j f .xmid/j > " D 0:005, and f .a.1// and f .xmid/ have opposite signs, we choose a.2/ D a.1/ D 0:0 and b.2/ D xmid D 1:0.

In the next iteration .i D 2/, we use xmid D a.2/Cb.2/2 D 0:5 and find the corresponding function value as f .xmid/ D C0:020353. Since j f .xmid/j > " D 0:005, and f .a.2// and f .xmid/ have opposite signs, we select a.3/ D a.2/ D 0:0 and b.3/ D xmid D 0:5.

In the third iteration .i D 3/, we use xmid D a.3/Cb.3/2 D 0:25 and find f .xmid/ D�0:376661. Since j f .xmid/j > " D 0:005, the process is repeated with i D 4. The results are summarized as follows:

64 Chapter 2 Solution of Nonlinear Equations

Iteration a.i/ b.i/ xmid D a.i/ C b.i/

2 f .xmid/ Is j f .xmid/j < "?

1 0.0 2.0 1.0 0.189491 no

2 0.0 1.0 0.5 0.020353 no

3 0.0 0.5 0.25 �0:376661 no 4 0.25 0.50 0.375 �0:141767 no 5 0.375 0.500 0.4375 �0:051935 no 6 0.4375 0.5000 0.46875 �0:013684 no 7 0.46875 0.50000 0.484375 0.003848 yes

It can be noted that the solution satisfies the convergence criterion, and hence, the root of Eq. (E1) can be taken as x ³ 0:484375. �

2.5 Newton–Raphson Method

The Newton–Raphson method, or simply the Newton’s method, is a well-known and most powerful method used for finding the root of the equation f .x/ D 0. The Newton’s method can be derived by considering the Taylor’s series expansion of the function f (x) about an arbitrary point x1 as

f .x/ D f .x1/ C .x � x1/ f 0.x1/ C 12! .x � x1/ 2 f 00.x1/ C Ð Ð Ð ; (2.10)

where the function, f , and its derivatives, f 0; f 00; : : : on the right-hand side of Eq. (2.10) are evaluated at x1. By considering only the first two terms in the expansion, we have

f .x/ ³ f .x1/ C .x � x1/ f 0.x1/: (2.11) In order to find the root of f .x/ D 0, we set f (x) equal to zero in Eq. (2.11) to obtain

f .x1/ C .x � x1/ f 0.x1/ D 0: (2.12) Since the higher order derivative terms were neglected in the approximation of f (x) in Eq. (2.11), the solution of Eq. (2.12) yields the next approximation to the root (instead of the exact root) as

x D x2 D x1 � f .x1/

f 0.x1/ ; (2.13)

where x2 denotes an improved approximation to the root. To further improve the root, we use x2 in place of x1 on the right-hand side of Eq. (2.13) to obtain x3. This

Section 2.5 Newton–Raphson Method 65

x

f (x)

f (x3)

f (x2)

f �(x3)

f �(x2)

f �(x1) f (x1)

f (x) Exact root

x4 x3 x2 x1 0

Figure 2.17 Newton’s method.

iterative procedure can be generalized as

xiC1 D xi � f .xi /

f 0.xi / : i D 1; 2; Ð Ð Ð : (2.14)

The procedure is shown graphically in Fig. 2.17 assuming a real root for the equation f .x/ D 0. If x1 is the initial guess for the root of f .x/ D 0, the point of intersection of the tangent to the curve at x1 with the x axis gives the next approximation to the root, x2. The convergence of the procedure to the exact root can also be seen in Fig. 2.17.

�

Example 2.8

Find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0 (E1) using Newton–Raphson method with the starting point x1 D 0:0 and the convergence criterion " D 10�5.

Solution

The derivative of f (x) can be found as

d f .x/

dx D 1:5

.1 C x2/2 � 6x2

.1 C x2/3 C 0:65

.1 C x2/ C 0:65

.1 C x2/ � 1:3x2

.1 C x2/2 ;

66 Chapter 2 Solution of Nonlinear Equations

which can be simplified as

f 0.x/ D d f .x/ dx

D 2:8 � 3:2x 2

.1 C x2/3 : (E2) The iterative process of Newton–Raphson method is given by Eq. (2.14):

xiC1 D xi � f .xi /

f 0.xi / : (E3)

At the starting point x1 D 0:0, we have x1 D 0:0; f .x1/ D �1:021017; and f 0.x1/ D 2:800000;

and hence,

x2 D x1 � f .x1/f 0.x1/ D 0:0 � .�1:021017/

2:800000 D 0:364649:

Since f .x2/ D �0:158413 and f 0.x2/ D 1:632745, Eq. (E3) gives x3 D x2 � f .x2/f 0.x2/ D 0:364649 �

.�0:158413/ 1:632745

D 0:461671: The values of f and f 0 at x3 are determined as f .x3/ D �0:021972 and f 0.x3/ D 1:186268, and Eq. (E3) yields

x4 D x3 � f .x3/f 0.x3/ D 0:461671 � .�0:021972/

1:186268 D 0:480193:

At x4, the value of f and f 0 are given by f .x4/ D �0:000742 and f 0.x4/ D 1:106571:

Equation (E3) yields the next approximation as

x5 D x4 � f .x4/f 0.x4/ D 0:480193 � .�0:000742/

1:106571 D 0:480864:

Since the value of f .x5/ D f .0:480864/ D �0:000001 satisfies the convergence criterion j f .x5/j D 0:000001 < 10�5; x5 is taken as the root of Eq. (E1). �

Notes

1. Newton’s method requires the derivative of the function, f 0 D d fdx . In many cases, such as polynomial equations, the derivative of f (x) can be obtained easily. However, in some problems, such as those involving transcendental functions, the differentiation of the function f (x) may be quite complicated. In some cases, the function f (x) may not be available in explicit form. In such cases, the derivative of f (x) can be obtained numerically by using a finite difference approximation. (See Chapter 7.)

2. Newton’s method converges very fast in most cases. However, it may not converge, if the initial guess .x1/ is very far from the exact root. Also, the method may not converge, if the value of the derivative is close to zero or varies substantially near the root. A few instances in which Newton’s method

Section 2.5 Newton–Raphson Method 67

f(x)

f(x)

f(x)

x

(a) f �(x�) � 0 (divergence).

(b) Derivative (slope) at x1 close to zero.

x

x

x� x4x2

x3 x1

x1

x1x2

f �(x2) f �(x1)

(c) Derivative (slope) varies rapidly resulting in oscillation between x1 and x2.

Exact root

Figure 2.18 Non-convergence of Newton’s method.

68 Chapter 2 Solution of Nonlinear Equations

does not converge are illustrated in Fig. 2.18. These difficulties can usually be corrected by starting from a new initial guess point. The following convergence criteria can be used for stopping the iterative procedure of Eq. (2.14):

jxi � xi�1j "; (2.15)

j xi � xi�1 xi

j "; xi 6D 0; (2.16)

and j f .xi /j ": (2.17)

3. Newton’s method can also be used for finding complex roots. In this case, we need to perform the computations using complex arithmetic with a complex number as the initial guess for the root.

�

Example 2.9

Find the root of the equation f .x/ D tan�1 x D 0 using the starting points x1 D 1:6 and x1 D 0:2.

Solution

The sequence of approximations generated by the Newton–Raphson method from the starting point x1 D 1:6 is given by x1 D 1:6; x2 D �2:0034213; x3 D 3:5509133; x4 D �14:090254; x5 D 285:20148; x6 D �127199:55; x7 D 2:5414801 ð 1010; : : : ; which can be seen to be diverging. On the other hand, when the starting point x1 D 0:2 (which is closer to the root xŁ D 0:0) is used, the Newton–Raphson method converged as follows: x1 D 0:2; x2 D �0:00529137; x3 D 0:98720193 ð 10�7, and x4 D 0:0.

This example indicates that, when f 00.xŁ/ D 0, the Newton–Raphson method might diverge, unless the starting point is sufficiently close to the root. �

�

Example 2.10

Find the root of the equation f .x/ D 5x3 � 39x2 C 22x � 60 D 0 using the starting points x1 D 0:95 and x1 D 0:0.

Solution

The Newton–Raphson method produced the following sequence of approximations from the first starting point: x1 D 0:95, x2 D 1:5569093, x3 D �1:0405447, x4 D 0:03627062, x5 D 0:75702357, x6 D 1:3309088, x7 D 3:0455060, x8 D 1:1724654, x9 D 2:0151870; : : : It can be seen that the method never converged. However, the use of the second starting point produced the root, xŁ D 5:0. The convergence can be

Section 2.5 Newton–Raphson Method 69

seen from the results: x1 D 0:0, x2 D 0:73170733, x3 D 1:3058259, x4 D 2:75444658, x5 D 1:4133221, x6 D 5:9946117, x7 D 5:2962012, x8 D 5:0381327, x9 D 5:0007577, and x10 D 5:0000005. �

�

Example 2.11

Find the root of the equation f .x/ D x3 � x �3 D 0 using the starting points x1 D 0:0 and x1 D 0:5.

Solution

The following sequence of approximations were generated by the Newton–Raphson method from x1 D 0:0 : x1 D 0:0, x2 D �3:0, x3 D �1:9615384, x4 D �1:1471759, x5 D �0:00657928, x6 D �3:0003891, x7 D �1:9618182, x8 D �1:1474303, x9 D �0:00725639, x10 D �3:0004730; : : : It can be seen that the method has not con- verged; rather the method resulted in a cyclic process with xiC4 ³ xi .i D 1; 2; 3; : : :/. In such a case, a different starting point is likely to result in the convergence of the process. The use of x1 D 0:5 led to the following sequence of approximations: x1 D 0:5, x2 D �13:0, x3 D �8:6778660, x4 D �5:7976274, x5 D �3:8737497, x6 D �2:5730152, x7 D �1:6472325, x8 D �0:83179486, x9 D 1:7189554, x10 D 1:6731508, x11 D 1:6717013, x12 D 1:6716999, and x13 D 1:6716999. This indicates the root as xŁ D 1:6716999. �

�

Example 2.12

Find the root of the equation f .x/ D sin 2x D 0 using the starting point x1 D 0:75.

Solution

The roots of f (x) nearest to the starting point x1 are 0.0 and ³2 ³ 1:5708. How- ever, the Newton–Raphson method converged to a root far from x1 as follows: x1 D 0:75; x2 D �6:3007; x3 D �6:2832, and x4 D �6:2832 ³ �2³ . This example illustrates that when f 0.x1/ is small (that is, when the tangent to the curve f (x) is nearly horizontal), the Newton–Raphson method might converge to a root farther than the one closer to the starting point. Thus, it is a good practice to verify that the root given by the Newton–Raphson method is the desired one. �

�

Example 2.13

Find the root of the equation f .x/ D x2 C x C5 D 0 using the starting points x1 D 1:0 and x1 D 1 C i .

Solution

When the first starting point was used, the Newton–Raphson method gave x1 D 1:0; x2 D �1:3333; x3 D 1:9333; x4 D �0:2594; x5 D �10:2492; x6 D �5:1310; x7 D

70 Chapter 2 Solution of Nonlinear Equations

�2:3026; : : : with no convergence. The exact roots of the equation are given by �0:5š2:1794i . By using the complex number as a starting point .x1 D 1 C i/, the Newton–Raphson method converged to the correct root, xŁ D �0:5 C 2:179449i , as follows: x1 D 1 C i; x2 D �0:8461539 C 1:230769i; x3 D �0:1701357 C 2:403620i , x4 D �0:4681635 C 2:171638i; x5 D �0:5001113 C 2:179229i , and x6 D �0:5C 2.179449i . �

2.6 Secant Method

The secant method is similar to the Newton’s method, but is different in that the derivative f 0 is approximated by using two consecutive iterative values of f . By using an approximation for the derivative f 0.xi / as

f 0.xi / ³ f .xi / � f .xi�1/ xi � xi�1

; (2.18)

Eq. (2.14) can be rewritten as

xiC1 D xi � f .xi /

f 0.xi / D xi � f .xi /

xi � xi�1

f .xi / � f .xi�1/ ½

; i D 2; 3; 4; : : : : (2.19)

Thus, Eq. (2.19) represents the general expression for the iterative process of the secant method. Notice that two initial guesses x1 and x2 are required to start the iterative process according to Eq. (2.19). The graphical illustration of the secant method and its convergence is shown in Fig. 2.19. It can be noted that geometrically

f (x)

f (x1)

f (x2)

f (x3)

x

Exact root

x1

f (x)

x4 x3 x2

0

Chord

Chord

Figure 2.19 Convergence of secant method.

Section 2.6 Secant Method 71

f(x)

x3 at infinity

x x2 x1

Figure 2.20 Nonconvergence of secant method.

the intersection of the chord joining the points .xi ; f .xi // and .xi�1; f .xi�1// with the x axis is taken as the next approximation in the secant method.

The iterative process can be stopped whenever any or all of the convergence criteria, given by Eqs. (2.15) to (2.17), are satisfied. The process may not converge, if f .x1/ ³ f .x2/ as shown in Fig. 2.20. In such a case, the next approximation .x3/ will be near infinity. Another difficulty might arise when two consecutive values of xi�1 and xi become very close, thereby causing f .xi�1/ and f .xi / also to be very close, and ultimately resulting in a significant round-off error in Eq. (2.19). This difficulty can be avoided by fixing or freezing the values of xi�1 and f .xi�1/ whenever j f .xiC1/j becomes smaller than a prescribed number, ". This means that Eq. (2.19) is modified as

xiCk D xiCk�1 � f .xiCk�1/

xiCk�1 � xi�1 f .xiCk�1/ � f .xi�1/

½ ; k D 2; 3; 4; : : : ; (2.20)

if j f .xiCk/j ". The following iterative process can be used to implement the secant method:

1. Start with two initial approximations x1 and x2 for the root of f .x/ D 0 and a small number " to test the convergence of the process. Set i D 2.

2. Find the new approximation, xiC1, as

xiC1 D xi � f .xi /.xi � xi�1/ f .xi / � f .xi�1/

:

3. Verify the convergence of the process. If

j f .xiC1/j "; stop the process by taking xiC1 as the root. Otherwise, update the iteration number as i D i C 1 and go to step 2.

�

Example 2.14

Find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

.1 C x2/ D 0

using 0.0 and 0.5 as the initial approximations with " D 10�5.

72 Chapter 2 Solution of Nonlinear Equations

Solution

At the initial approximate points, the function values are given by

f .x1/ D f .0:0/ D �1:02101707and f .x2/ D f .0:5/ D 0:02035314: Since j f .0:0/j > " and j f .0:5/j > ", we proceed to find the next approximation using Eq. (2.19):

x3 D x2 � f .x2/.x2 � x1/f .x2/ � f .x1/

D 0:5 � .0:02035314/.0:5 � 0:0/ 0:02035314 � .�1:02101707/ D 0:49022764:

The value of f .x3/ is found to be

j f .x3/j D j f .0:49022764/j D j0:01015013j > "; and hence, we compute the next approximate value as

x4 D x3 � f .x3/.x3 � x2/f .x3/ � f .x2/

D 0:49022764 � .0:01015013/.0:49022764 � 0:5/ 0:01015013 � 0:02035314 D 0:48540229:

The value of f .x4/ is found as

j f .x4/j D j f .0:48540229/j D j0:00496507j > ": Hence, we proceed to the next iteration. The results of the iterative process are as follows:

Value of i Point .xi / f .xi / Is j f .xi /j "? xiC1 given by Eq. (2.19) 1 0.0 �1.02101707 no — 2 0.5 0.02035314 no 0.49022764

3 0.49022764 0.01015013 no 0.48540229

4 0.48540229 0.00496507 no 0.48305339

5 0.48305339 0.00240600 no 0.48191774

6 0.48191774 0.00116086 no 0.48137063

7 0.48137063 0.0055802 no 0.48110777

8 0.48110777 0.00026876 no 0.48098111

9 0.48098111 0.00012833 no 0.48092049

10 0.48092049 0.00006276 no 0.48089087

11 0.48089087 0.00002861 no 0.48087740

12 0.48087740 0.00001460 no 0.48087072

13 0.48087072 0.00000721 yes —

�

Section 2.7 Regula Falsi Method 73

2.7 Regula Falsi Method

The regula falsi method (also known as the method of false position or the method of linear interpolation) is based on linear interpolation. The method is similar to the bisection method, in that the initial interval (a, b) is generated so as to bracket a root of the equation f .x/ D 0. However, instead of continuously bisecting the interval, the linear interpolation of f (x) passing through the two end points of the current interval of uncertainty is used to find a new approximation for the root. Thus, if x1 D a and x2 D b denote two points on either side of the exact root, a linear equation (straight line) is assumed through the points [x1; f .x1/] and [x2; f .x2/] as

y.x/ D px C q; (2.21) where y D f .x1/ at x D x1 and y D f .x2/ at x D x2. This yields

y.x/ D

f .x2/ � f .x1/ x2 � x1

½ x C

² x2 f .x1/ � x1 f .x2/

x2 � x1

¦ : (2.22)

The point of intersection of this straight line, y(x), with the x axis, x3, can be found by setting y(x), given by Eq. (2.22), to zero:

y.x3/ D

f .x2/ � f .x1/ .x2 � x1/

½ x3 C

² x2 f .x1/ � x1 f .x2/

x2 � x1

¦ D 0: (2.23)

This gives the new approximation to the root .x3/ as

x3 D x2 � f .x2/.x2 � x1/f .x2/ � f .x1/ : (2.24)

It can be seen that x3 will be closer to the exact root, xŁ, than either x1 or x2. Equation (2.24) can be generalized as

xiC1 D xi � f .xi /.xi � xi�1/f .xi / � f .xi�1/ : (2.25)

For the next iteration, xiC1 replaces xi if f .xiC1/ f .xi / > 0 or xi�1 if f .xiC1/ f .xi�1/ > 0. The convergence of the process is shown in Fig. 2.21. The iterative process can be stopped whenever the convergence criteria

j f .xi /j " (2.26) and

jxi � xi�1j " (2.27) are satisfied. In these equations, " is a small number on the order of 10�3 to 10�6.

Notes

1. In some cases, the convergence of the regula falsi method may be very slow as indicated in Fig. 2.21. Here one of the end points of interpolation .x1/ does not change from the original end point (such a point is called a stagnation point). Assuming that x1 is the stagnation point and xi is the current approximation

74 Chapter 2 Solution of Nonlinear Equations

f(x)

f(x2)

f(x1)

f(x3) f(x4)

f(x5) Exact root

x x5 x1x4x3

x2 x�

Third interpolation

Second interpolation

First interpolation

Third approximation to root

First approximation to root Second approximation to root

Figure 2.21 Slow convergence of regula falsi method.

f(x2)

f(x1)

f(x)

f(x3) f(xi)

xi�1 x

Stagnation point

x1

xix2 x3

x�

f(xi�1)

Figure 2.22 Improvement of stagnation point. (xiC1 is determined as middle point of x1 and xi ; xiC1 D x1Cxi2 .)

to the root, the convergence of the process can be improved by taking the next value, xiC1, as

� x1Cxi 2

Ð . This procedure is illustrated in Fig. 2.22.

2. Equation (2.25) can be seen to be identical to Eq. (2.19) of the secant method. However, the regula falsi method converges in cases where the secant method might diverge. Since the estimates xi and xiC1 in the regula falsi method are selected such that the values f .xi / and f .xiC1/ have opposite signs, the method converges all the time (Fig. 2.23(a)). In the secant method, on the other hand, the current estimates, xi and xiC1, may not bracket the minimum (as shown in Fig. 2.23(b)). Hence, the method may diverge in some cases.

Section 2.8 Fixed Point Iteration or Successive Substitution Method 75

x1 x

f(x)

f(x)

f(x)

f(x)

x2x3

x�

Iteration 1

x1 x x2x3

x�

Iteration 2

(a) Regula falsi method (converges).

x1 x x2x3

x�

Iteration 1

x x2x3x4 x

�

Iteration 2

(b) Secant method (diverges).

x1

Figure 2.23 Comparison of regula falsi and secant methods.

2.8 Fixed Point Iteration or Successive Substitution Method

In this method, the equation f .x/ D 0 is rewritten in the form x D g.x/; (2.28)

and an iterative procedure is adopted using the relation

xiC1 D g.xi /; i D 1; 2; 3; Ð Ð Ð ; (2.29) where a new approximation to the root, xiC1, is found using the previous approxima- tion, xi (x1 denotes the initial guess). The iterative process can be stopped whenever the convergence criterion

jxiC1 � g.xiC1/j " (2.30) is satisfied, where " is a small number on the order of 10�3 to 10�6. The method is very simple; however, it may not always converge with an arbitrarily chosen form of the function g(x). The condition to be satisfied for convergence to the correct root is given by [2.14]

jg0.x/j < 1 (2.31) in the neighborhood of the correct root. It has been observed that convergence of the process is asymptotic if 0 < g0.x/ < 1 and oscillatory if �1 < g0.x/ < 0.

76 Chapter 2 Solution of Nonlinear Equations

f2 � g(x)

f2 � g(x) f2 � g(x)

f1�x

f2�g(x)

f1�x

f1�x f1�x

x x1 x3 x1x5 x4 x2

x�

(Correct root) (Correct root)

(Correct root) (Correct root)

x2 x3 x� x

(a) �1 g � (x) 0 ; Convergent. (b) 0 g � (x) 1 ; Convergent.

x1 x � x�

x

(c) g � (x) �1 ; Divergent. (d) g � (x) � 1 ; Divergent.

x1x2 x2 x3 x4x3 x

Figure 2.24 Fixed point iteration.

In addition, the convergence will be faster as jg0.x/j ! 0 in the neighborhood of the correct root. The convergence of the process in different cases is illustrated in Fig. 2.24.

�

Example 2.15

Find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0 (E1)

using the fixed point iteration with x1 D 0:0 and " D 10�5.

Section 2.9 Determination of Multiple Roots 77

Solution

Equation (E1) can be rearranged as

x D g.x/; where

g.x/ D 13 30

.1 C x2/2 tan�1 �

1 x

� � 13

30 x.1 C x2/:

Thus, the iterative process can be expressed as

xiC1 D g.xi /: At the starting point, x1 D 0:0, the value of g can be found as g.0:0/ D 0:606535. Since x1 6D g.x1/, we set x2 D g.x1/ D 0:606535 and find g.x2/ D 0:472035. Again x2 6D g.x2/, and hence we proceed to the next iteration with x3 D 0:472035. The results of the iterative process are as follows:

Iteration xi g.xi / Is jxi � g.xi /j " ? number (i)

1 0.100000 0.606535 no

2 0.606535 0.472035 no

3 0.472035 0.481914 no

4 0.481914 0.480744 no

5 0.480744 0.480878 no

6 0.480878 0.480863 no

7 0.480863 0.480864 yes

It can be seen that the specified convergence criterion is satisfied at the end of seventh iteration; as such, the root of Eq. (E1) is taken as x7 D 0:480863. �

2.9 Determination of Multiple Roots

A function f (x) is said to have a root xŁ of multiplicity p > 1 if

f .x/ D .x � xŁ/pg.x/; (2.32) where g(x) is continuous at x D xŁ with g.xŁ/ 6D 0. If p is an integer, Eq. (2.32) gives

f .xŁ/ D f 0.xŁ/ D Ð Ð Ð D f .p�1/.xŁ/ D 0; f .p/.xŁ/ 6D 0: (2.33) The function f (x) does not change sign at xŁ, if xŁ is an even multiple root .p D even number/. Since the function value f (x) and its derivative f 0.x/ are zero at a multiple root xŁ, the Newton–Raphson and secant methods, which contain derivatives in the denominator, cannot be applied conveniently to find the root xŁ.

78 Chapter 2 Solution of Nonlinear Equations

It was shown [2.18] that the Newton–Raphson and secant methods converge only linearly for multiple roots. However, a slightly modified formula, given by Eq. (2.34), has been shown to lead to quadratic convergence,

xiC1 D xi � p f .xi /

f 0.xi / ; (2.34)

where p denotes the multiplicity of the root (p D 2 for double root, p D 3 for triple root, etc.). Equation (2.34) is not very practical since its use requires a prior knowledge of the multiplicity of the root.

Another approach has been suggested [2.18] for modifying Newton–Raphson method to find multiple roots of a function. In this approach, we consider the function

g.x/ D f .x/ f 0.x/

: (2.35)

The function g(x) can be shown to have all its roots at the same locations as those of f (x). Differentiation of g(x) gives

g0.x/ D f 0.x/ f 0.x/ � f .x/ f 00.x/

f f 0.x/g2 : (2.36)

The use of g(x) in place of f (x) in Eq. (2.14) leads to the desired formula:

xiC1 D xi � g.xi /g0.xi / : (2.37)

�

Example 2.16

Find the roots of the equation

f .x/ D .x � 2/2.x � 5/ D x3 � 9x2 C 24x � 20 D 0 (E1) starting from x1 D 0:5 using Eqs. (2.14), (2.34) with p D 2, and (2.37).

Solution

Noting that f 0.x/ D 3x2 � 18x C 24 (E2)

and f 00.x/ D 6x � 18; (E3)

Eqs. (2.14), (2.34), and (2.37) can be expressed as

xiC1 D xi �

x3i � 9x2i C 24xi � 20 3x2i � 18xi C 24

! ; (E4)

xiC1 D xi � 2

x3i � 9x2i C 24xi � 20 3x2i � 18xi C 24

! ; (E5)

Section 2.10 Bairstow’s Method 79

and

xiC1 D xi � (

.x3i � 9x2i C 24xi � 20/.3x2i � 18xi C 24/ .3x2i � 18xi C 24/2 � .x3i � 9x2i C 24xi � 20/.6xi � 18/

)

D xi �

3x5i � 45x4i C 258x3i � 708x2i C 936xi � 480 3x4i � 36x3i C 162x2i � 312xi C 216

! : (E6)

Using x1 D 0:5, the sequence of approximations produced by Eqs. (E4), (E5), and (E6) are shown in the following table.

Iteration Eq. (2.14) or Eq. (2.34) or Eq. (2.37) or

i (E4) (E5) (E6)

1 0.5 0.5 0.5

2 1.1428571 1.7857143 2.1578946

3 1.5285714 1.9930871 2.0044441

4 1.7492974 1.9998698 1.5044441

5 1.8699937 2.0047598 2.0298657

6 1.9336742 — 1.9989179

7 1.9664814 —

8 1.9831409

9 1.9915178

10 1.9958467

11 1.9979852

12 1.9990886

13 1.9994373

14 2.0000024

It can be seen that the modified Newton–Raphson methods, Eqs. (2.34) and (2.37), are more efficient in finding the multiple root compared with the standard New- ton–Raphson method, Eq. (2.14). �

2.10 Bairstow’s Method

Bairstow’s method is an iterative method in which a quadratic factor of the polynomial is found in each stage. When the procedure is applied repeatedly

80 Chapter 2 Solution of Nonlinear Equations

to the deflated polynomials, all the quadratic factors of the polynomial can be determined. Consider an nth degree polynomial equation,

f .x/ D an xn C an�1xn�1 C an�2xn�2 C Ð Ð Ð C a2x2 C a1x C a0: (2.38) By factoring an arbitrary quadratic term, f (x) can be expressed as

f .x/ D .x2 C Þx C þ/.bnxn�2 C bn�1xn�3 C Ð Ð Ð C b4x2 C b3x C b2/ C R.x/; (2.39) where R(x) is the remainder. The remainder term can be found to be of the form

R.x/ D b1x C b0: (2.40) By multiplying the quadratic term and the (n-2)nd degree polynomial in Eq. (2.39) and equating the coefficients of like powers of x , we obtain the following system of equations:

an D bn; an�1 D bn�1 C Þbn; an�2 D bn�2 C Þbn�1 C þbn;

:::;

a2 D b2 C Þb3 C þb4; a1 D b1 C Þb2 C þb3;

and a0 D b0 C þb2: (2.41)

Equations (2.41) can be solved to find the coefficients bn; bn�1; : : : ; b1; b0 as bn D an;

bn�1 D an�1 � Þbn; bn�2 D an�2 � Þbn�1 � þbn;

:::;

b2 D a2 � Þb3 � þb4; b1 D a1 � Þb2 � þb3;

and b0 D a0 � þb2: (2.42)

If we require .x2 C Þx C þ/ to be an exact quadratic factor of f (x), then the remainder, R(x), given by Eq. (2.40) should be zero. This requires that

b0 D 0; b1 D 0: (2.43) It can be seen from Eq. (2.42) that b0 and b1 are functions of Þ and þ. Thus, if NÞ andNþ denote the correct values of Þ and þ, Eq. (2.43) gives

b0. NÞ; Nþ/ D 0; b1. NÞ; Nþ/ D 0: (2.44)

Section 2.10 Bairstow’s Method 81

If Þ and þ denote initial approximations of NÞ and Nþ, we seek their improved values as Þ C dÞ and þ C dþ, which drive Eqs. (2.44) closer to zero. By expanding b0 and b1 in Taylor’s series about Þ and þ, we have

b0.Þ C dÞ; þ C dþ/ D b0.Þ; þ/ C 1Þ �

@b0 @Þ

� C 1þ

� @b0 @þ

� (2.45)

and

b1.Þ C dÞ; þ C dþ/ D b1.Þ; þ/ C 1Þ �

@b1 @Þ

� C 1þ

� @b1 @þ

� ; (2.46)

where the partial derivatives are evaluated at Þ and þ, and 1Þ and 1þ are given by

1Þ D NÞ � Þ and

1þ D Nþ � þ: (2.47) By setting Eqs. (2.45) and (2.46) to zero (assuming that Þ C dÞ ³ NÞ and þ C dþ ³ Nþ), we obtain

1Þ

� @b0 @Þ

� C 1þ

� @b0 @þ

� D �b0.Þ; þ/

and

1Þ

� @b1 @Þ

� C 1þ

� @b1 @þ

� D �b1.Þ; þ/: (2.48)

It can be seen that the evaluation of the partial derivatives of b0 and b1 in Eqs. (2.48) requires a recursive calculation of the partial derivatives of all the bi .i D n; n � 1; : : : ; 1; 0/ in Eq. (2.42):

@bn @Þ

D 0; @bn @þ

D 0;

@bn�1 @Þ

D �bn � Þ @bn @Þ

D �bn ; @bn�1

@þ D �Þ @bn

@þ D 0;

@bn�2 @Þ

D �bn�1 � Þ @bn�1 @Þ

� þ @bn @Þ

D �bn�1 C Þbn ; @bn�2

@þ D �Þ @bn�1

@þ � bn � þ @bn

@þ D �bn ;

:::

@b2 @Þ

D �b3 � Þ @b3 @Þ

� þ @b4 @Þ

;

@b2 @þ

D �Þ @b3 @þ

� b4 � þ @b4 @þ

;

82 Chapter 2 Solution of Nonlinear Equations

@b1 @Þ

D �b2 � Þ @b2 @Þ

� þ @b3 @Þ

;

@b1 @þ

D �Þ @b2 @þ

� b3 � þ @b3 @þ

;

@b0 @Þ

D �þ @b2 @Þ

;

@b0 @þ

D �b2 � þ @b2 @þ

: (2.49)

Here all the partial derivatives are evaluated at .Þ; þ/. Once the partial derivatives of b0 and b1 with respect to Þ and þ are evaluated, Eqs. (2.48) can be solved for 1Þ and 1þ as

1Þ D b1

� @b0 @þ

� � b0

� @b1 @þ

� �

@b0 @Þ

@b1 @þ

� @b1 @Þ

@b0 @þ

� (2.50)

and

1þ D b0

� @b1 @Þ

� � b1

� @b0 @Þ

� �

@b1 @þ

@b0 @Þ

� @b0 @þ

@b1 @Þ

� ; (2.51)

Once 1Þ and 1þ are known, new approximations to the correct values of NÞ and Nþ are computed as Þ C 1Þ and þ C 1þ, respectively.

The quadratic factor .x2CÞxCþ/ of the polynomial f (x) and the corresponding roots are computed iteratively as follows:

1. Choose initial approximations Þ and þ to NÞ and Nþ, respectively. 2. Compute the coefficients b1 and b0 using Eqs. (2.42).

3. Compute the partial derivatives of b1 and b0 with respect to Þ and þ using Eqs. (2.49).

4. Solve Eqs. (2.50) and (2.51) for 1Þ and 1þ.

5. Compute the new approximations for NÞ and Nþ as Þ C 1Þ and þ C 1þ. 6. By using the new values of Þ and þ computed in step 5, repeat steps 2 through

5 until convergence occurs (i.e., until the values of 1Þ and 1þ are found to be very small).

7. Using the converged values of Þ and þ for NÞ and Nþ, the roots of f (x) from the quadratic factor are computed as

x2 C NÞx C Nþ D 0; x D � NÞ š p

. NÞ/2 � 4 Nþ 2

: (2.52)

Section 2.10 Bairstow’s Method 83

Next the reduced polynomial of degree .n � 2/ is considered to obtain another quadratic factor. The procedure is continued until all the roots of f .x/ D 0 are found.

The Bairstow’s method can be used to determine the real, complex, or repeated roots of a polynomial equation. However, often it is difficult to select good initial approximations of Þ and þ. Another disadvantage with the method is that the roots found may not be very accurate.

�

Example 2.17

Find the roots of the polynomial equation

f .x/ D x4 � 10x3 C 35x2 � 50x C 24 D 0

using Bairstow’s method. Assume the initial approximations of NÞ and Nþ as Þ D �2:0 and þ D 1:0.

Solution

The coefficients of the polynomial equation are given by

a4 D 1:0; a3 D �10:0; a2 D 35:0; a1 D �50; and a0 D 24:0:

The coefficients b0; b1; Ð Ð Ð can be computed from Eqs. (2.42) as

b4 D 1:0; b3 D �8:0; b2 D 18:0; b1 D �6:0; and b0 D 6:0:

Equations (2.49) can be used to compute the partial derivatives of bi with respect to Þ and þ as

@bi @Þ

; i D 0; 1; 2; 3; 4 : �6:0; �5:0; 6:0; �1:0; 1:0

and @bi @þ

; i D 0; 1; 2; 3; 4 : �17:0; 6:0; �1:0; 0:0; 0:0:

The values of 1Þ and 1þ can be determined from Eqs. (2.50) and (2.51) as

1Þ D �0:54545456 and 1þ D 0:54545456:

The new approximations to NÞ and Nþ are given by

Þ D �2:54545456; and þ D 1:54545456:

Since 1Þ and 1þ are not small, we proceed to the next iteration with the new values of Þ D �2:54545456 and þ D 1:54545456. The values of Þ and þ found in different iterations are indicated in the following table:

84 Chapter 2 Solution of Nonlinear Equations

Iteration (i) Þ þ

— �2.0 1.0 1 �2.545455 1.545455 2 �2.848951 1.848951 3 �2.974668 1.974668 4 �2.999086 1.999086 5 �2.999984 1.000086 6 �3.000010 2.000008

From the converged values of Þ D �3:000010 and þ D 2:000008, the roots of the quadratic factor can be obtained from Eq. (2.52) as x1 D 1:000002 and x2 D 1:999982. The coefficients of the remaining polynomial are found to be

a2 D 1:0; a1 D �7:000016; and a0 D 12:000075:

Since the order of this remaining polynomial is 2, its roots can be determined readily as x3 D 3:000028 and x4 D 3:999900. The roots found can be compared with the exact values: x1 D 1; x2 D 2; x3 D 3, and x4 D 4. �

2.11 Muller’s Method

This method is based on quadratic interpolation and can find the roots of any equation f .x/ D 0. Let p1; p2, and p3 denote three approximations to the root with the corresponding values of f given by f1 D f .p1/, f2 D f .p2/ and f3 D f .p3/, respectively. Then a quadratic equation is assumed to pass through the three points. The proper root of the quadratic equation is then used to define another quadratic equation for the next approximation (or iteration). It can be seen that Muller’s method can be considered a generalization of the regula falsi method where a quadratic equation is used instead of a linear equation for interpolation. Since a higher degree equation is used for the interpolation, Muller’s method converges faster than the secant method. The method can be used to find both real and complex roots of a given equation.

To define the approximation quadratic, a new variable t is defined for conve- nience as

t D x � p3; (2.53) so that

t1 D p1 � p3; t2 D p2 � p3; (2.54)

Section 2.11 Muller’s Method 85

and t3 D p3 � p3 D 0:

The quadratic equation is assumed as

h.t/ D at2 C bt C c: (2.55) Since h(t) passes through the three points, .t1; f1/; .t2; f2/, and .t3; f3/, we have

f1 D at21 C bt1 C c; (2.56) f2 D at22 C bt2 C c; (2.57)

and f3 D at23 C bt3 C c: (2.58)

Equations (2.56) to (2.58) can be solved to obtain

c D f3; (2.59)

b D t 2 2 . f1 � f3/ � t21 . f2 � f3/

.t1t22 � t2t21 / ; (2.60)

and

a D t2. f1 � f3/ � t1. f2 � f3/ .t2t21 � t1t22 /

: (2.61)

By setting h.t/ D 0, the roots of Eq. (2.55) can be found as

t D �b š p

b2 � 4ac 2a

: (2.62)

The numerator and denominator of Eq. (2.62) are multiplied by b š p

b2 � 4ac, and Eq. (2.62) is rewritten as

t D �2c b š

p b2 � 4ac

: (2.63)

The roots are expressed in the form of Eq. (2.63) to avoid loss of significant digits in subtraction. The root having the smallest absolute value is selected in Eq. (2.63) for stability. The root can be expressed in terms of x as

t D t4 D x � p3 or

x D p4 D p3 C t4 D p3 � 2c b š

p b2 � 4ac

: (2.64)

For the next approximation, the two values out of p1; p2, and p3 that are closest to p4 are chosen as p1 and p2, and p4 is set equal to p3. In Fig. 2.25, for example, p2; p3, and p4 are chosen as p1; p2, and p3, respectively, for the next approximation (that is, p1, which lies farthest to p4 is discarded). Note that only one function evaluation is needed (at the new p3 or the previous p4) in each new iteration.

86 Chapter 2 Solution of Nonlinear Equations

x

f

f(x) h(t)

f1�f(t1) f2�f(t2)

f3�f(t3�0)

p1 p2 p3 t�t4

p4 t�t3�0t�t1 t�t2

Figure 2.25 Selection of points for quadratic interpolation.

The iterative process of Muller’s method can be described as follows:

1. Start with three approximations p1; p2, and p3 to the root of the equation f .x/ D 0 and a small number " to test the convergence of the method. Choose p1 D �1, p2 D 0, and p3 D 1, if no better approximations are available. Set i D 3 (i D 3 indicates that the next approximation to the root will be piC1 D p4).

2. Define ti�2 D pi�2 � pi ; ti�1 D pi�1 � pi , and ti D pi � pi D 0. 3. Find the function values

fi�2 D f .pi�2/; fi�1 D f .pi�1/; and fi D f .pi /: 4. Find the coefficients of the interpolating quadratic function (using relations

similar to Eqs. (2.59) to (2.61)) as

a D ti�1. fi�2 � fi / � ti�2. fi�1 � fi / ti�1t2i�2 � ti�2t2i�1

; (2.65)

b D t 2 i�1. fi�2 � fi / � t2i�2. fi�1 � fi /

ti�2t2i�1 � ti�1t2i�2 ; (2.66)

and c D fi : (2.67)

5. Compute the root of the interpolating quadratic equation as

t D �2c b š

p b2 � 4ac

(2.68)

by choosing the sign in the denominator that gives the smaller absolute value of t .

6. Find the new approximation to the root as

piC1 D pi � 2c

b š p

b2 � 4ac : (2.69)

Section 2.11 Muller’s Method 87

7. Determine fiC1 D f .piC1/ and test the convergence of the method: j fiC1j ": (2.70)

If convergence is achieved, stop the process by taking piC1 as the root of the equation f .x/ D 0. Otherwise, set the two values out of pi ; pi�1, and pi�2 that are closest to piC1 as pi and pi�1, update the iteration number as i D i C 1, and go to step 2.

�

Example 2.18

Find the root of the equation

f .x/ D x3 � 26x2 C 175x � 150 D 0 with the initial approximations p1 D �2; p2 D 0, and p3 D 3 and the convergence parameter " D 10�4 using Muller’s method.

Solution

The values of f (x) corresponding to p1; p2, and p3 are given by

f1 D f .�2:0/ D �352:0; f2 D f .0:0/ D �150:0; and f3 D f .3:0/ D 168:0: The values of t1; t2, and t3 can be determined as

t1 D p1 � p3 D �1:0 � 3:0 D �4:0; t2 D p2 � p3 D 0:0 � 3:0 D �3:0; and t3 D p3 � p3 D 3:0 � 3:0 D 0:0:

The coefficients of the approximating quadratic can be computed, using Eqs. (2.59) to (2.61), as

a D �24:0; b D 34:0; and c D 168:0: The value of t , given by Eq. (2.63), can be found as t D �2:03059673, which corresponds to p4 D p3 C t D 3:0 � 2:03059673 D 0:06940327 and f4 D f .p4/ D �3:87673950. Since j f4j > ", we proceed to the next iteration with the new p1 D 0:0; p2 D 3:0, and p3 D 0:96940327.

The corresponding values of the function f (x) are given by

f1 D f .0:0/ D �150:0; f2 D f .3:0/ D 168:0; and f3 D f .0:96940327/ D �3:87673950:

The values of t1; t2, and t3 can be determined as

t1 D p1 � p3 D 0:0 � 0:96940327 D �0:96940327; t2 D p2 � p3 D 3:0 � 0:96940327 D 2:03059673;

and t3 D p3 � p3 D 0:0:

The coefficients of the approximating quadratic can be found from Eqs. (2.59) to (2.61) as a D �22:03062439; b D 129:37876892, and c D �3:87673950. The value

88 Chapter 2 Solution of Nonlinear Equations

of t can be found, using Eq. (2.63), as t D 0:03011874, which yields p4 D p3 C t D 0:99952203 and f4 D f .p4/ D �0:06024170. Since j f4j > ", we continue with the iterative process.

For the next iteration, we take p1 D 3:0; p2 D 0:96940327, and p3 D 0:99952203. These correspond to f1 D f .p1/ D 168:0; f2 D f .p2/ D �3:87673950; f3 D f .p3/ D �0:06024170; t1 D p1 � p3 D 2:00047779; t2 D p2 � p3 D �0:03011876, and t3 D p3 � p3 D 0:0. The coefficients of the approximating quadrat- ic are given by Eqs. (2.59) to (2.61) as a D �21:03071594; b D 126:08151245, and c D �0:06024170. These result in a value of t D 0:00047784, which corresponds to p4 D p3 C t D 0:99999988 and f4 D f .p4/ D �0:00001526. Since j f4j < " D 10�5, we stop the iterative process and take the root as x D 0:99999988. This value can be compared with the exact value of 1; the other two roots being x D 10 and x D 15. �

2.12 Newton–Raphson Method for Simultaneous Nonlinear Equations

Many engineering analyses require the solution of a set of simultaneous nonlinear algebraic equations. The Newton–Raphson method can be extended to solve a system of nonlinear equations. Consider the following system of nonlinear equations:

f1.x1; x2; : : : ; xn/ D 0; f2.x1; x2; : : : ; xn/ D 0;

:::

fn.x1; x2; : : : ; xn/ D 0: (2.71)

If the functions f1; f2; : : : ; fn are expanded in Taylor’s series about an arbi- trary (guess or approximate) solution .x1; x2; : : : ; xn/ and if only the linear terms are retained, we can obtain the values of the functions at the exact solution . Nx1; Nx2; : : : ; Nxn/ as

f1. Nx1; Nx2; : : : ; Nxn/ ³ f1.x1; x2; : : : ; xn/ C �

@ f1 @x1

� . Nx1 � x1/ C

� @ f1 @x2

� . Nx2 � x2/ C

Ð Ð Ð C �

@ f1 @xn

� . Nxn � xn/; (2.72)

:::

fn. Nx1; Nx2; : : : ; Nxn/ ³ fn.x1; x2; : : : ; xn/ C �

@ fn @x1

� . Nx1 � x1/ C

� @ fn @x2

� . Nx2 � x2/ C

Ð Ð Ð C �

@ fn @xn

� . Nxn � xn/; (2.73)

Section 2.12 Newton–Raphson Method for Simultaneous Nonlinear Equations 89

where the partial derivatives are evaluated at .x1; x2; : : : ; xn/. Equations (2.72) to (2.73) provide corrections to the assumed solution .x1; x2; : : : ; xn/. By setting Eqs. (2.72) to (2.73) to zero and denoting 1xi D Nxi � xi .i D 1; 2; : : : ; n/ as the corrections to the variables xi .i D 1; 2; : : : ; n/, we can rewrite Eqs. (2.72) to (2.73) as2

66666666666666666664

@ f1 @x1

@ f1 @x2

@ f1 @xn

@ f2 @x1

@ f2 @x2

@ f2 @xn

Ð Ð Ð Ð Ð Ð

Ð Ð Ð Ð Ð Ð

@ fn @x1

@ fn @x2

@ fn @xn

3 77777777777777777775

8>>>>>>>>>>>>< >>>>>>>>>>>>:

1x1

1x2

Ð Ð Ð

1xn

9>>>>>>>>>>>>= >>>>>>>>>>>>;

D

8>>>>>>>>>>>>< >>>>>>>>>>>>:

� f1 � f2

Ð Ð Ð

� fn

9>>>>>>>>>>>>= >>>>>>>>>>>>;

; (2.74)

where the partial derivatives @ fi @xj

and the functions fi are evaluated at the current approximation .x1; x2; : : : ; xn/. The matrix of first partial derivatives of fi , on the left- hand side of Eq. (2.74), is known as the Jacobian of fi . Equations (2.74) denote a sys- tem of linear equations in the unknowns 1x1; 1x2; : : : ; 1xn . By solving Eqs. (2.74), we find 1xi .i D 1; 2; : : : ; n/ and determine the next approximate solution as

x .iC1/1 D x .i/1 C 1x .i/1 ; x .iC1/2 D x .i/2 C 1x .i/2 ;

::: (2.75)

x .iC1/n D x .i/n C 1x .i/n ; where the superscripts i and i C1 denote the current and next iteration, respectively. The iterative process can be stated as follows:

1. Start with an initial guess solution .x .1/1 ; x .1/ 2 ; : : : ; x

.1/ n / and a small number " to

test the convergence of the process. Set the iteration number, i D 1. 2. Evaluate the function values:

f j .x .i/ 1 ; x

.i/ 2 ; : : : ; x

.i/ n /; j D 1; 2; : : : ; n: (2.76)

3. Find the partial derivatives of the functions f j as

@ f j .x .i/ 1 ; x

.i/ 2 ; : : : ; x

.i/ n /

@xk ; j D 1; 2; : : : ; n; k D 1; 2; : : : ; n: (2.77)

4. Solve the Eqs. (2.74) and find 1x .i/j ; j D 1; 2; : : : ; n.

90 Chapter 2 Solution of Nonlinear Equations

5. Find the new solution as

x .iC1/j D x .i/j C 1x .i/j ; j D 1; 2; : : : ; n: (2.78)

6. Evaluate the function values

f j .x .iC1/ 1 ; x

.iC1/ 2 ; : : : ; x

.iC1/ n /; j D 1; 2; : : : ; n; (2.79)

and test for the convergence of the process:

j f j .x .iC1/1 ; x .iC1/2 ; : : : ; x .iC1/n /j "; j D 1; 2; : : : ; n: (2.80)

If the convergence criteria of Eq. (2.80) are satisfied, stop the process by taking the solution as .x .iC1/1 ; x

.iC1/ 2 ; : : : ; x

.iC1/ n /. Otherwise, update the iteration

number as i D i C 1 and go to step 3.

�

Example 2.19

Find the solution of the equations

f .x1; x2/ D x21 C x22 � 8x1 � 4x2 C 11 D 0 (E1) and

f2.x1; x2/ D x21 C x22 � 20x1 C 75 D 0 (E2) by taking the starting point as .x1 D 2; x2 D 4/ and the value of " as 10�5.

Solution

The partial derivatives of f1 and f2 can be found as @ f1 @x1

D 2x1 � 8; @ f1 @x2

D 2x2 � 4 and

@ f2 @x1

D 2x1 � 20; @ f2 @x2

D 2x2:

At the starting point .x .1/1 D 2; x .1/2 D 4/, the function values are given by f1.2; 4/ D 22 C 42 � 8.2/ � 4.4/ C 11 D �1

and f2.2; 4/ D 22 C 42 � 20.2/ C 75 D 55:

Since j f1j > " and j f2j > ", we compute the Jacobian as

[J ]j .x .1/1 ;x

.1/ 2 /

D 2 4 2x1 � 8 2x2 � 4

2x1 � 20 2x2

3 5 þþþþþþ .2;4/

D 2 4 �4 4

�16 8

3 5

Section 2.12 Newton–Raphson Method for Simultaneous Nonlinear Equations 91

and solve the equations2 4 �4 4

�16 8

3 5 8< : 1x

.1/ 1

1x .1/2

9= ; D

8< : � f1� f2

9= ; D

8< : 1�55

9= ;

to obtain 1x .1/1 D 7:125 and 1x .1/2 D 7:375:

This yields the new point as

x .2/1 D x .1/1 C 1x .1/1 D 2 C 7:125 D 9:125 and

x .2/2 D x .1/2 C 1x .1/2 D 4 C 7:375 D 11:375:

The values of the functions f1 and f2 at .x .2/ 1 D 9:125; x .2/2 D 11:375/ are given by

f1 D 105:15625 and f2 D 105:15625. Since j f1j > " and j f2j > ", the Jacobian is computed as

[J ]j .x .2/1 ;x

.2/ 2 /

D 2 4 10:25 18:75

�1:75 22:75

3 5 :

The solution of the equations2 4 10:25 18:75

�1:75 22:75

3 5 8< : 1x

.2/ 1

1x .2/2

9= ; D

8< : � f1� f2

9= ; D

8< : �105:15625�105:15625

9= ;

is given by 1x .2/1 D �1:581295; 1x .2/2 D �4:743892:

The new approximate solution can be found as

x .3/1 D x .2/1 C 1x .2/1 D 9:125 � 1:581295 D 7:543705 and

x .3/2 D x .2/2 C 1x .2/2 D 11:375 � 4:743892 D 6:631108: Since

j f1.x .3/1 ; x .3/2 /j D j5:141006j D 5:141006 > " and

j f2.x .3/1 ; x .3/2 /j D j5:141037j D 5:141037 > "; we proceed to the next iteration. The results of the iterative process are given in the following table (note that the iterative process converged in eight iterations): �

92 Chapter 2 Solution of Nonlinear Equations

Iteration x.i/1 x .i/ 2 f1.x

.i/ 1 ; x

.i/ 2 / f2.x

.i/ 1 ; x

.i/ 2 / [J ]j.x.i/1 ;x

.i/ 2 /

1x1 1x2

(i)

1 2.0 4.0 �1.0 55.0 2 4 �4 4

�16 8

3 5 7.125 7.375

2 9.125 11.375 105.15625 105.15625

2 4 10:25 18:75

�1:75 22:75

3 5 �1.581295 �4.743892

3 7.543705 6.631108 25.00502 25.00499

2 4 7:087410 9:262217

�4:912590 13:262217

3 5 �0.717011 �2.151025

4 6.826694 4.480083 5.141006 5.141037

2 4 5:653387 4:960167

�6:346613 8:960167

3 5 �0.250366 �0.751102

5 6.576328 3.728981 0.626846 0.626831

2 4 5:152657 3:457962

�6:847343 7:457962

3 5 �0.040374 �0.121117

6 6.535954 3.607864 0.016280 0.016296

2 4 5:071909 3:215729

�6:928091 7:215729

3 5 �0.001105 �0.003319

7 6.534849 3.604545 0.000019 0.000015

2 4 5:069698 3:209089

�6:930302 7:209089

3 5 �0.000002 �0.000004

8 6.543847 3.604541 0.000004 0.000000 — — —

2.13 Unconstrained Minimization

It can be seen that the root of the equation f .x/ D 0 is same as the minimum of the function

F.x/ D f 2.x/: (2.81) Similarly the solution of the simultaneous linear or nonlinear equations

f1.x1; x2; : : : ; xn/ D 0; f2.x1; x2; : : : ; xn/ D 0;

:::;

fn.x1; x2; : : : ; xn/ D 0 (2.82) is same as the unconstrained minimum of the function

F.x1; x2; : : : ; xn/ D f 21 .x1; x2; : : : ; xn/ C f 22 .x1; x2; : : : ; xn/ C Ð Ð Ð C f 2n .x1; x2; : : : ; xn/: (2.83)

The necessary and sufficient conditions for the minimum of F(x), given by Eq. (2.81), are

d F.xŁ/ dx

D 0 (2.84)

Section 2.14 Convergence of Methods 93

and d2 F.xŁ/

dx2 > 0; (2.85)

where xŁ is the minimum of F(x) or the root of the equation f .x/ D 0. Similarly, the necessary and sufficient conditions for the minimum of F.x1; x2; : : : ; xn/, given by Eq. (2.82), are

@ F.xŁ1 ; x Ł 2 ; : : : ; x

Ł n /

@xi D 0; i D 1; 2; ::; n (2.86)

and

[H.xŁ1 ; x Ł 2 ; : : : ; x

Ł n /] D

" @2 F.xŁ1 ; x

Ł 2 ; : : : ; x

Ł n /

@xi @xj

# D positive definite: (2.87)

Although Eqs. (2.84) and (2.85) or (2.86) and (2.87) can be solved to find the minimum of F , it is more convenient to find the minimum of F(x) or F.x1; x2; : : : ; xn/ through numerical methods. The numerical methods of finding the solution of simultaneous equations through the minimization of the function F(x) or F.X1; x2; : : : ; xn/ are presented in Chapter 12.

2.14 Convergence of Methods

2.14.1 Bisection Method

The interval at the end of i th bisection, hiC1, can be expressed as

hiC1 D b.iC1/ � a.iC1/ D 12 hi D 1 22

hi�1 D 12iC1 h0 D 1

2iC1 .b.1/ � a.1//: (2.88)

Since the root of f .x/ D 0 can lie anywhere in the interval hiC1, the error, after ith bisection ei , can be expressed as

jei j hiC1 D 1

2iC1 h0: (2.89)

The bisection method is independent of the nature of the function f (x) and always converges at the same rate. The inequality given by Eq. (2.89) is called an a priori error bound. By taking logarithms on both sides of the inequality (2.89), we obtain

log10 jei j log10 h0 � .i C 1/ log10 2; (2.90) which can be rearranged to obtain

i ²

log10 h0 � log10 jei j log10 2

¦ � 1: (2.91)

Thus, the inequality (2.90) gives an a priori bound on the number of bisections or iterations required to achieve a specified accuracy, jei j.

94 Chapter 2 Solution of Nonlinear Equations

�

Example 2.20

Determine the number of iterations (bisections) required in a problem with h0 D 0:05, if an accuracy of 0.0001 is desired.

Solution

In this example, h0 D 0:5. Since the accuracy desired is 0.0001, we have " D 0:00005. Hence, Eq. (2.90) gives

i ²

log10 0:5 � log10 0:00005 log10 2

¦ � 1

²

log10 10;000 log10 2

¦ � 1

12:28771238: By rounding the number of bisections to an integral number, we find that at most 13 bisections are required to achieve the specified accuracy. �

2.14.2 Newton’s Method

The Taylor’s series expansion can be used to find the value of f (x) at xiC1 as

f .xiC1/ D f .xi / C .xiC1 � xi / f 0.xi / C 12 .xiC1 � xi / 2 f 00.xi /; xi ¾ xiC1: (2.92)

If xiC1 D xŁ is the root of f .x/ D 0, Eq. (2.92) gives

f .xŁ/ D 0 D f .xi / C .xŁ � xi / f 0.xi / C 12 .x Ł � xi /2 f 00.¾/: (2.93)

Newton’s iterative formula, Eq. (2.14), can be rewritten as

f .xi / C .xiC1 � xi / f 0.xi / D 0: (2.94) Subtraction of Eq. (2.94) from (2.93) yields

.xŁ � xiC1/ f 0.xi / C 12 .x Ł � xi /2 f 00.¾/ D 0: (2.95)

Denoting the error in iteration i; ei , as

ei D xŁ � xi : (2.96) Eq. (2.95) can be expressed as

eiC1 f 0.xi / C 12e 2 i f

00.¾/ D 0: (2.97) As the method converges, both xi and ¾ can be replaced by the actual root xŁ so that Eq. (2.96) can be rewritten as

eiC1 D � f 00.xŁ/

2 f 0.xŁ/ e2i : (2.98)

Section 2.14 Convergence of Methods 95

Equation (2.98) indicates that the error in an iteration is proportional to the square of the error in the previous iteration. This means that the number of correct decimal places in the numerical solution nearly doubles with each iteration. This type of convergence is termed quadratic convergence, and hence Newton’s method is considered to be a quadratically convergent method.

�

Example 2.21

Show that the Newton–Raphson method is quadratically convergent using the results of Example 2.8.

Solution

The results of Example 2.8 can be summarized as follows:

i xi Error, Nei D xŁ � xi Percent error 1 0.0 0.480864 100.0

2 0.364649 0.116215 24.17

3 0.461671 0.019193 3.99

4 0.480193 0.000671 0.14

5 0.480864 ³ 0:0 ³ 0:00

Equation (2.98) gives the relation between errors in two consecutive itera- tions as

eiC1 D � f 00.xŁ/

2 f 0.xŁ/ e2i : (E1)

The first and second derivatives of f (x) considered in Example 2.8 are given by

f 0.x/ D 2:8 � 3:2x 2

.1 C x2/3 (E2)

and

f 00.x/ D �23:2x C 12:8x 2

.1 C x2/4 : (E3)

Since x1 D 0:0; f 0.x1/ D 2:8, and f 00.x1/ D 0:0, Eq. (E1) gives e2 D .0/ e21 D 0. Hence we start the computations with x2. At x2 D 0:364649; f 0.x2/ D 1:6327441 and f 00.x2/ D �4:7577443, so that e3 D 1:4569779 e22. Using e2 D 0:116215, we obtain e3 D 0:0196778, which can be seen to be close to the true error, Ne3 D 0:019193. For the next iteration, x3 D 0:461671; f 0.x3/ D 1:186268917; f 00.x3/ D �4:363599799, and Eq. (E1) gives, with e3 D 0:019193; e4 D 1:839211892 e23 D 0:000677513, which compares very well with the true error of Ne4 D 0:000671. For the next iteration

96 Chapter 2 Solution of Nonlinear Equations

x4 D 0:480193, f 0.x4/ D 1:106573381; f 00.x4/ D �4:239961723, and Eq. (E1) gives, with e4 D 0:000671, e5 D 1:915806848 e24 D 0:000000863, which can be considered to be approximately equal to zero. Thus, the results prove the quadratic convergence of the Newton–Raphson method. �

2.14.3 Fixed Point Iteration Method

The iterative procedure gives xiC1 D g.xi /: (2.99)

If xŁ denotes the true solution, then

xŁ D g.xŁ/: (2.100) Subtraction of Eq. (2.99) from (2.100) gives

xŁ � xiC1 D g.xŁ/ � g.xi /: (2.101) Dividing both sides of Eq. (2.101) by xŁ � xi yields

xŁ � xiC1 xŁ � xi D

g.xŁ/ � g.xi / xŁ � xi : (2.102)

From the mean value theorem, the right-hand side of Eq. (2.102) can be rewritten as

g.xŁ/ � g.xi / xŁ � xi D g

0.¾/; xŁ ¾ xi ; (2.103)

where g0.¾/ denotes the derivative of g(x) at x D ¾ , which lies somewhere between x and xi . Equations (2.102) and (2.103) yield

eiC1 D g0.¾/ei ; (2.104) where ei and eiC1 denote the errors in iterations i and i C 1, respectively. Note that the method converges only if g0.¾/ < 1, since in this case, the error in iteration i C 1 will be smaller than the error in iteration i . If g0.¾/ < 1, the error decreases with successive iteration. Hence the method is said to exhibit linear convergence.

2.14.4 Secant Method

The error in the secant method can be shown to follow the relation [2.17]

jeiC1j D cjei jp; (2.105) where

c D þþþþ f 00.xŁ/2 f 0.xŁ/

þþþþ (2.106) and p D 1:618. Thus, the secant method is considered to have order of convergence p D 1:618.

Section 2.15 Choice of the Method 97

2.15 Choice of the Method

The graphical method can be used to find the approximate, real roots of nonlinear equations. As such, the method can be used to find good starting points for use with other more accurate and efficient methods. The incremental search method can be used to find the real roots of any function. This method cannot find complex roots. The method fails to find the root when the function f (x) does not undergo a sign change at the root. Although the method is not efficient, it can be used to find an approximate value of the root for use as a starting point for other more efficient methods. The bisection method, although slow, always finds the root. If f (x) can be evaluated quickly, the bisection method is recommended. The accuracy improves with each iteration in this method. The bisection method cannot be used to find the complex roots of an equation.

If the function f (x) is continuous with continuous derivatives, the New- ton–Raphson method can be used to find the root efficiently. However, the method may not converge if the value of the derivative of f is close to zero. The method may also not converge if the initial guess is not good. The secant method is simi- lar to the Newton–Raphson method where the derivative is approximated by the chord. This method is preferred over the Newton–Raphson method, when the evaluation of f 0.x/ is difficult. The regula falsi method uses linear interpolation of the function in finding the root. The Newton–Raphson, secant, and regula falsi methods require smaller number of iterations to find a root. Although the secant method gives better results than the regula falsi method, it may not converge in some cases. The Newton–Raphson, secant, and regula falsi methods can find com- plex roots provided the initial approximations are given as complex numbers. Note that a complex root cannot be found by giving a real number as the initial guess. The secant and regula falsi methods require two initial guess values of the root, while the Newton–Raphson method requires only one initial guess value of the root.

The fixed point iteration is simple and efficient, if the equation f .x/ D 0 can be rewritten as x D g.x/. This method may not converge if jg0.x/j > 1. Although Muller’s method requires three initial approximations for the root, it is particularly efficient in finding the complex roots of equations. Since the method is based on a quadratic approximation of f (x), it converges faster than the secant method.

For a polynomial, if all the roots are required, Bairstow’s method is recom- mended. After finding a quadratic factor, Bairstow’s method is to be applied to the deflated polynomial. In this method, if the locations of some of the roots are known, these roots must be found first.

The Newton–Raphson method is very efficient in solving systems of nonlinear equations. The method may not converge if the initial approximation is not good. Since the method requires the first derivatives of the functions, we need to use either the exact expressions or the finite difference equivalents for the derivatives of the functions involved.

98 Chapter 2 Solution of Nonlinear Equations

2.16 Use of Software Packages

All the software packages have built-in programs (functions) for finding roots of nonlinear equations, roots of polynomial equations, and solutions of simultaneous nonlinear equations. Usually, all the roots of a polynomial equation are determined by the program without requiring initial guess values for the roots. However, most programs used for finding the roots of nonlinear equations require input in the form of lower and upper bounds (range) or an initial guess for the root. In addition, some symbolic programs can find the solution in the form of a mathematical expression.

2.16.1 MATLAB

To find the roots of a polynomial equation use the function roots(p) where p is a row vector of coefficients of the polynomial in descending order of the power of the variable. The roots of the polynomial are returned in the form of a column vector. To find the roots of a nonlinear equation, the function fzero(y, x1) is to be used. Here y defines the nonlinear function and x1 indicates the initial estimate (starting value) of the root. Several illustrative examples follow:

�

Example 2.22

The roots of f .x/ D x12 � 2 D 0 are as follows >> roots ([1 0 0 0 0 0 0 0 0 0 0 0 -2]) ans =

-1.0595 -0.9175 + 0.5297i -0.9175 - 0.5297i -0.5297 + 0.9175i -0.5297 - 0.9175i 0.0000 + 1.0595i 0.0000 - 1.0595i 0.5297 + 0.9175i 0.5297 - 0.9175i 1.0595 0.9175 + 0.5297i 0.9175 - 0.5297i

�

Example 2.23

The roots of f .x/ D tan�1 x D 0. Note: To find the roots of a general nonlinear equation, the lower bound, increment for search and upper bound on x are to be given before using the function fzero. In addition, the value of x in the vicinity of the desired root is to be given as an argument of fzero.

Section 2.16 Use of Software Packages 99

>> x=-10:0.01:10; >> f=’atan(x)’; >> root=fzero(f,1.6) Zero found in the interval: [-0.448, 3.0482].

root =

1.5040e-20

>> root=fzero(f,0.2) Zero found in the interval: [-0.056, 0.38102].

root =

-5.9260e-24

>>

�

Example 2.24

The roots of f .x/ D sin 2x D 0 are

>> x=-5:0.01:5; >> y4=’sin(2*x)’; >> x4=fzero(y4,0.75) x4 =

9.2059e-023 >> x4=fzero(y4,2.75) x4 =

3.1416 >> x4=fzero(y4,-2.75) x4 =

-3.1416

�

Example 2.25

The roots of

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 1 x

C 0:65x .1 C x2/ D 0 :

>> x=-2:0.01:2; >> y2=’1.5*x./(1+x.^2).^2-0.65*(atan(1./x))+0.65*x./(1+x.^2)’; >> fplot(y2,[-2 2]) >> xlabel(‘x’),ylabel(‘y=f(x)’), >> title(‘y=1.5*x/(1+x^2)^2-0.65* (atan(1/x))+0.65*x/(1+x^2)’) >> x2=fzero(y2,-0.1)

100 Chapter 2 Solution of Nonlinear Equations

x2 =

2.4831e-016

>> x2=fzero(y2,0.1)

x2 =

-2.4831e-016

>> x2=fzero(y2,0.8)

x2 =

0.4809

>> x2=fzero(y2,-0.8)

x2 =

-0.4809

�

Example 2.26

Solve the following simultaneous nonlinear equations:

f .x; y/ D x2 C y2 � 8x � 4y C 11 D 0

and g.x; y/ D x2 C y2 � 20x C 75 D 0:

(a) By minimizing F D f 2.x; y/ C g2.x; y/ with initial estimates of .x D 2; y D 4/, and .x D 1; y D 1/. Note: The function F is defined as f 3 in an m-file.

Listing of f3.m:

function f=f3(x)

f=(x(1)^2+x(2)^2-8*x(1)-4*x(2)+11)^2+(x(1)^2+x(2)^2-20*x(1) +75)^2;

>> [x3]=fmins(‘f3’,[2,4]); >> x3

Section 2.16 Use of Software Packages 101

x3 =

5.0652 -0.8045

>> [x3]=fmins(‘f3’,[1,1]); >> x3

x3 =

5.0652 -0.8045

>>

(b) By using the function solve. This gives all the solutions of the equations.

Note: A blank space is to be left on both sides of the equality sign in defining the equations.

>> [x,y] = solve(‘x^2+y^2-8*x-4*y+11=0’,‘x^2+y^2-20*x+75 = 0’)

x =

[ 29/5+3/10*6^(1/2)] [ 29/5-3/10*6^(1/2)]

y =

[ 7/5+9/10*6^(1/2)] [ 7/5-9/10*6^(1/2)]

>>

2.16.2 MAPLE

The solve command finds the all the possible roots of a nonlinear equation, including polynomial equations. This command provides symbolic solutions of equations, if they exist. The command fsolve uses Newton–Raphson method to find the roots of a general nonlinear equation. For a nonlinear equation, it gives one root, while the command provides all the real roots for a polynomial equation. In all the roots of a polynomial equation, including real and complex ones, the word complex is to be included in the fsolve command. The following examples illustrate the procedures.:

�

Example 2.27

Roots of f .x/ D x12 � 2 D 0:

102 Chapter 2 Solution of Nonlinear Equations

" > eqn1:=x 1̂2-2;

eqn1 :D x12 � 22 66666666666664

> fsolve(eqn1 =0, x, complex);

fx D �1:059463094g; fx D �:9175219541 � :5297315472 I g; fx D �:9175219541 C :5297315472 I g; fx D �:5297315472 � :9175219541 I g; fx D :5297315472 C :9175219541 I g; fx D �1:059463094 I g; fx D 1:059463094 I g; fx D :5297315472 � :9175219541 I g; fx D :5297315472 C :9175219541 I g; fx D :9175219541 � :5297315472 I g; fx D :9175219541 C :5297315472 I g; fx D 1:059463094g2

66666666664

> evalf(solve(eqn1=0, x));

fx D 1:059463094g; fx D �1:059463094g; fx D 1:059463094 I g; fx D �1:059463094 I g; fx D :9175219535 C :5297315470 I g; fx D �:9175219535 � :5297315470 I g; fx D :5297315475 � :9175219545 I g; fx D �:5297315475 C :9175219545 I g; fx D :9175219535 � :5297315470 I g; fx D �:9175219535 C :5297315470 I g; fx D :5297315475 C :9175219545 I g; fx D �:5297315475 � :9175219545 I g

�

Example 2.28

Roots of f .x/ D tan�1x D 0:" > eqn:=arctan(x);

eqn :D arctan.x/" > fsolve(eqn=0, x);

fx D 0g

�

Example 2.29

Roots of

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0 :

(a) By giving different ranges for the root:

> eq:=1.5*x/(1.0+x^2.0)^2.0-0.65*arctan(1.0/x)+0.65*x/(1+x^2) =0:

> fsolve(eq, x, -1.0..1.0); -.4808644853

Section 2.16 Use of Software Packages 103

If the range does not contain a root, the program returns just the equation and the range given:

> eq:=1.5*x/(1.0+x^2.0)^2.0-0.65*arctan(1.0/x)+0.65*x/(1+x^2) =0:

> fsolve(eq, x, -0.48..0.48);

x x fsolve(1.5 ------------ - .65 arctan(1.01/x)+ .65 ------ = 0.

2.0 2.0 2 (1.0 + x ) 1 + x

x, -.48 .. .48)

h >2 664 > f:=1.5*x/(1+x 2̂) 2̂-0.65*(arctan(1/x))+0.65*x/(1+x 2̂):

> fsolve(f,x,0..1);

:4808644853h >

(b) Without giving a range for the root:2 4> eqn2:=1.5*x/(1+x 2̂) 2̂-0.65*(arctan(1/x))+0.65*x/(1+x 2̂);

eqn2 :D 1:5 x .1 C x2/2 � :65arctan

� 1 Px �

C :65 x 1 C x2"

> fsolve(eqn2=0,x);

fx D �:4808644853g

�

Example 2.30

Roots of f .x/ D x4 � 12x3 C 49x2 � 78x C 40 D 0: > eqn:=x^4.0-12.0*x^3.0+49.0*x^2.0-78.0*x+40.0=0;

4.0 3.0 2.0 eqn := x - 12.0 x + 49.0 x - 78.0 x + 40.0 = 0

> solve(eqn, {x}); {x = 5.}, {x = 1.}, {x = 2.}, {x = 4.}

> poly :=x^4-12.*x^3+49.*x^2-78.*x+40.:fsolve(poly,x,maxsols=4 ); 1., 2.000000000, 4.000000000, 5.000000000

solve(x^4-12.*x^3+49.*x^2-78.*x+40.,x); 5., 1., 2., 4.

104 Chapter 2 Solution of Nonlinear Equations

roots(x^4-12*x^3+49*x^2-78*x+40); [[1, 1], [2, 1], [4, 1], [5, 1]]

�

Example 2.31

Solution of simultaneous nonlinear equations:

f .x; y/ D x2 C y2 � 8x � 4y C 11 D 0 and

g.x; y/ D x2 C y2 � 20x C 75 D 0

(a) By giving ranges for x and y:

A different range gives a different answer:h > f:=x 2̂+y 2̂-8*x-4*y+11=0:h > g:=x 2̂+y 2̂-20*x+75=0:" > fsolve(f,g,x,y,x=0..100,y=-29.8..100);

fx D 6:534846923; y D 3:604540769g" > fsolve(f,g,x,y,x=0..100,y=-29.9..100);

fx D 5:065153077; y D �:8045407685g" > fsolve(f,g,x,y);

fx D 5:065153077; y D �:8045407685gh >

(b) Without giving ranges for variables:

" > fsolve(x 2̂+y 2̂-8*x-4*y+11=0,x 2̂+y 2̂-20*x+75=0);

fx D 5:065153077; y D �:8045407685gh >

2.16.3 MATHCAD

To find the roots of polynomial equations, one can use the polyroots (v) command. Here v denotes the vector of coefficients of the polynomial starting with the constant term. All the roots of the equation are returned as a vector. To find the root of a nonlinear equation, the root(f(x),x) command can be used. Here f (x) is the nonlinear function and x is the variable. By giving a guess value of the root, the program returns the root closest to the guess value. To find the solution of several nonlinear equations, we need to use the appropriate solve blocks. The procedures are illustrated with the following examples.

Section 2.16 Use of Software Packages 105

�

Example 2.32

Roots of f .x/ D x12 � 2 D 0 (V1 denotes the vector of polynomial coefficients): x12 � 2

V1 :D

0 BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB@

�2 0

0

0

0

0

0

0

0

0

0

0

1

1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCA

–1.059

0

–0.918 + 0.53i

polyroots(V1)=

–0.918 – 0.53i

–0.53– 0.918i

–0.53 +0.918i

–3.841.10–9 + 1.059i

–1.059i

0.53– 0.918i

0.53 +0.918i

0.918 + 0.53i

0

1

2

3

4

5

6

7

8

9

0.918 – 0.53i

1.059

10

11

�

Example 2.33

The roots of f .x/ D tan�1 x D 0: with different initial guess values:

106 Chapter 2 Solution of Nonlinear Equations

x :D 1:6 root(atan(x),x) D �5:335 ð 10�5

x :D 0:2 root(atan(x),x) D 6:804 ð 10�5

�

Example 2.34

Roots of f .x/ D sin 2x D 0: with different initial guess values: x :D 0:75 root(sin(2Ðx),x) D �6:283 x :D 0:0 root(sin(2Ðx),x) D 0

�

Example 2.35

Roots of f .x/ D 1:5x .1Cx2/2 �0:65 tan�1 x C 0:65x1Cx2 D 0 using different initial guess values:

x :D �1:2

root

"" .1:5 Ð x/ð

.1 C x2/2Ł #

� 0:65 Ð atan �

1 x

� C

.0:65 Ð x/ .1 C x2/

½ ; x

# D �11:222

x :D 0:8

root

"" .1:5 Ð x/ð

.1 C x2/2Ł #

� 0:65 Ð atan �

1 x

� C

.0:65 Ð x/ .1 C x2/

½ ; x

# D 12:001

x :D �0:8

root

"" .1:5 Ð x/ð

.1 C x2/2Ł #

� 0:65 Ð atan �

1 x

� C

.0:65 Ð x/ .1 C x2/

½ ; x

# D �12:001

�

Example 2.36

Solution of simultaneous nonlinear equations:

f .x; y/ D x2 C y2 � 8x � 4y C 11 D 0 and

g.x; y/ D x2 C y2 � 20x C 75 D 0:

Section 2.17 Computer Programs 107

x :D 2 y :D 4 Given

x2 C y2 � 8 Ð x � 4 Ð y C 11 D 0 x2 C y2 � 20 Ð x C 75 D 00 @xval

yval

1 A :D Find(x,y)

xval D 6:535 yval D 3:605

2.17 Computer Programs

2.17.1 Fortran Programs

Program 2.1 Secant method

The input data required are as follows:

X1,X2 D initial (guess) values of the root MITER D maximum number of iterations permitted (usual value: 100) EPS D convergence criterion (usual value: 0.00001) Subprogram to evaluate the function f .x/ in the following form:

FUNCTION F(X) F=... RETURN END

Illustration (Example 2.5)

The listing and output of Program 2.1 are as follows:

PROGRAM FOR SECANT METHOD DATA X1,X2,MITER,EPS/0.1,0.5,100,0.00001/ I=0 PRINT 10

10 FORMAT (‘SOLUTION BY SECANT METHOD’,/) Y=X1 YY=X2 Z=F(Y) ZZ=F(YY) DO WHILE(1)

I=I+1

108 Chapter 2 Solution of Nonlinear Equations

IF (I .GT. MITER) Then PRINT 90

90 FORMAT (/, ‘NO CONVERGENCE IN MITER ITERATIONS’) STOP END IF X3=YY-ZZ*(YY-Y)/(ZZ-Z) F3=F(X3) PRINT 30, I, Y, Z, YY, ZZ

30 FORMAT(‘ITER=’,I3,’x1=’,E12.6,’F(x1)=’,E12.6,/,10x, & & ‘ x2 = ’, E12.6,’ F (x2) = ’, E12.6)

PRINT 40, X3, F3 40 FORMAT (9X, ‘NEWX = ’, E12.6, 1X, ‘NEWF = ’, E12.6,/)

IF (ABS (F3) .LE. EPS) THEN PRINT 60, X3, F3

60 FORMAT (‘CONVERGED SOLUTION: X = ’, E15.8,’ F = ’, E15.8) STOP ELSE

Y=YY YY=X3 Z=ZZ ZZ=F3 END IF

END DO END

FUNCTION F(X) F=1.5*X/((1.0+X*X)**2)-0.65*ATAN(1.0/X)+0.65*X/(X*X+1.0) RETURN END

SOLUTION BY SECANT METHOD

ITER = 1 X1 = 0.100000E+00 F(X1) = -.744832E+00 X2 = 0.500000E+00 F(X2) = 0.203533E-01

NEWX = 0.489360E+00 NEWF = 0.922492E-02

ITER = 2 X1 = 0.500000E+00 F(X1) = 0.203533E-01 X2 = 0.489360E+00 F(X2) = 0.922492E-02

NEWX = 0.480541E+00 NEWF = -.357777E-03

ITER = 3 X1 = 0.489360E+00 F(X1) = 0.922492E-02 X2 = 0.480541E+00 F(X2) = -.357777E-03

NEWX = 0.480870E+00 NEWF = 0.587106E-05

Section 2.17 Computer Programs 109

CONVERGED SOLUTION: X = 0.48086980E+00 F = 0.58710575E-05

Program 2.2 Bairstow’s method

The input data required are

N D degree of the polynomial NP1 D N C 1 NM1 D N � 1 ALP,BET D initial approximations to NÞ and Nþ, respectively MITER D maximum number of iterations permitted (usual value: 100) EPS D convergence criterion (usual value: 0.0001) A(I),I D 1; 2; : : : ; NP1 D coefficients of the polynomial (A(1): constant, A(NP1): coefficient of xŁŁN)

Illustrative Example

Find the roots of the equation

f .x/ D x4 � 8x3 C 18x2 � 6x C 6 D 0: Using N D 4; ALP D �2:0; BET D 1:0; MITER D 100; EPS D 10�4; A(1) D 6; A(2) D �6; A(3) D 18; A(4) D �8, and A(5) D 1, the output of Program 2.2 follows:

ITERATION NO: 1 B(I): 6.00000000 -6.00000000 18.00000000 -8.00000000 1.00000000

ITERATION NO: 2 B(I): 1.62283897 -1.62283993 14.47933960 -7.45454550 1.00000000

ITERATION NO: 3 B(I): 0.37398720 -0.37398624 12.77804947 -7.15104675 1.00000000

ITERATION NO: 4 B(I): 0.05259323 -0.05259418 12.12726974 -7.02532578 1.00000000

ITERATION NO: 5 B(I): 0.00181770 -0.00181770 12.00453663 -7.00090694 1.00000000

ITERATION NO: 6 B(I): 0.00000572 -0.00000477 12.00000286 -7.00000095 1.00000000

QUADRATIC EQUATION: SOLUTION IS 2.000000 1.000000

COEFFICIENTS OF REMAINING POLYNOMIAL 12.000003 -7.000001 1.000000

QUADRATIC EQUATION: SOLUTION IS 4.000001 3.000000

Program 2.3 Newton–Raphson method for solving a system of nonlinear equations:

The following input data are required

110 Chapter 2 Solution of Nonlinear Equations

N D number of equations NP1 D N C 1 X(1); X(2); : : : ; X(N) D initial (guess) values for the N unknown variables DX(1); DX(2); : : : ; DX(N) D increments of X(1),X(2),: : :,X(N) for use in finite difference derivatives (usual value: 0.01 for each) MITER D maximum number of iterations permitted (usual value: 50) EPS D convergence criterion (usual value: 0.00001) Subprogram to evaluate the functions fi E.X/ in the following form:

SUBROUTINE FUN(X,F,N) DIMENSION X(N),F(N) F(1)=... F(2)=... . . . F(N)=... RETURN END

Illustration (Example 2.19)

Using N D 2; X(1) D 2; X(2) D 4; DX(1) D DX(2) D 0:01; MITER D 50, and EPS D 10�5, the output of Program 2.3 is as follows:

ITER = 1 VALUES OF X (I): 0.20000000E+01 0.15000000E+02 VALUES OF F (I): 0.40000000E+01 -0.43890562E+01 NEW X (I): 0.12145693E+01 0.25270457E+01 NEW F (I): 0.28611388E+01 -0.18417971E+01

ITER = 2 VALUES OF X (I): 0.12145693E+01 0.28611388E+01 VALUES OF F (I): 0.25270457E+01 -0.18417971E+01 NEW X (I): 0.59310120E+00 0.22624547E+01 NEW F (I): 0.47047043E+00 -0.54713690E+00

ITER = 3 VALUES OF X (I): 0.59310120E+00 0.47047043E+00 VALUES OF F (I): 0.22624547E+01 -0.54713690E+00 NEW X (I): 0.27988195E+00 0.22411087E+01 NEW F (I): 0.10090208E+00 -0.81864953E-01

ITER = 4 VALUES OF X (I): 0.27988195E+00 0.10090208E+00 VALUES OF F (I): 0.22411087E+01 -0.81864953E-01 NEW X (I): 0.20796204E+00 0.22276907E+01 NEW F (I): 0.58541298E-02 -0.34757853E-02

ITER = 5 VALUES OF X (I): 0.20796204E+00 0.58541298E-02 VALUES OF F (I): 0.22276907E+01 -0.34757853E-02 NEW X (I): 0.20435423E+00 0.22267201E+01

Section 2.17 Computer Programs 111

NEW F (I): 0.42915344E-04 -0.12516975E-04 ITER = 6 VALUES OF X (I): 0.20435423E+00 0.42915344E-04

VALUES OF F (I): 0.22267201E+01 -0.12516975E-04 NEW X (I): 0.20433748E+00 0.22267120E+01 NEW F (I): 0.00000000E+00 -0.11920929E-06

SOLUTION CONVERGED IN 6 ITERATIONS SOLUTION IS: X(I) ARE: 0.20433748E+00 0.22267120E+01

F(I) ARE: 0.00000000E+00 -0.11920929E-06

2.17.2 C Programs

Program 2.4 Newton–Raphson method The following input data is to be given interactively by the user:

X1 D initial (guess) value of the root DX D increment of x (step size) to be used in finite difference formula (usual value: 0.01) MITER D maximum number of iterations permitted (usual value: 100) EPS D convergence criterion (usual value: 0.00001) Subprogram to evaluate the function f (x). Illustration (Example 2.8) The following is a program listing and output:

/* Newton-Raphson Method - C Program */ #include <stdio.h> #include <stdlib.h> #include <math.h>

/* X1, X2 - current and next x value F1, F2 - current and next f value DX - percentage increment of x I, MITER - number of iteration and maximum number FD1 - current derivative of f value */

void func1 (X, F) float X, *F; {

*F=1.5*X/pow((1.0+X*X),2)-0.65*atan(1.0/X)+0.65*X/(X*X+1.0); return;

} void func2 (X, DX, FD) float X, DX, *FD; {

float FR, FP, Y; void func1 (); func1 (X, &FR); Y=(1.0+DX)*X; if (fabs (X) > = 0.001)

112 Chapter 2 Solution of Nonlinear Equations

{ func1 (Y, &FP); *FD=(FP-FR)/(DX*X);

} else {

Y=DX; func1 (Y, &FP); *FD=(FP-FR)/DX;

} return;

} main () {

int I, MITER; float X1, X2, F1, F2, DX, EPS, FD1; void func1 (), func2 (); printf (‘‘Initial data for Newton-Raphson Method\n’’); printf (‘‘\nMaximum number of Iterations?\n’’); scanf (‘‘%d’’, &MITER); printf (‘‘Tolerance of result?\n’’); scanf (‘‘%f’’, &EPS); printf (‘‘Initial guess value of X?\n’’); scanf (‘‘%f’’, &X1); printf (‘‘Increamental percentage of X?\n’’); scanf (‘‘%f’’, &DX); F2=1.0E10; func1 (X1, &F1); printf (‘‘\nIte No. X(I) F(I) FD(I) X(I + 1) F(

I+1)\n’’); for (I = 1; (fabs (F2) >EPS) &&(I}<=MITER); I++) {

func2 (X1, DX, &FD1); X2 = X1 - F1/FD1; func1 (X2, &F2); printf (‘‘%d %12.6f %12.6f %12.6f %12.6f %12.6f\n’’, I,

X1, F1, FD1, X2, F2); X1=X2; F1=F2;

} if (I>MITER)

printf(‘‘\nSolution not converged in MITER iterations’’); else

printf(‘‘\nSolution converged in \%d iterations’’, I - 1); printf(‘‘\nSolution is: X = \%12.6F F = \%12.6F$’’, X1, F1);

}

Initial data for Newton-Raphson Method

Maximum number of Iterations? 20 Tolerance of result? 0.0001

References and Bibliography 113

Initial guess value of X? 0.5 Incremental percentage of X? 0.1

Ite No. X(I) F(I) FD(I) X(I+1) F(I+1) 1 0.500000 0.290353 4.149514 0.430027 0.161431 2 0.430027 0.161431 3.798167 0.387525 0.047333 3 0.387525 0.047333 3.609621 0.374412 0.012016 4 0.374412 0.012016 3.555357 0.371032 0.002904 5 0.371032 0.002904 3.541673 0.370212 0.000693 6 0.370212 0.000693 3.538372 0.370016 0.000165 7 0.370016 0.000165 3.537585 0.369970 0.000039

Solution converged in 7 iterations Solution is: X= 0.369970 F = 0.000039

Program 2.5 Fixed–point iteration method The following input data are required:

X1 D initial (guess) value of the root MITER D maximum number of iterations permitted (usual value: 100) EPS D convergence criterion (usual value: 0.00001)

Subprogram to evaluate g(x) and g’(x) where f .x/ D x � g.x/ D 0. Illustration (Example 2.15)

The output of Program 2.5 is follows:

SOLUTION BY FIXED POINT ITERATION METHOD

I X(I) F(I) GD(I) 1 0.10000000 0.60653551 0.16134542 2 0.60653551 0.47203515 0.43700819 3 0.47203515 0.48191504 0.43697774 4 0.48191504 0.48074364 0.43867080 5 0.48074364 0.48087844 0.43848403 6 0.48087844 0.48086287 0.43850571

Solution converged in 6 iterations Solution: X = 0.48087844 F = 0.48086287 GPRIME = 0.43850571

REFERENCES AND BIBLIOGRAPHY 2.1. S. S. Rao, Mechanical Vibrations, 3rd edition, Addison-Wesley, Reading, MA, 1995. 2.2. E. Volterra and E. C. Zachmanoglou, Dynamics of Vibrations, Charles E. Merrill

Books, Columbus, OH, 1965.

114 Chapter 2 Solution of Nonlinear Equations

2.3. J. R. Howell and R. O. Buckius, Fundamentals of Engineering Thermodynamics, 2d edition, McGraw-Hill, New York, 1992.

2.4. R. W. Fox and A. T. McDonald, Introduction to Fluid Mechanics, 4th edition, John Wiley, New York, 1992.

2.5. J. E. Shigley and C. R. Mischke, Mechanical Engineering Design, 5th edition, McGraw- Hill, New York, 1989.

2.6. R. L. Norton, Design of Machinery, McGraw-Hill, New York, 1992. 2.7. W. H. Hayt and J. E. Kemmerly, Engineering Circuit Analysis, McGraw-Hill, New

York, 1986. 2.8. A. Constantinides, Applied Numerical Methods with Personal Computers, McGraw-

Hill, New York, 1987. 2.9. B. Carnahan, H. A. Luther, and J. O. Wilkes, Applied Numerical Methods, John Wiley,

New York, 1969. 2.10. L. Bairstow, ‘‘Investigations Relating to the Stability of the Aeroplane,’’ Reports and

Memoranda No. 154, Advisory Committee of Aeronautics, 1914. 2.11. F. L. Hitchcock, ‘‘An Improvement on the G. C. D. Method for Complex Roots,’’ J.

Math and Phys, Vol. 23, pp. 69–74, 1944. 2.12. D. Muller, ‘‘A Method for Solving Algebraic Equations Using an Automatic Comput-

er,’’ Math Tables Aids Comput, Vol. 10, pp. 208–215, 1956. 2.13. W. L. Frank, ‘‘Finding Zeros of Arbitrary Functions,’’ J. Assoc Comput Math, Vol. 5,

pp. 154–165, 1958. 2.14. S. Nakamura, Applied Numerical Methods with Software, Prentice Hall, Englewood

Cliffs, NJ, 1991. 2.15. F. M. White, Heat and Mass Transfer, Addison-Wesley, Reading, MA, 1988. 2.16. M. G. Salvadori, Numerical Methods in Engineering, Prentice-Hall, Englewood Cliffs,

NJ, 1952. 2.17. K. E. Atkinson, An Introduction to Numerical Analysis, 2d edition, John Wiley, New

York, 1989. 2.18. A. Ralston and P. Rabinowitz, A First Course in Numerical Analysis, 2d edition,

McGraw-Hill, New York, 1978. 2.19. J. H. Wilkinson, ‘‘The Evaluation of Zeros of Ill-Conditioned Polynomials,’’ Numer.

Math., Vol. 1, pp. 150–180, 1959.

Review Questions 115

REVIEW QUESTIONS

The following questions along with corresponding answers are available in an interactive format at the Companion Web site at http://www.prenhall.com/rao.

2.1. Give short answers: 1. What is the difference between a polynomial equation and a transcendental

equation? 2. How many roots can an nth degree polynomial equation have? 3. How many roots can a transcendental equation have? 4. Can you use the incremental search method to find the complex roots of an

equation? 5. Can you use the Newton–Raphson method to find the roots of a nondifferentiable

function? 6. What is the difference between the secant and the regula falsi methods? 7. What is the basic principle used in Bairstow’s method? 8. What is the similarity between the secant and Muller’s methods? 9. What is the role of Taylor’s series expansion in the development of New-

ton–Raphson method? 10. What is the limitation of the fixed point iteration method in finding the roots of

an equation? 2.2. Answer true or false to each of the following:

1. A polynomial with real coefficients can have all complex roots. 2. A polynomial with complex coefficients can have all real roots. 3. The root of the equation f .x/ D x3 � sin x D 0 can be found by finding the

intersection point of the two graphs f1.x/ D x3 and f2 D sin x . 4. The graphical method of finding the roots of an equation is computationally

simple. 5. If f .xi / f .xiC1/ < 0, the root of f (x) is guaranteed to lie between xi and xiC1. 6. The incremental method can be used to find both real and complex roots of f (x). 7. The incremental method can find the root x1 even when f .x1/ touches the x axis

tangentially. 8. A complex root can be determined even by giving a real value as the initial guess

in Newton’s method. 9. The regula falsi method is the same as the method of linear interpolation.

(continued, page 116)

116 Chapter 2 Solution of Nonlinear Equations

10. The fixed point iteration method and the method of successive substitutions are the same.

11. The Newton–Raphson method requires the inversion of the n ð n Jacobian matrix.

12. If the initial guesses are real numbers, Muller’s method cannot find the complex roots.

2.3. Define the following: 1. Jacobian of fi .x1; x2; : : : ; xn/. 2. stagnation point 3. interval of uncertainty 4. deflated polynomial 5. bisection method

2.4. Fill in the blanks with the proper words: 1. If the graph of f (x) touches the x axis at x1, then x1 denotes a root.

(single/multiple) 2. The method requires two initial guesses to find a root. (Newton’s/secant) 3. The secant method is similar to the method. (bisection/Newton’s) 4. The regula falsi method is similar to the method in bracketing a root.

(bisection/secant). 5. Muller’s method can be considered to be more general than the method.

(regula falsi/Newton’s) 2.5. Match the following:

1. Þ cosh 100 Þ

� Þ � 1 D 0 (a) Polynomial equation 2. 2Þ5 � 3Þ3 C 5Þ2 C 4 D 0 (b) General nonlinear equation 3. 4x

2 C 20 ln x � x � 1 D 0 (c) Transcendental equation

2.6. Match the following:

1. Fixed point iteration method (a) requires no derivatives of f (x)

2. Newton’s method (b) requires the equation to be written in the form x D g.x/ to find the root

3. Bisection method (c) requires the first derivative of f (x)

4. Muller’s method (d) requires two initial guess values for the root

5. Secant method (e) requires three initial guess values for the root

(continued, page 117)

Review Questions 117

2.7. Match the following:

1. Regula falsi method (a) successive substitution method

2. Fixed point iteration method (b) accelerated step size can be used

3. Muller’s method (c) quadratic factors are found

4. Bairstow’s method (d) linear interpolation method

5. Newton–Raphson method (e) quadratic interpolation method

6. Incremental search method (f) intersection of chord and x axis is sought

7. Secant method (g) based on Taylor’s series expan- sion of f (x)

2.8. Select the most appropriate answer out of the multiple choices given: 1. The values of x and the corresponding values of the function, f .x/ D x3 � 4x C 2,

in the bisection method are x D 0:0; f .0/ D 2:0; x D 1:0; f .1:0/ D �1:0; x D 0:5; f .0:5/ D 0:125; and x D 0:75; f .0:75/ D �0:578125. The next value of x is (a) 0.625 (b) 0.875 (c) 1.250.

2. In the incremental search method using an accelerated step size, the function values are evaluated at x D 0:00, 0.05, 0.15, and 0.35. The next function value is to be evaluated at (a) 0.55 (b) 0.75 (c) 0.45.

3. The forward finite difference approximation of the derivative of f .x/ D x3�4xC2 at x D 1 with a step size (h) of 0.1 is given by (a) �0.99 (b) �1.29 (c) �0.69.

4. The backward finite difference approximation of the derivative of f .x/ D x3 � 4x C 2 at x D 1 with a step size (h) of 0.1 is given by (a) �0.99 (b) �1.29 (c) �0.69.

5. The central finite difference approximation of the derivative of f .x/ D x3�4x C2 at x D 1 with a step size (h) of 0.1 is given by (a) �0.99 (b) �1.29 (c) �0.69.

6. Newton’s iteration, xiC1 D xi � f .xi /f 0.xi / , applied to the function f .x/ D x3 � 4x C 2 with x1 D 1:0 gives x2 as (a) 0.0 (b) 2.0 (c) 1.0.

7. The secant method, xiC1 D xi � f .xi / xi �xi�1f .xi /� f .xi�1/ , applied to the function f .x/ D x3 � 4x C 2 with x1 D 0:0 and x2 D 1:0 gives x3 as (a) 0.6667 (b) 0.3333 (c) 1.3333.

(continued, page 118)

118 Chapter 2 Solution of Nonlinear Equations

8. The regula falsi method, xiC1 D xi � f .xi / xi �xi�1f .xi /� f .xi�1/ , applied to the function f .x/ D ex � sin x with x1 D 0:0 and x2 D 1:0 gives x3 as (a) 3.1405 (b) �2.1405 (c) �1.1405.

9. The fixed point iteration, xiC1 D g.xi /, applied to the function f .x/ D x3 �4x C2 or x D .x3 C 2/=4 with x1 D 0:0 gives x2 as (a) 0.5 (b) 2.0 (c) 1.0.

10. Newton’s method, xiC1 D xi � f .xi /f 0.xi / , applied to the function f .x/ D x2 � 3x C 3 with x1 D i gives x2 as (a) 2:0 � 3:0i (b) 0:9231 � 0:3846i (c) 0:9231 C 0:6154i .

11. The Jacobian of the functions f1 D x21 � x2 C 2 and f2 D x31 � x22 at x1 D 0 and x2 D 1 is given by

(a)

2 4 0 �1

0 �2

3 5 (b)

2 4 2 �1

3 0

3 5 (c)

2 4 0 �1

0 0

3 5 :

12. The Newton–Raphson method, applied to the functions f1 D x1 � x2 and f2 D x21 C x22 � 4 with x .1/1 D 0 and x .1/2 D 1 gives x .2/1 and x .2/2 as (a) 2.5, 1.5 (b) 2.5, 2.5 (c) �2:5; �0:5.

13. If x1 D 1 is a double root of the polynomial f .x/ D x4 � 7x3 C 17x2 � 17x C 6, the function f1.x/ that gives f .x/ D f1.x/.x � 1/2 is (a) x2 � 2x C 1 (b) x2 � 5x C 6 (c) x2 � 5x � 6.

Problems 119

PROBLEMS Section 2.1

2.1. For a curved beam subjected to bending, such as a crane hook lifting a load, the location of the neutral axis .rn/ is given by (see Fig. 2.26)

4rn.2R � p

4R2 � d2/ D d2; where R is the radius of the centroidal axis and d is the diameter of the cross section (assumed to be circular) of the curved beam. Use a graphical procedure to find the value of d for which rn D 4d when R D 10 in.

2.2. When two spheres in contact, such as ball bearings, transmit a load P , the normal stresses developed at a point along the z axis are given by (see Fig. 2.11)

¦x

pmax D ¦y

pmax D �

2 641 � � z

a

tan�1

1� z a

3 75 .1 C ¼/ C 1

2 �

1 C . z a

/2 ;

where ¦x and ¦y are the normal stresses in the x and y directions, respectively, pmax is the maximum pressure at the contact point, and a is the radius of the deformed flat circular surface around the contact point. Using a graphical procedure, find the value of . za / at which

¦x pmax

attains a value of �0.5. Assume the value of Poisson’s ratio .¼/ as 0.3.

2.3. The design of a mechanical component requires the maximum principal stress to be less than the material strength. For a component subjected to arbitrary loads, the principal stresses are given by the solution of the cubic equation

¦ 3 � I1¦ 2 C I2¦ � I3 D 0 (E1) where

I1 D ¦x C ¦y C ¦z; (E2)

R R

d

rn

(b)(a)

d

Figure 2.26 Crane hook.

120 Chapter 2 Solution of Nonlinear Equations

I2 D ¦x ¦y C ¦x¦z C ¦y¦z � − 2xy � − 2yz � − 2zx ; (E3) and

I3 D ¦x ¦y¦z C 2−xy−yz−zx � ¦x − 2yz � ¦y− 2zx � ¦z− 2xy; (E4) with ¦x ; ¦y; ¦z denoting the normal stresses acting along the x, y, z directions, and −xy; −yz; −zx indicating the shear stresses acting in the xy, yz, zx planes, respectively, as shown in Fig. 2.27. Find the principal stresses in a machine component for the following data using a graphical procedure: ¦x D 3 ksi, ¦y D 0; ¦z D 0; −xy D 1 ksi, −yz D 2 ksi, and −zx D 1 ksi.

Section 2.2

In Problems 2.4 to 2.21, find an approximate solution using a graphical procedure. 2.4. The stress induced in a column subjected to an eccentric load (as shown in

Fig. 2.28) is given by the secant formula [2.5]

P

A D Syc

1 C �ec

k2

sec

l

2k

r P

AE

! ;

where P is the axial load applied, A is the cross-sectional area of the column, Syc is the yield stress of the material in compression, e is the eccentricity of the load, l is the length of the column, E is Young’s modulus, c is the distance of the outermost fiber from the neutral axis of the column, and k is the radius of gyration of the cross section of the column. Find the value of PA for the following data: Syc D 40; 000 psi, eck2 D 0:2; lk D 50, and E D 30 ð 106 psi.

2.5. The specific volume of a gas (v) can be expressed as a function of its pressure p and temperature T , using van der Waals equation, as [2.3]

p D K T v � b �

a

v2 ; (E1)

y

z

x

x�

y�

z�

zx�

xz�

xy�

zy�

yz�

yx�

Figure 2.27 Stresses acting on an element.

Problems 121

PA

O P

e

y

�

x

l

Figure 2.28 Column subjected to eccentric load.

where K is the gas constant, and a and b are constants. Find the specific volume of methane gas at a pressure of 100 atm and temperature 0Ž F with a D 4787:4 ft4 � lbf/lbm2; b D 0:0496ft3= lbm, and K D 96:35 ft � lbf/lbm �Ž R.

2.6. If two plates are fastened by a bolted joint, the stiffness of the fastened members or plates .km/ is given by

km D 0:577³ Ed 2 ln

� 5 ²

0:577l C 0:5d 0:577l C 2:5d

¦� ;

where E is Young’s modulus of the fastened members, d is the diameter of the bolt, and l is the thickness of the fastened members. Find the value of d that corresponds to a value of km D 7:5 ð 106 lb/in when E D 12 ð 106 psi and l D 1:5 in.

2.7. The normal stress induced at the inner fiber of a torsional helical spring is given by

¦i D (

4C2 � C � 1 4C.C � 1/

) Mc

I ;

where I D ³d464 ; c D d2 ; C D Dd , M is the bending moment, D is the mean coil diameter, and d is the wire diameter. Find the value of C that corresponds to a stress of ¦i D 55 ð 103 psi, when M D 5 lb-in and d D 0:1 in.

2.8. The angular position of the output link .�4/ of a four-bar linkage corresponding to any specified angular position of the input link .�2/ can be computed using the Freudenstein’s equation [2.6]

k1 cos �4 � k2 cos �2 C k3 D cos.�2 � �4/; (E1) where

k1 D da ; k2 D d

c ; and k3 D a

2 � b2 C c2 C d2 2ac

: (E2)

122 Chapter 2 Solution of Nonlinear Equations

Input link (crank)

Coupler link

Fixed link (ground)

Output link

a

d

c

b

2� 4�

Figure 2.29 Four-bar linkage.

a, b, c and d are the link lengths shown in Fig. 2.29. Find the value of �4 for the following data: �2 D 30Ž; a D 1; b D 2; c D 4, and d D 5.

2.9. The natural frequencies of radial vibration of an elastic sphere .!/ are given by [2.2]

tan.ka/ ka

D 4G 4G � .½ C 2G/.ka/2 ; (E1)

where

½ D vE .1 C ¹/.1 � 2¹/ D Lame’s constant

and G D E2.1C¹/ D shear modulus, a D radius of the sphere, ² D density, and

k2 D

!2²

½ C 2G

! :

When v is equal to 0.25, which makes ½ equal to G, Eq. (E1) reduces to

tan ka ka

D 1 1 � 3

4 k2a2

:

Find the first natural frequency of vibration (value of !) of a sphere with the following data. a D 10 in, ¹ D 0:25; ²g D 0:28 lb/in3, and E D 30 ð 106 psi.

2.10. The variation of temperature with time along the thickness of a slab, which is suddenly immersed in a hot fluid, requires the solution of the transcendental equation

þ tan þ D h0l k

; (E1)

where h0 is the convection coefficient, k is the thermal conductivity, and 2 l is the thickness of the slab. (See Fig. 2.30.) Determine the smallest root (value of þ) of Eq. (E1) for a slab with h0lk D 0.

2.11. An open box is to be constructed from a rectangular piece of cardboard of size 3000 ð 5000. If four square pieces are to be cut from the corners of the cardboard

Problems 123

x

k

Fluid, ho

l l Figure 2.30 A slab of thickness 2l (infinite size in yz–plane).

x

x

x

x

x x

x

x

50 �

30 �

Figure 2.31 Open box from rectangular cardboard.

to obtain a height of x 00 (see Fig. 2.31), determine the value of x that gives a box of volume 3000 in3.

2.12. The natural frequencies of a 3 degree of freedom spring–mass system, shown in Fig. 2.32, are given by the equation

Þ3 � 5Þ2 C 6Þ � 1 D 0; (E1) where Þ D m!2k , m is the mass, ! is the natural frequency, and k is the spring stiffness. Find the roots of Eq. (E1).

2.13. The determination of the temperature variation with time in the radial direction of a sphere of radius R and thermal conductivity k, when suddenly immersed in a hot fluid with a heat transfer coefficient h0, requires the solution of the transcendental equation

1 � þ cot þ D h0 R k

: (E1)

Determine the first root of Eq. (E1) for a sphere with h0 Rk D 1:0.

124 Chapter 2 Solution of Nonlinear Equations

m

x1

k

m

x2

k

m

x3

k

Figure 2.32 Spring-mass system.

P

y(x)

y

x

l

x Figure 2.33 Transverse deflection of column.

2.14. The transverse deflection of a long column, fixed at one end and pinned at the other end (Fig. 2.33), under axial load is given by

y.x/ D 0:77y0

1 � cos ax l

C 1 a

sin ax

l � x

l

½ ;

where a D 4:49, l is the length of the column, and y0 is a constant. Find the value of

� x l

Ð at which y.x/ D 0:5 y0.

2.15. The velocity of a nonNewtonian fluid flowing in a circular tube (u) can be expressed as

u

umean D �

3n C 1 n C 1

�" 1 �

� r R

nC1 n

# ;

Problems 125

Length of plate v (Suction velocity)

Porous plate

u�

�

A

y

u(y)

B

Figure 2.34 Flow past a porous flat plate.

where umean is the mean velocity of the fluid, r is radial distance from the center of the tube, R is the radius of the tube, and n is a constant whose value depends on the fluid (for example, n D 1 for a Newtonian fluid, n D 3 for a dilatant fluid, and n D 13 for a pseudoplastic fluid). Determine the value of n of a fluid for which uumean D 0:9 at rR D 0:8.

2.16. When water flows steadily past a porous flat plate as shown in Fig. 2.34, the variation of velocity (u) with the distance y at section AB is given by [2.4]

u

u1 D 3

� y Ž

� 2

� y Ž

1:5 :

Find the value of y Ž

at which uu1 attains a value of 0.5.

2.17. The discharge coefficient (C) of an orifice (shown in Fig. 2.35) is defined as

C D actual mass flow rate theoretical mass flow rate

and is given by [2.4]

C D 0:5959 C 0:0312þ2:1 � 0:184þ8 C 91:71þ 2:5

Re0:75 ; (E1)

where þ denotes the ratio of diameters,

þ D d D

;

and Re indicates Reynold’s number. Find the value of þ that yields a discharge coefficient of C D 0:5959 at a Reynold’s number of Re D 104.

2.18. The ratio of stagnation pressure .p0/ to static pressure (p) in a compressible fluid flow can be expressed as [2.4]:

p0 p

D 1 C k 2

M2 ²

1 C 1 4

M2 C �

2 � k 24

� M4 ¦

;

where k is the ratio of specific heats of the fluid and M is the Mach number.

126 Chapter 2 Solution of Nonlinear Equations

D

d

Flow

Figure 2.35 Flow through orifice.

Find the value of the Mach number at which the ratio of pressures, �

p0 p

, will

be equal to 1.2 when the fluid is air for which k D 1:4. 2.19. A set of three mutually perpendicular axes in a body for which the products

of inertia about these axes are zero are known as the principal axes of inertia of the body. They are given by the solution of the following determinantal equation: þþþþþþþþþ

Ixx � I Ixy Ixz Iyx Iyy � I Iyz Izx Izy Izz � I

þþþþþþþþþ D 0:

Find the principal moments of inertia of a body for which Ixx D 12; Iyy D 16; Izz D 16; Ixy D 6; Iyz D 2, and Ixz D �6. Hint: Ii j D Iji for any body.

2.20. The vertical velocity of a motorcycle due to a road bump is given by (see Fig. 2.36)

v.t/ D Xe��!n t .��!n sin !d t C !d cos !d t/; where X is the maximum displacement of the motor cycle, � is the damping constant of the suspension, !n is the undamped natural frequency of the system, !d D !n

p 1 � � 2 is the damped natural frequency of the system, and

t is time. Determine the time, t , at which the velocity of the motorcycle attains a value of 1 m/s for the following data: X D 0:4550 m; � D 0:4037, and !n D 3:4338 rad/s.

2.21. The schematic diagram of a large cannon is shown in Fig. 2.37 [2.1]. When the gun is fired, high-pressure gases accelerate the projectile inside the barrel to a high velocity. The reaction force pushes the gun barrel in the opposite direction of the projectile. Since it is desirable to bring the gun barrel to rest in the shortest time without oscillation, it is made to translate backward against a critically damped spring–damper system called the recoil mechanism. In a particular case, the recoil distance (x) at time t is given by

x D vte�!n t ;

Problems 127

Figure 2.36 Motor cycle over a road bump.

Gun barrel

Recoil mechanism (spring and damper)

Projectile

Figure 2.37 Cannon with a projectile.

where v is the initial recoil velocity, !n is the natural frequency of the recoil mechanism, and t is the time. Find the time corresponding to a recoil distance of 0.1 m when v D 4:8626 m/s and !n D 4:4721 rad/s.

Section 2.3

2.22. Find the root of the equation (see Example 2.2)

f .x/ D 1:1288p x

� tanh.1:3288px/ D 0

using incremental search method with x1 D 1:0 and 1x .1/ D 0:1.

128 Chapter 2 Solution of Nonlinear Equations

2.23. Find the root of the equation (see Problem 2.3)

f .x/ D x3 � 3x2 � 6x C 8 D 0

using incremental search method with x1 D 0:0 and 1x .1/ D 0:5. 2.24. Find the interval a x b such that f (a) and f (b) have opposite signs for

the function f .x/ D e2x � 3x � 8

using incremental search procedure. 2.25. Find the solution of Problem 2.15 using incremental search method with

n1 D 0:1 and 1n.1/ D 0:1.

Section 2.4

2.26. Find the root of the equation (see Example 2.3)

f .x/ D p

4x2 � 0:04 C 0:24 sin�1 �

0:1 x

� � 10:055752 D 0

using bisection method with a D 2; b D 10, and " D 0:001. 2.27. Find the solution of Problem 2.1 using bisection method with a D 0:0; b D 5:0,

and " D 0:001. 2.28. Find the solution of Problem 2.7 using bisection method with the approxima-

tions 5.0 and 12.0 for C , and " D 0:001.

Section 2.5

2.29. Find the roots of the equation

f .x/ D x3 � 0:2589x2 C 0:02262x � 0:001122 D 0

using Newton–Raphson method with x1 D 0:0 and " D 0:001. 2.30. Find the solution of Problem 2.2 using Newton–Raphson method with the

starting point � z

a

1

D 0:0 and " D 0:001. 2.31. Find the solution of Problem 2.8 using Newton–Raphson method with the

starting value of �4 as 90Ž and " D 0:001.

Problems 129

yn

� Figure 2.38 Flow through trapezoidal channel.

2.32. Find the solution of Example 2.5 using Newton–Raphson method with the starting value of v D 0:0 and " D 0:001.

Section 2.6

2.33. Find the root of the equation f .x/ D 8x3 C x2 C 8x � 3 D 0 using secant method with 0.0 and 0.6 as the initial approximations for x and " D 0:001.

2.34. Find the solution of Problem 2.6 using secant method with 0.0 and 1.0 as the initial approximations for d and " D 0:001.

2.35. When water flows through a trapezoidal channel section (Fig. 2.38), the cross- sectional area of the fluid (A) and the wetted perimeter of the fluid (P) are related as

P D A yn

� yn cot � C 2ynsin � :

Find the value of � for which P D 3 m, when A D 1 m2 and yn D 0:5 m with 0Ž and 60Ž as the initial approximations for � and " D 0:001.

Section 2.7

2.36. Find the solution of the equation f .½/ D 0:8 � tan ½ �

1�0:08 tan ½ 0:08Ctan ½

D 0 using

regula falsi method with 10Ž and 30Ž as the initial approximations for ½ and " D 0:001.

2.37. Find the solution of Problem 2.16 using regula falsi method with the initial approximations 0.0 and 1.0 for y

Ž and " D 0:001.

2.38. Find the root of the equation f .x/ D ln x with the initial approximations x1 D 0:6 and x2 D 6:0 using (a) secant method and (b) regula falsi method. Discuss your observations.

130 Chapter 2 Solution of Nonlinear Equations

Section 2.8

2.39. Find the root of the equation

f .x/ D x � 1 x2 C 9 D 0

using fixed point iteration with x1 D 3:0 and " D 0:001. 2.40. Find the root of the equation

f .x/ D x4 � 31x3 C 305x2 � 1025x C 750 D 0 using fixed point iteration with x1 D 20:0 and " D 0:001.

2.41. A person plans to deposit $ 100 a month into a bank account for 10 years. Find the interest rate needed to accumulate a total amount of $ 20,000, including interest at the end of 10 years. Use fixed-point iteration by taking the initial value of the interest rate as 0.05 and " D 0:001. Hint: Formula is S D R

n .1Ci/n�1

i

o , where S D future sum, R D uniform

payments, n D number of payments, and i D interest rate. 2.42. The fixed point iteration method can be extended to find the solution of several

nonlinear equations. If f1.x; y/ D 0 and f2.x; y/ D 0 denote the equations to be solved, they are rewritten in the form:

f1.x; y/ D 0 �! x D g1.x; y/ and

f2.x; y/ D 0 �! y D g2.x; y/: If .x1; y1/ denotes the initial approximation to the solution, then the fixed-point iteration can be used to find the next approximation as

xiC1 D g1.xi ; yi / and

yiC1 D g2.xi ; yi /: Using this technique, solve the following system of equations:

y2 D 8x and

.x � 1/2 C .y � 2/2 � 4 D 0:

Section 2.9

2.43. Find the root of multiplicity two of the following equation using the starting point, x1 D 1:0:

f .x/ D x4 � 27x3 C 265x2 � 1125x C 1750 D 0:

Problems 131

2.44. Find the root of multiplicity three of the following equation using the starting point, x1 D 2:0:

f .x/ D x4 � 8x3 C 18x2 � 16x C 5 D 0: 2.45. Find the root of muliplicity four of the following equation using x1 D 0:0:

f .x/ D x4 � 8x3 C 24x2 � 32x C 16 D 0:

Section 2.10

2.46. Find the roots of the equation

x4 � 32x3 C 244x2 � 20x � 1200 D 0

using Bairstow’s method. Assume the initial approximations of NÞ and Nþ as �1 and 1, respectively.

2.47. Find the roots of the equation

x3 � 7x2 C 12x � 10 D 0

using Bairstow’s method. Assume the initial values of NÞ and Nþ as �1 and 1, respectively.

Section 2.11

2.48. Solve Problem 2.17 using Muller’s method with the starting values þ1 D 0:0; þ2 D 0:2, and þ3 D 0:8.

2.49. Use Muller’s method to find a root of the equation

f .x/ D ex � x2 D 0

using the initial approximations x1 D �1; x2 D 0, and x3 D 1. 2.50. Use Muller’s method to find a root of the equation

f .x/ D x4 � 31x3 C 305x2 � 1025x C 750 D 0

using the starting points x1 D 4; x2 D 6, and x3 D 12.

Section 2.12

2.51. Find the Jacobian matrix for the following functions at the points (i) x D 1; y D 1, and (ii) x D 0; y D 0:

f1.x; y/ D 10x3 � y and f2.x; y/ D xy � 1:

132 Chapter 2 Solution of Nonlinear Equations

2.52. Find the solution of the following equations using Newton–Raphson method with the starting point (0.0, 0.0):

10x3 � y D 0 and xy � 1 D 0: 2.53. Find the solution of the following equations using Newton–Raphson method

with the starting point (1.0, 1.0):

cosh x � y D 0 and xy � 1 D 0: 2.54. Find the solution of the following equations using Newton–Raphson method

with the starting point (0.1, 0.1):

10x2 C y2 � 10 D 0 and x2 C y2 � 4x � 2y C 1 D 0:

Section 2.13

2.55. Consider the simultaneous equations:

f1.x; y/ D x2 � y D 0 (E1) and

f2.x; y/ D y2 � x D 0: (E2) Construct a function F(x, y) whose minimum will yield the solution of Eqs. (E1) and (E2). Derive the necessary and sufficient conditions corresponding to the minimum of F(x, y).

Section 2.14

2.56. Prove that the results of Example 2.7 satisfy the convergence relationship given by Eq. (2.91).

2.57. Show that the secant method converges with an order of convergence 1.618 using the results of Example 2.14.

2.58. Show that the fixed point iteration converges linearly using the results of Example 2.15.

Section 2.16

2.59. Find the solution of Problem 2.36 using Maple. 2.60. Solve Problem 2.18 using Maple. 2.61. Find the solution of the equations f1.x; y/ D 2ex C xy � 2 D 0 and f2.x; y/ D

sin.x y/ C 2x C y � 1 D 0 using Maple. 2.62. Solve Problem 2.36 using Matlab. 2.63. Solve Problem 2.18 using Matlab.

Problems 133

2.64. Find the solution of the equations given in Problem 2.61 using Matlab. 2.65. Solve Problem 2.36 using Mathcad. 2.66. Solve Problem 2.18 using Mathcad. 2.67. Find the solution of the equations given in Problem 2.61 using Mathcad.

Section 2.17

2.68. Find the solution of Problem 2.8 using PROGRAM 2.1 (Secant method) with X1 D 0:0, X2 D 5:0, MITER D 100, and EPS D 0:001.

2.69. Solve Problem 2.18 using PROGRAM 2.1 (Sécant method) with X1 D 0:1, X2 D 2:0, MITER D 100, and EPS D 0:0001.

2.70. Solve Problem 2.40 using PROGRAM 2.2 (Bairstow’s method). Assume the initial values of NÞ and Nþ as 0.0 and 2.0 respectively.

2.71. Solve Problem 2.44 using PROGRAM 2.2 (Bairstow’s method) with NÞ D 1:0 and Nþ D 2:0.

2.72. Using PROGRAM 2.3 (Newton–Raphson method), find the solution of the following equations:

cosh x � y D 0 and xy � 1 D 0: Assume the initial values of x and y as 1.0 and 1.0, respectively.

2.73. Solve Problem 2.61 using PROGRAM 2.3. 2.74. Find the solution of Problem 2.19 using PROGRAM 2.4 (Newton–Raphson

method) with X1 D 0:1; DX D 0:01; MITER D 100, and EPS D 0:001. 2.75. Find the solution of Problem 2.26 using PROGRAM 2.4 (Newton–Raphson

method) with X1 D 0:1; DX D 0:01; MITER D 100, and EPS D 0:001. 2.76. Solve Problem 2.14 using PROGRAM 2.5 (Fixed point iteration) with X1 D

0:1 MITER D 100, and EPS D 0:0001. 2.77. Solve Problem 2.22 using PROGRAM 2.5 (Fixed point iteration) with X1 D

0:01 for þ; MITER D 100, and EPS D 0:001.

General

2.78. Write a computer program to implement the incremental search method for finding the root of the equation, f .x/ D 0, using x1 as the starting point and 1x .1/ as the initial step length. Use this program to find the root of the equation

f .x/ D 1:5x .1 C x2/2 � 0:65 tan

�1 �

1 x

� C 0:65x

1 C x2 D 0

using x1 D 0:0 and 1x .1/ D 0:1.

134 Chapter 2 Solution of Nonlinear Equations

x

P

P

l

y

z

2b

Figure 2.39 Contact stress between cylinders.

2.79. Write a computer program to implement the regula falsi method for finding the root of the equation, f .x/ D 0, using x1 and x2 as the starting points. Using this program, find the root of the function given in Problem 2.78 with x1 D 0:0 and x2 D 1:0.

2.80. Write a computer program to find the root of the equation, f .x/ D 0, using Muller’s method. Using this program, find the roots of the polynomial equation

f .x/ D x4 � 12x3 C 49x2 � 78x C 40 D 0:

2.81. When two cylinders in contact, such as roller bearings, transmit a load P as shown in Fig. 2.39, the normal stress developed at a point along the z-axis is given by

¦y

pmax D �

0 B@2 � 1

1 C � z

b

2 1 CA r

1 C � z

b

2 C 2 � z b

;

where pmax is the maximum pressure along the contact line and 2b is the width of the deformed flat rectangular surface around the contact line. Find the value of

� z b

Ð at which ¦ypmax D �0:5 using a graphical procedure.

2.82. The efficiency of a square-threaded screw .�/ is given by

� D tan ½ �

1 � ¼ tan ½ tan ½ C ¼

� ;

where ¼ is the coefficient of friction and ½ is the lead angle of the threads. Find the lead angle ½ which gives an efficiency of 0.8 when ¼ D 0:08.

Problems 135

2.83. The natural frequencies of vibration of a circular plate with a clamped edge .½/ are given by [2.2]

tan � ½a � ³

4 .n C 1/

D

1 � Q .8½a/2

C R .8½a/3

� S .8½a/4

1 � P .8½a/

C Q .8½a/2

� S .8½a/4

;

with

P D 2.4n2 � 1/; Q D .4n2 � 1/.4n2 � 9/;

R D 2 3

.4n2 � 1/.4n2 � 9/.4n2 � 13/;

and

S D 1 6

.4n2 � 1/.4n2 � 9/.4n2 � 25/;

where ¹ is the Poisson’s ratio, a is the radius of the plate, and n is the number of nodal lines that define the mode shape. Find the natural frequency (value of ½a) corresponding to n D 1 for a plate with ¹ D 0:3.

2.84. The thermal efficiency of a uniform fin .�/, considering convection through the tip, is given by [2.15]

� D

0 B@ sinh

l

½ C Þ cosh l

½

cosh l

½ C Þ sinh l

½

1 CA 0 B@ ½

l C A P

1 CA ;

with

½ D �

k A

h1 P

� 1 2

and

Þ D �

h1 A k P

� 1 2

;

where l D length; A D cross-sectional area, P D perimeter of the cross section, k D thermal conductivity, and h1 D heat transfer coefficient of the fin (Fig. 2.40). If the fin is made of aluminum with a square cross section with l D 0:1m; k D 240W=.m �Ž C/, and h1 D 9W=.m2 �Ž C/, determine the necessary cross-sectional dimensions of the fin to achieve an efficiency of 0.95.

136 Chapter 2 Solution of Nonlinear Equations

l

Tip

Figure 2.40 Uniform fin.

2.85. The equation governing the natural frequencies of a free-free beam with continuous material distribution is given by [2.1]

cos þnl cosh þnl D 1; where the nth natural frequency, !n , is given by

!n D .þnl/2 �

E I

² Al4

� 1 2

and n D 1; 2; : : : ;

where l is the length, E is Young’s modulus, I is the area moment of inertia of the cross section, ² is the density, and A is the area of cross section of the beam. Find the fundamental (smallest) natural frequency of vibration of the beam, þ1l.

2.86. The pressure-volume (p-v) diagram of an ideal air-standard dual-cycle heat engine is shown in Fig. 2.41. In this engine, the piston starts compressing the air isentropically at state 1 and continues to state 2. Fuel is injected into the air at state 2. Because of the high compression ratio used, the temperature of the compressed air will be sufficiently high to ignite the fuel without any ignition source. Thus, combustion occurs during the constant volume process of fuel injection and continues at constant pressure while the piston moves to reach state 4. The combustion process ends at state 4, and expansion starts at point 4 and continues to point 5 isentropically. At state 5, the exhaust valve opens and the pressure drops to the initial pressure at state 1. The cycle then repeats. The efficiency of the engine .�/ can be expressed as [2.3]

� D 1 �

" 1

rk�11

#" p3 p2

� r1 r2

�k � 1

#

p3 p2

� 1 C k �

p3 p2

�²� r1 r2

� � 1

¦ ; where r1 is the volume compression ratio,

r1 D v1

v2 ;

r2 is the volume expansion ratio

r2 D v5 v4

;

Problems 137

Pressure (p)

3 4

2

5

1

Isentropic compression

Isentropic expansion

Volume (v) Figure 2.41 p � v diagram of the dual cycle.

p2 and p3 are the pressures at states 2 and 3, respectively, and k is the ratio of specific heats of air. If the efficiency of a dual-cycle heat engine is known to be 0.651, find the value of r1 when k D 1:4; r2 D 11:82, and p3p2 D 1:3657.

2.87. The speed of wave propagation of a surface wave in a fluid flowing in a channel is given by [2.4].

v D s

g½

2³ tanh

� 2³y

½

� ;

where v is the speed of the wave, ½ is the length of the wave, y is the depth of the fluid and g is the acceleration due to gravity. Find the value of ½y for which

vp gy D 0:99 assuming the fluid to be water.

2.88. Consider the simultaneous equations

f1.x; y/ D y2 � 2y � 8x � 23 D 0 (E1) and

f2.x; y/ D y2 C 4x2 C 4y � 8x � 8 D 0 (E2) (a) Plot the functions f1.x; y/ D 0 and f2.x; y/ D 0 and find an approximate

solution of the problem. Using the approximate solution as a starting point, find the solution of Eqs. (E1) and (E2) by applying the Newton–Raphson method.

(b) Construct a function F.x; y/ whose minimum will yield the solution of Eqs. (E1) and (E2). Derive the necessary and sufficient conditions corresponding to the minimum of F.x; y/.

138 Chapter 2 Solution of Nonlinear Equations

PROJECTS 2.1. The thermal efficiency of a uniform fin .�/, with an insulated tip, is given

by [2.15]

� D tanh l ½

l ½

with

½ D �

k A

h1 P

� 1 2

;

where l D length; A D cross-sectional area, P D perimeter of the cross section, k D thermal conductivity, and h1 D heat transfer coefficient of the fin. If the fin is made of aluminum with a square cross section with l D 0:1 m, k D 240 W=.m-ŽC/, and h1 D 9 W=.m2-ŽC/, determine the necessary cross- sectional dimensions of the fin to achieve an efficiency of 0.95 using the following methods:

(a) Graphical method (b) Bisection method (c) Newton–Raphson method (d) Regula Falsi method (e) Muller’s method

2.2. The natural frequencies of a beam, fixed at one end .x D 0/, and carrying a mass M at the other end .x D l/, are given by the frequency equation

1 C 1 cos þnl cosh þnl

� Þþnl.tan þnl � tanh þnl/ D 0;

where

Þ D mass attached at the end of the beam mass of the beam

D M ² Al

is the mass ratio,

þn D p

!n

l

( ² Al4

E I

) 1 4

;

l is the length, E is the Young’s modulus I is the area moment of inertia of the cross section, ² is the density and A is the cross-sectional area of the beam. (a) Determine the fundamental frequency of vibration of the beam, þ1l, for

Þ D 0:0, 1.0, 10.0, 100.0, and 1000.0. (b) Plot a graph to show the variation of þ1l with Þ.

Projects 139

r1 r5

r2

r3

r4

Straight line generated

�

1�

2�

3� P

Figure 2.42 Four - bar straight - line mechanism.

(c) Determine how each of the parameters of the beam need to be changed to maximize the fundamental frequency of vibration of the beam.

2.3. A four-bar straight-line mechanism is shown in Fig. 2.42. For specified values of the input link length r1 and the input angle �1, the coupler angle �2 and the follower angle �3 can be computed by solving the following vector loop equations:

r1 cos �1 C r2 cos �2 � r3 cos �3 � r4 D 0 and

r1 sin �1 C r2 sin �2 � r3 sin �3 D 0: (E1) The x and y components of the position of the coupler point P can be determined as

xP D r1 cos �1 C r5 cos.�2 C �/ and

yP D r1 sin �1 C r5 sin.�2 C �/: (E2) Solve Eqs. (E1) using Newton–Raphson method, find the (x, y) coordinates of P using Eqs. (E2), and plot the position of P for the following data:

r1 D 0:963; r2 D 0:764; r3 D 0:528; r4 D 1:815; r5 D 0:778; � D �89:65Ž; and 0 �1 90Ž:

2.4. Discuss the influence of small errors in the coefficients ai of the polynomial equation

f .x/ D an xn C an�1xn�1 C Ð Ð Ð C a1x C a0 D 0: on the accuracy of the computed roots of the equation. Note: If the new coefficients are given by bi D ai C "i ; i D 0; 1; 2; : : : n, determine the new roots yi D xi C Ži ; i D 1; 2; : : : ; n, where "i and Ži are small numbers and xi are the roots of the original polynomial.

140 Chapter 2 Solution of Nonlinear Equations

(a) Find the roots of the equation

f .x/ D a2x2 C a1x C a0 D 0 using a2 D 1; a1 D 5; a0 D 6:24; "2 D 0:01; "1 D �0:03, and "0 D 0:02.

(b) Find the roots of the equation [2.19]

f .x/ D 20Y

iD1 .x � i/ C 2�23x19 D 0

and discuss the results.