To Find All Roots of a Regular Polynomial

A regular polynomial is one that contains only integer powers of x. The Bairstow (or Bairstow-Lin) method finds all roots, both real and imaginary, of a regular polynomial with real coefficients. The method involves the successive extraction of quadratic factors from the original polynomial of degree n and subsequent reduced polynomials of degree n-2, AM and so on. The quadratic formula is then used to obtain pairs of roots, either real or complex, from the quadratic factors. If the degree of the polynomial is odd, then the remainder, after extracting quadratic factors, will be a linear factor, yielding the final root directly.

The calculation proceeds as follows. For the polynomial y = a„xn + £V,Jtn-1 + ... + ajx + aQ (8-8)

performing synthetic division by a trial quadratic x2+px + q (8-9)

yields a quotient and a remainder.

y = (x2 + px + q) (bnx"'2 + b„.jx"-3 + ... + b2) + (Rx + S) (8-10)

If {x2 + px + q) is an exact divisor, then the remainder (Rx + S) will be zero. Our task therefore is to find the values of p and q that make (Rx + S) equal to zero. This will make (x2 + rx + s) a quadratic factor of the polynomial.

Examination of the process of synthetic division reveals that there is a correspondence between the coefficients of the two preceding forms of the polynomial:

bn.k = an_k -pb„-k+\ - qbn-k+2 (k = 2,3, ... ,n-\) (8-14)

If the polynomial has been normalized so that a„ = 1, then the equations are simplified somewhat.

The trial quadratic will be a factor of the polynomial if the remainder is zero, that is, R = S = 0. Since R and S are functions of p and q:

we need to find the values of p and q that make R and S equal to zero. We will do this by means of a two-dimensional analog of the Newton-Raphson method. If p* and q* are the desired solution, then the solution can be expressed as a Taylor series

op Sq and

ignoring terms other that the first, since as we approach the correct answer the higher terms become negligible. The preceding result in two equations in two unknowns, which can be solved to obtain

To find the partial derivatives 5R/bp, etc, we could follow the usual procedure of making a small change in p to find the corresponding change in b. Instead, we will calculate the partial derivatives using analytical expressions. Differentiating the expressions 8-11 to 8-14 with respect to p yields the following:

dp " dp db„_7 , db„i db„ cn-2=~~ = -b -p^^-q—?- (8-27)

dp dp dp dbn db7

op dp

Equations 8-25 to 8-29 can be written in the form cn=0 (8-30)

The simultaneous equations to be solved are c2Ap + c3Aq = -¿i

A q:

-6,

c3

-b0

c2

c2

c3

cl

c2

c2

Cl

-b0

c2

c3

C1

The procedure for calculating the roots therefore is as follows: with initial estimates of p and q (zero or one can be used), calculate the values of b] and cs. Use these values to calculate Ap and Aq, and correct the initial values. Continue until convergence is reached. Obtain the two roots by use of the quadratic formula. Use the result of synthetic division of the polynomial as the new polynomial, and repeat the process. Continue until the polynomial is of order one or zero.

The VBA code is shown in Figure 8-28. The portion of the code that performs the Bairstow calculation is based on code found in Shoup, T. E., Numerical Methods for the Personal Computer, Prentice-Hall, 1983.

This procedure contains code, not found in other procedures in this book, that allows the macro to accept a polynomial equation as a reference to a cell that contains a formula or as a reference to a cell that contains a formula as text. The procedure also handles an implicit reference.

Option Explicit

Function Bairstow(equation, reference)

'Obtains the coefficients of a regular polynomial (maximum order = 6).

'Polynomial is a cell formula.

'Polynomial can contain cell references or names.

'Poynomial can be text.

'Reference can be a cell reference or a name.

Dim J As Integer, N As Integer

Dim p1 As Integer, p2 As Integer, p3 As Integer

Dim expnumber As Integer, ParenFlag As Integer

Dim R As Integer, C As Integer

Dim FormulaText As String, RefText As String, NameText As String Dim char As String, term As String ReDim A(6)

' GET equation EITHER AS CELL FORMULA OR AS TEXT. If Application.IsText(equation) Then

FormulaText = equation 'If in quotes, remove them. If Asc(Left(FormulaText, 1)) = 34 Then _ FormulaText = Mid(FormulaText, 2, Len(FormulaText) -1) Else

FormulaText = equation.Formula End If

If Left(FormulaText, 1) = "=" Then FormulaText = Mid(FormulaText, 2, 1024) FormulaText = Application.ConvertFormula(FormulaText, xlA1, xlA1, _ xIAbsolute)

FormulaText = Application.Substitute(FormulaText,"","") 'remove all spaces

'GET THE NAME CORRESPONDING TO reference NameText =""

On Error Resume Next 'Handles case where no name has been assigned NameText = reference.Name.Name On Error GoTo 0

NameText = Mid(NameText, lnStr(1, NameText,"!") + 1)

'HANDLE CASE WHERE reference IS A RANGE 'by finding cell in same row or column as cell containing function. If reference.Rows.Count > 1 Then R = equation.Row

Set reference = lntersect(reference, Range(R & ":" & R)) Elself reference.Columns.Count > 1 Then C = equation.Column

Set reference = lntersect(reference, Range(C & ":" & C)) _

End If

RefText = reference.Address

'PARSE THE FORMULA INTO TERMS 'pointers: p1, beginning; p2, end of string.

FormulaText = FormulaText & add extra character for parsing

ParenFlag = 0 'Keep track of left and right parentheses For J = 1 To Len(FormulaText) char = Mid(FormulaText, J, 1) If char = "(" Then ParenFlag = ParenFlag + 1 If char = ")" Then ParenFlag = ParenFlag -1

If ((char = "+" Or char = "-") And ParenFlag = 0) Or J = Len(FormulaText) Then term = Mid(FormulaText, p1, J - p1)

term = Application.Substitute(term, NameText, RefText)

•GET THE EXPONENT AND COEFFICIENT FOR EACH TERM 'p3: location of reference in term.

If lnStr(1, term, RefText & "A") Then 'function returns zero if not found These are the xA2 and higher terms p3 = lnStr(1, term, RefText & "A") expnumber = Mid(term, p3 + Len(RefText) + 1,1) term = Left(term, p3 -1) 'term is now the coefficient part Elself lnStr(1, term, RefText) Then 'This is the x term p3 = lnStr(1, term, RefT ext) expnumber = 1

term = Left(term, p3 -1) 'term is now the coefficient part Else

'This is the constant term expnumber = 0 End If

If term = "" Then term = "=1 " 'If missing, Evaluate will require a string.

If term = "+" Or term = "-" Then term = term & "1"

If Right(term, 1) = "*" Then term = Left(term, Len(term) -1)

A(expnumber) = Evaluate(term)

End If

Next J

'RESIZE THE ARRAY

Next

ReDim Preserve A(N) ReDim Root(1 To N, 1)

•REDUCE POLYNOMIAL SO THAT FIRST COEFF = 1 For J = 0 To N: A(J) = A(J) / A(N): Next

'SCALE THE POLYNOMIAL, IF NECESSARY '<code to be added later>

Call EvaluateByBairstowMethod(N, A, Root) Bairstow = Root()

End Function

Sub EvaluateByBairstowMethod(N, A, Root)

'Code adapted from Shoup, "Numerical Methods for the Personal Computer".

Dim M As Integer, I As Integer, J As Integer, IT As Integer

Dim P As Double, Q As Double, delP As Double, delQ As Double

Dim denom As Double, S1 As Double

Dim tolerance As Double

ReDim B(N + 2), C(N + 2) tolerance = 0.000000000000001 M = N

If M = 1 Then Root(M, 0) = -A(0): Call Sort(Root, N): Exit Sub

If Abs(delP) < tolerance And Abs(delQ) < tolerance Then Exit For For J = 0 To M B(M - J) = A(M - J) + P * B(M - J + 1) + Q * B(M - J + 2) C(M - J) = B(M - J) + P * C(M - J + 1) + Q * C(M - J + 2) Next J

denom = C(2) A 2 - C(1) * C(3) delP = (-B(1) * C(2) + B(0) * C(3)) / denom delQ = (-C(2) * B(0) + C(1) * B(1)) / denom P = P + delP Q = Q + delQ Next IT

S1 = PA2 + 4*Q If S1 < 0 Then 'Handle imaginary roots Root(M, 0) = P/2: Root(M, 1) = Sqr(-S1) 12 Root(M -1, 0) = P / 2: Root(M -1,1) = -Sqr(-S1) / 2 Else

'Handle real roots Root(M, 0) = (P + Sqr(S1)) / 2 Root(M -1, 0) = (P - Sqr(S 1)) / 2 End If

For I = M To 0 Step -1: A(l) = B(l + 2): Next M = M - 2 Wend End Sub i. __, , , , , , , I i , . , . . . I j. 4.4.4.4.4.^4..^^+ 4.

Sub Sort(Root, N)

'SORT ROOTS IN ASCENDING ORDER

Dim I As Integer, J As Integer_

Dim tempO As Double, tempi As Double

If Root(l, 0) > Root(J, 0) Then tempO = Root(l, 0): tempi = Root(l, 1) Root(l, 0) = Root(J, 0): Root(l, 1) = Root(J, 1) Root(J, 0) = tempO: Root(J, 1) = tempi End If Next J Next I End Sub

Figure 8-28. VBA code for the Bairstow custom function, (folder 'Chapter 08 Examples', workbook 'Bairstow', module 'BairstowFn')

The syntax of the Bairstow function is

Bairstow(eguai/'on,reference)

Equation is a reference to a cell that contains the formula of the function, reference is the cell reference of the argument to be varied (the jc value of F(x)).

The Bairstow function is an array function. To return the roots of a polynomial of order N, you must select a range of cells 2 columns by N rows, enter the function and then press CONTROL+SHIFT+ENTER.

Figure 8-29 shows a chart of the polynomial y = x3 - 0.003 lx2 + 2.3 x 10"sx + 5 x 10-9

Figure 8-29. A regular polynomial with one real root and two imaginary roots, (folder 'Chapter 08 Examples', workbook 'Bairstow', sheet 'Example')

The function has one real root and a pair of imaginary roots. Figure 8-30 shows a portion of the spreadsheet in which the Bairstow custom function is used to obtain the roots of the function.

A

B I

25

Roots:

26

Real part

imaginary part

27

-0.0D109Q

o

28

0.002095

-0.000447311

29

0.002095

0.000447311

Figure 8-30. Calculation of all roots (real and imaginary) of a regular polynomial by the Bairstow custom function, (folder 'Chapter 08 Examples', workbook 'Bairstow', sheet 'Example 2')

Figure 8-30. Calculation of all roots (real and imaginary) of a regular polynomial by the Bairstow custom function, (folder 'Chapter 08 Examples', workbook 'Bairstow', sheet 'Example 2')

The formula

=A2A3-0.0031 *A2A2+0.000000023*A2+0.000000005 was entered in cell B2 and the Bairstow custom function {=Bairstow(B2,A2)}

in cells A27:B29. The real part of the root is in the left cell and the imaginary part in the right cell. Note that, since the custom function handles only polynomials with real coefficients, the complex roots (if any) occur in conjugate pairs.

+1 -1

Responses

Post a comment