(11C) Bairstow's Method

[Thomas Klemm]

Bairstow's Method

Links

Wikipedia
MathWorld
Numerical Recipes in Fortran 77
9.5 Roots of Polynomials

Algorithm

We are looking for a quadratic factor T(x) of the polynomial P(x). In general there will be a linear remainder R(x) after the division:

\(P(x) = Q(x) \cdot T(x) + R(x)\)

\(P(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0\)

\(T(x) = x^2 + p x + q\)

\(Q(x) = b_n x^{n-2} + b_{n-1} x^{n-3} + \cdots + b_4 x^2 + b_3 x + b_2\)

\(R(x) = b_1 ( x + p ) + b_0\)

If we set \(b_{n+1}\equiv b_{n+2}\equiv 0\), then \(\forall k\leq n\):

\(b_k = a_k - b_{k+1} p - b_{k+2} q\)

The two coefficients of the remainder R(x) depend on p and q:

\(b_0 = b_0(p, q)\)

\(b_1 = b_1(p, q)\)

To make T(x) a factor of P(x) the remainder R(x) must vanish. Let's assume that we already have p and q close to the solution. We want to know how to change these values to get a better solution. Therefore we set the first order Taylor series expansion to 0:

\(b_0(p + \delta p, q + \delta q) \approx b_0(p, q) + \frac{\partial b_0}{\partial p} \delta p + \frac{\partial b_0}{\partial q} \delta q = 0\)

\(b_1(p + \delta p, q + \delta q) \approx b_1(p, q) + \frac{\partial b_1}{\partial p} \delta p + \frac{\partial b_1}{\partial q} \delta q = 0 \)

It turns out that the partial derivatives of \(b_k\) with respect to p and q can be calculated with a similar recurrence formula:

\(c_{k+1} = - \frac{\partial b_k}{\partial p}\)

\(c_{k+2} = - \frac{\partial b_k}{\partial q}\)

\(c_k = b_k - c_{k+1} p - c_{k+2} q\)

Thus we have to solve:

\(\begin{matrix}
c_1 \delta p + c_2 \delta q &= b_0 \\
c_2 \delta p + c_3 \delta q &= b_1
\end{matrix}\)


Using Linear Regression

We have to solve a 2-dimensional linear equation:

\(\begin{bmatrix}
c_1 & c_2 \\
c_2 & c_3
\end{bmatrix}\cdot\begin{bmatrix}
\delta p \\
\delta q
\end{bmatrix}=\begin{bmatrix}
b_0 \\
b_1
\end{bmatrix}\)

But since the matrix is symmetric we can use the built-in function L.R. for linear regression.

We know that the following equation is solved for A and B:

\(\begin{bmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{bmatrix}\cdot\begin{bmatrix}
A \\
B
\end{bmatrix}=\begin{bmatrix}
\sum xy \\
\sum y
\end{bmatrix}\)

Cramer's Rule gives the solutions:

\(A=\frac{\begin{vmatrix}
\sum xy & \sum x \\
\sum y & n
\end{vmatrix}}{\begin{vmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{vmatrix}}=\frac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x)^2}\)

\(B=\frac{\begin{vmatrix}
\sum x^2 & \sum xy \\
\sum x & \sum y
\end{vmatrix}}{\begin{vmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{vmatrix}}=\frac{\sum y\sum x^2-\sum x\sum xy}{n\sum x^2-(\sum x)^2}\)

For this to work we have to fill the registers R₀ - R₅ with the corresponding values of the equation:

\(\begin{matrix}
n & \rightarrow & R_0 \\
\sum x & \rightarrow & R_1 \\
\sum x^2 & \rightarrow & R_2 \\
\sum y & \rightarrow & R_3 \\
\sum y^2 & \rightarrow & R_4 \\
\sum xy & \rightarrow & R_5
\end{matrix}\)

This mapping leaves us with the following linear equation:

\(\begin{bmatrix}
R_2 & R_1 \\
R_1 & R_0
\end{bmatrix}\cdot\begin{bmatrix}
Y \\
X
\end{bmatrix}=\begin{bmatrix}
R_5 \\
R_3
\end{bmatrix}\)

Example:

\(\begin{bmatrix}
2 & -3 \\
-3 & 5
\end{bmatrix}\cdot\begin{bmatrix}
2 \\
-3
\end{bmatrix}=\begin{bmatrix}
13 \\
-21
\end{bmatrix}\)

CLEAR Σ
2 STO 2
-3 STO 1
5 STO 0
13 STO 5
-21 STO 3
L.R.

Y: 2
X: -3


Registers

\(\begin{matrix}
R_0 & c = c_j \\
R_1 & {c}′ = c_{j+1} \\
R_2 & {c}'' = c_{j+2} \\
R_3 & b = b_i \\
R_4 & \\
R_5 & {b}′ = b_{i+1} \\
R_6 & \text{index} = 9.fff \\
R_7 & p \\
R_8 & q \\
R_9 & a_n \\
R_{.0} & a_{n-1} \\
R_{.1} & a_{n-2} \\
R_{.2} & \cdots \\
R_{.3} & \cdots \\
\end{matrix}\)

Program

Code:

001 - 42,21,11   LBL A                             032 - 44,40, 7   STO + 7
002 -   42  32   CLEAR Σ                           033 -       34   x<>y
003 -    45  6   RCL 6                             034 - 44,40, 8   STO + 8
004 -   44  25   STO I                             035 -   43  26   →P
005 - 42,21, 0   LBL 0                             036 -   43  34   RND
006 -    45  3   RCL 3                             037 -   43  30   x≠0
007 -    45  1   RCL 1                             038 -   22  11   GTO A
008 -    44  2   STO 2                             039 -    45  6   RCL 6
009 -    45  8   RCL 8                             040 -        2   2
010 -       20   ×                                 041 -       26   EEX
011 -       30   -                                 042 -        3   3
012 -    45  0   RCL 0                             043 -       16   CHS
013 -    44  1   STO 1                             044 -       30   -
014 -    45  7   RCL 7                             045 -    44  6   STO 6
015 -       20   ×                                 046 -   44  25   STO I
016 -       30   -                                 047 -        0   0
017 -    44  0   STO 0                             048 -       36   ENTER
018 -   45  24   RCL (i)                           049 - 42,21, 1   LBL 1
019 -    45  5   RCL 5                             050 -   45  24   RCL (i)
020 -    45  8   RCL 8                             051 -    45  8   RCL 8
021 -       20   ×                                 052 -   43  33   R↑
022 -       30   -                                 053 -       20   ×
023 -    45  3   RCL 3                             054 -       30   -
024 -    44  5   STO 5                             055 -    45  7   RCL 7
025 -    45  7   RCL 7                             056 -   43  33   R↑
026 -       20   ×                                 057 -       20   ×
027 -       30   -                                 058 -       30   -
028 -    44  3   STO 3                             059 -   44  24   STO (i)
029 -    42  6   ISG                               060 -    42  6   ISG
030 -    22  0   GTO 0                             061 -    22  1   GTO 1
031 -   42  49   L.R.                              062 -   43  32   RTN

Description

Initialization k = n
Lines 002 - 004
\((b, {b}′, c, {c}′, {c}'') \leftarrow (0, 0, 0, 0, 0)\)
\(I \leftarrow\) index

Polynomial division
Lines 005 - 030
At the end of the loop \(b = b_0\) but \(c = c_1\)!
That's why the second division is performed before the first.

Iteration k+1 → k
Lines 006 - 017
\((c, {c}′, {c}'') \leftarrow (b - c p - {c}′ q, c, {c}′)\)

Lines 018 - 028
\((b, {b}′) \leftarrow (a_k - b p - {b}′ q, b)\)

Solve the linear equation
Line 031
The trick using linear regression only works since the matrix is symmetric.
\(\begin{bmatrix}
c & {c}' \\
{c}' & {c}''
\end{bmatrix}\cdot\begin{bmatrix}
\delta p \\
\delta q
\end{bmatrix}=\begin{bmatrix}
b \\
{b}'
\end{bmatrix}\)

Find improved values for p and q
Lines 032 - 034
\(p \leftarrow p + \delta p\)
\(q \leftarrow q + \delta q\)

Stop Criterion
Lines 035 - 037
\(\left | (\delta p, \delta q) \right |< \varepsilon\)
\(\varepsilon = 10^{-n}\) if display is set to FIX n

Polynomial Division
Lines 039 - 061
We have to repeat the division here since we didn't keep \(b_k\) in lines 018 - 028.


Quadratic Solver

This program solves the quadratic equation: \(T(x)=x^2+px+q=0\)

Code:

063 - 42,21,12  LBL B
064 -    45  7  RCL 7
065 -        2  2
066 -       16  CHS
067 -       10  ÷
068 -       36  ENTER
069 -       36  ENTER
070 -    43 11  x^2
071 -    45  8  RCL 8
072 -       30  -
073 -       11  SQRT
074 -       30  -
075 -       34  x<>y
076 -    43 36  LSTx
077 -       40  +
078 -    43 32  RTN

Just be aware that this program can't find complex roots. Instead an Error 0 will be displayed.
However it's easy to find the complex solutions. Just use:

CHS
\(\sqrt{x}\)

The solutions then are: \(Y \pm iX\)


Example

\(P(x)=2x^5-9x^4+15x^3+65x^2-267x+234=0\)

Insert coefficients

Code:

2      STO 9
-9     STO .0
15     STO .1
65     STO .2
-267   STO .3
234    STO .4

Initialization

Code:

9.014  STO 6
1      STO 7
       STO 8

Run program

Code:

GSB A        -52.0000

RCL 7          1.5000
RCL 8         -4.5000

RCL 9          2.0000
RCL .0       -12.0000
RCL .1        42.0000
RCL .2       -52.0000

Code:

GSB B          1.5000
x<>y          -3.0000

Conclusion

\(2x^5-9x^4+15x^3+65x^2-267x+234=\)
\((x^2+1.5x-4.5)(2x^3-12x^2+42x-52)\)

Solutions

For \(x^2+1.5x-4.5=0\):
\(x_1=1.5\)
\(x_2=-3\)

Initialize guess

Code:

1      STO 7
       STO 8

Run program again

Code:

GSB A         -4.0000

RCL 7         -4.0000
RCL 8         13.0000

RCL 9          2.0000
RCL .0        -4.0000

Code:

GSB B          Error 0
←             -9.0000
CHS            9.0000
√x             3.0000
x<>y           2.0000

Code:

RCL .0        -4.0000
RCL 9          2.0000
÷             -2.0000
CHS            2.0000

Conclusion

\(2x^3-12x^2+42x-52=\)
\((x^2-4x+13)(2x-4)\)

Solutions

For \(x^2-4x+13=0\):
\(x_3=2+3i\)
\(x_4=2-3i\)

For \(2x-4=0\):
\(x_5=2\)

Summary

Factors

\(2x^5-9x^4+15x^3+65x^2-267x+234=\)
\((x^2+1.5x-4.5)(x^2-4x+13)(2x-4)=\)
\((x-1.5)(x+3)(x^2-4x+13)2(x-2)=\)
\((2x-3)(x+3)(x^2-4x+13)(x-2)\)

Solutions

\(x_1=1.5\)
\(x_2=-3\)
\(x_3=2+3i\)
\(x_4=2-3i\)
\(x_5=2\)