Pao Y.C. Engineering Analysis: Interactive Methods and Programs with FORTRAN, QuickBASIC, MATLAB, and Mathematica

Подождите немного. Документ загружается.

Cramer’s rule can then be applied to obtain x and y as

(22)

and

(23)

where u, v, and their derivatives are to be evaluated at (x

). Equations 22 and 23

are to be continuously applied to adjust the guessing values of (x

) until both

u(x

) and v(x

) are negligibly small.

Program NewRaphG has been developed by use of the iterative equations 22

and 23. Both QuickBASIC and FORTRAN versions of this program are listed

below along with a sample application of solving the intercepts of two ellipses,

f(x,y) = (x/3)

+ (y/4)

–1 and g(x,y) = (x/4)

+ (y/3)

–1.

QUICKBASIC VERSION

FORTRAN VERSION

∆xuvvuuvuv

yyxyyx

=− +

()

−

()

,,,,,,

∆yuvvuuvuv

xxxyyx

=+ −

()

−

()

,,,,,,

Sample Application

MATLAB APPLICATIONS

Here, we provide a m ﬁle called NewRaphG.m as a companion of the FOR-

TRAN and QuickBASIC versions:

For using this function, the problem to be solved needs to be deﬁned by creating

two m ﬁles, in which the equations involved and the expressions for their ﬁrst

derivatives should be spelled out. In case of solving the sample problem used in

FORTRAN and QuickBASIC versions, ﬁrst we may deﬁne the equation as:

Next, the expressions for their ﬁrst derivatives may then be written as:

To solve this problem, the interactive application of MATLAB proceeds as

follows (some displays have been rearranged for saving spaces):

Notice that the initial values are taken as X(1) = –2 and X(2) = 2, a tolerance

of 10

–5

and the iteration is limited to 20 trials. The solutions are found after four

Newton-Raphson trials when the sum of the absolute values of the two equations is

equal to 3.7361x10

–8

MATHEMATICA APPLICATIONS

Mathematica applies the Newton’s method in its function FindRoot which can

be applied for solving a polynomial, or, transcendental equation, and also for multiple

equations. We illustrate its applications by using the examples discussed earlier.

First, the root near X = 1.75 for a third-order polynomial is found:

In[1]: = FindRoot[{X^3–6*X^2 + 11*X–6 = = 0}, {X,1.75}]

Out[1] = {X -> 2.}

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

a transcendental equation described inside the ﬁrst pair of braces:

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

Out[2] = {T -> 1.09911}

For solving two simultaneous transcendental equations, two examples are pre-

sented below. The ﬁrst is to ﬁnd one of the intercepts of two ellipses and the second

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

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

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

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

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

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

3.4 PROGRAM BAIRSTOW — BAIRSTOW’S METHOD FOR

FINDING POLYNOMIAL ROOTS

Program Bairstow is developed for ﬁnding the roots of polynomials based on the

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

Let a nth-order polynomial be denoted as:

(1)

Notice that the highest term x

has a coefﬁcient equal to 1; otherwise the entire

equation must be normalized by dividing by that coefﬁcient. The Bairstow’s method

consists of ﬁrst selecting a trial divider D(x) = x

+ d

x + d

, and to obtain the

quotient Q(x) = x

N–2

+ q

N–3

+ q

N- 4

+ ••• + q

N–4

+ q

N–3

x + q

N–2

and a remainder

R(x) = r

x + r

. The objective is to continuously adjust the values of d

and d

until

both values of r

and r

are sufﬁciently small. It is apparent that both r

and r

are

dependent of d

and d

. Taylor’s series expansions of r

and r

can be written as:

(2)

and

(3)

where

and so on.

The adjustments d

and d

are to be calculated so as to make the left-hand

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

be small enough (if the guessed values of d

and d

values are sufﬁciently close to

those which make both r

and r

equal to zero) so that the second and higher derivative

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

(4)

and

(5)

Px x ax ax a x a x a

NN N

NNN

()

=+ + +…+ + +

−−

−−1

rd dd d rdd r dd d

rddd

11 12 2 112

12 1

12 2

()

+…

∆∆ ∆

∆

,,,

rd dd d rdd r dd d

rddd

21 12 2 212

12 1

12 2

()

+…

∆∆ ∆

∆

,,,

rrdr rd

dd1

, ≡∂ ∂ ≡∂ ∂

rdddrddd rdd

dd1

12 1

12 2 112

,,,

()

=−

()

∆∆

rdddrddd rdd

dd2

12 1

12 2 212

,,,

()

=−

()

∆∆

Cramer’s rule can then be applied to obtain d

and d

as:

(6)

and

(7)

where r

, r

, and their partial derivatives are to be evaluated at (d

). Equations 6

and 7 are to be continuously applied to adjust the guessing values of (d

) until

both r

) and r

) are negligibly small.

To calculate the adjustments d

and d

based on Equations 6 and 7, we need

to ﬁnd the partial derivatives ∂r

/∂d

, ∂r

/∂d

, ∂r

/∂d

, and ∂r

/∂d

. These derivatives

are, however, depend on the d

and d

, and the coefﬁcients q’s in the quotient Q(x).

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

are:

(8,9)

and

(10)

It can also be shown that the coefﬁcients in the remainder R(x) are:

(11,12)

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

extended to N and if the remainder is redeﬁned as:

(13)

That is, r

is renamed as q

N–1

and r

is equal to d

N–1

+ q

. As a consequence,

we need to replace r

and r

in Equations 6 and 7 by q

N–1

and q

. For calculation of

the adjustments d

and d

, Equation 10 should be used for q

N–1

and q

and to

derive their partial derivatives respect to d

and d

. Since all q’s are functions of d

and d

, to derive the partial derivatives of the last two q’s we must ﬁnd the partial

derivatives for all q’s starting with q

. From Equations 8 to 10, we can have:

(14,15)

(16)

∆drrrrrrrr

dddddd

11212

221221

=− +

()

−

()

,,,,,,

∆drrrrrrrr

dddddd

11212

1 1 12 21

=− +

()

−

()

,,,,,,

qad q aqdd

111 2 2112

=− =− − and

qaqdqd k N

kkk k

=− − = …−

−−1

34 2, , , ,for

ra qdqd ra qd

NN N NN112132 2 22

=− − =−

−− − −

and

Rx x d q q

()

−11

∂∂=− ∂∂=qd qd

11 12

10, ,

∂∂=−∂∂

()

−=−q d qddqdq

2 1 1111 11

(17)

and for k = 3,4,…,N

(18)

(19)

It can be concluded from the above results that:

(20)

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

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

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

changed to b’s, d

and d

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

to represent the derivatives of q’s.)

(1) Specify the values of N, a

through a

, and a tolerance .

(2) Assume an initial guessing values for d

and d

for the divider D(x).

(3) Calculate the coefﬁcients q

through q

N–2

for the quotient Q(x) using

Equations 8 to 10.

(4) Also use Equation 10 to calculate the coefﬁcients q

N–1

and q

for the

remainder R(x).

(5) Test the absolute values of q

N–1

and q

. If they are both less than , two

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

order of P(x), N, is to be reduced by 2, and q

through q

N–2

are to become

through a

N–2

, respectively, and return to Step 2. This looping continues

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

can be calculated.

(6) If the absolute value of either q

N–1

or q

is greater than ε, calculate the

partial derivatives of q

with respect to d

, c’s using Equations 14, 16, and

18 for k = 3,4,…,N. The derivatives of q

with respect to d

are already

available due to Equation 20.

(7) Use Equations 6 and 7 to calculate the adjustments d

and d

, noticing

that r

and r

are to be replaced by q

N–1

and q

, respectively. The iteration

is resumed by returning to Step 3.

Both QuickBASIC and FORTRAN versions of the program Bairstow coded

following the steps described above are to be presented next.

∂∂=−∂∂

()

−=−qd qdq

22 121

∂∂=−∂ ∂

()

−−∂ ∂

()

−−−

qd q ddq q dd

kk kk

∂∂=−∂∂

()

−−∂∂

()

+−−

qd qddq qdd

kkkk1

221

∂∂=∂∂ =…−

qdqd k N

kk1

12 1 , , ,for

QUICKBASIC VERSION

Sample Application

As an example, the polynomial P(x) = x

4–

+ 13x

2–

19x + 10 = 0 is solved by

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

screen is:

The quotient in this case is a quadratic equation:

FORTRAN VERSION

Qx x i x i x i x x

()

=−+

()

[]

−−

()

[]

=−

()

−

()

=−+=12 12 1 2 2 5 0