(15C) Bairstow's Method
02-25-2022, 01:20 AM
Post: #1
 Thomas Klemm Senior Member Posts: 1,814 Joined: Dec 2013
(15C) Bairstow's Method
Bairstow's Method

Links

Wikipedia
MathWorld

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 R2 - R7 with the corresponding values of the equation:

$$\begin{matrix} n & \rightarrow & R_2 \\ \sum x & \rightarrow & R_3 \\ \sum x^2 & \rightarrow & R_4 \\ \sum y & \rightarrow & R_5 \\ \sum y^2 & \rightarrow & R_6 \\ \sum xy & \rightarrow & R_7 \end{matrix}$$

This mapping leaves us with the following linear equation:

$$\begin{bmatrix} R_4 & R_3 \\ R_3 & R_2 \end{bmatrix}\cdot\begin{bmatrix} Y \\ X \end{bmatrix}=\begin{bmatrix} R_7 \\ R_5 \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 4
-3 STO 3
5 STO 2
13 STO 7
-21 STO 5
L.R.

Y: 2
X: -3

Registers

$$\begin{matrix} R_0 & p \\ R_1 & q \\ R_2 & c = c_j \\ R_3 & {c}′ = c_{j+1} \\ R_4 & {c}'' = c_{j+2} \\ R_5 & b = b_i \\ R_6 & \\ R_7 & {b}′ = b_{i+1} \\ R_8 & \text{index} = 9.fff \\ 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 } f LBL A 002 {       42 32 } f ∑ 003 {       45  8 } RCL 8 004 {       44 25 } STO I 005 {    42 21  0 } f LBL 0 006 {       45  5 } RCL 5 007 {       45  3 } RCL 3 008 {       44  4 } STO 4 009 {    45 20  1 } RCL × 1 010 {          30 } − 011 {       45  2 } RCL 2 012 {       44  3 } STO 3 013 {    45 20  0 } RCL × 0 014 {          30 } − 015 {       44  2 } STO 2 016 {       45 24 } RCL (i) 017 {       45  7 } RCL 7 018 {    45 20  1 } RCL × 1 019 {          30 } − 020 {       45  5 } RCL 5 021 {       44  7 } STO 7 022 {    45 20  0 } RCL × 0 023 {          30 } − 024 {       44  5 } STO 5 025 {    42  6 25 } f ISG I 026 {       22  0 } GTO 0 027 {       42 49 } f L.R. 028 {    44 40  0 } STO + 0 029 {          34 } x↔y 030 {    44 40  1 } STO + 1 031 {       43  1 } g →P 032 {       43 34 } g RND 033 {    43 30  0 } g TEST x≠0 034 {       22 11 } GTO A 035 {       45  8 } RCL 8 036 {           2 } 2 037 {          26 } EEX 038 {           3 } 3 039 {          16 } CHS 040 {          30 } − 041 {       44  8 } STO 8 042 {       44 25 } STO I 043 {           0 } 0 044 {          36 } ENTER 045 {    42 21  1 } f LBL 1 046 {       45 24 } RCL (i) 047 {       45  1 } RCL 1 048 {       43 33 } g R⬆ 049 {          20 } × 050 {          30 } − 051 {       45  0 } RCL 0 052 {       43 33 } g R⬆ 053 {          20 } × 054 {          30 } − 055 {       44 24 } STO (i) 056 {    42  6 25 } f ISG I 057 {       22  1 } GTO 1 058 {       43 32 } g RTN

Description

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

Partial Derivatives
Lines 005 - 026
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 - 015
$$(c, {c}′, {c}'') \leftarrow (b - c p - {c}′ q, c, {c}′)$$

Lines 016 - 024
$$(b, {b}′) \leftarrow (a_k - b p - {b}′ q, b)$$

Solve the linear equation
Line 027
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 028 - 030
$$p \leftarrow p + \delta p$$
$$q \leftarrow q + \delta q$$

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

Polynomial Division
Lines 035 - 057
We have to repeat the division here since we didn't keep $$b_k$$ in lines 016 - 024.

Quadratic Solver

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

Code:
059 {    42 21 12 } f LBL B 060 {       45  0 } RCL 0 061 {           2 } 2 062 {          16 } CHS 063 {          10 } ÷ 064 {          36 } ENTER 065 {          36 } ENTER 066 {       43 11 } g x² 067 {    45 30  1 } RCL − 1 068 {          11 } √x̅ 069 {          30 } − 070 {          34 } x↔y 071 {       43 36 } g LSTx 072 {          40 } + 073 {       43 32 } g 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 8 1      STO 0        STO 1

Alternatively use:
Code:
MATRIX 1

Run program

Code:
GSB A        -52.0000 RCL 0          1.5000 RCL 1         -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 0        STO 1

Alternatively use:
Code:
MATRIX 1

Run program again

Code:
GSB A         -4.0000 RCL 0         -4.0000 RCL 1         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-2)(x+3)(x^2-4x+13)$$

Solutions

$$x_1=1.5$$
$$x_2=2$$
$$x_3=-3$$
$$x_4=2+3i$$
$$x_5=2-3i$$

Attached File(s)
bairstow-15c.zip (Size: 1.95 KB / Downloads: 3)
 « Next Oldest | Next Newest »

 Messages In This Thread (15C) Bairstow's Method - Thomas Klemm - 02-25-2022 01:20 AM RE: (15C) Bairstow's Method - rprosperi - 02-25-2022, 01:58 AM RE: (15C) Bairstow's Method - Thomas Klemm - 02-25-2022, 07:04 AM RE: (15C) Bairstow's Method - Albert Chan - 02-27-2022, 12:40 PM

User(s) browsing this thread: 1 Guest(s)