Weierstrass-Durand-Kerner Method for the HP-41

Overview
 

 1°)  Polynomials with Real Coefficients & All Roots Real
 2°)  Complex Polynomials ( 2 programs )

-This iterative method finds simultaneously all the roots of a polynomial !

-To solve  p(x) = xⁿ + a_n-1 x^n-1 + ..... + a₁ x + a₀ = 0   it takes initial approximations  r₁ , ...... , r_n

-The next approximations are:   r'_k = r_k - p(r_k) / ¶_j#k( r_k - r_j )                                                                 ¶ = product
-The iterations stop when the difference between successive approximations are enough small.
-In the following programs, this is controlled by the display format.

1°)  Polynomials with Real Coefficients & All Roots Real
 

Data Registers:              R00 to R04: temp                                     ( Registers Rbb thru Ree+nn are to be initialized before executing "WDK" )

                                      •  Rbb = a_n-1   •  Rbb+1 = a_n-2  ...........   •  Ree-1 = a₁   •  Ree = a₀     coefficients of the polynomial        ( bb > 04 )

                                      •  Ree+1 = r₁   •  Ree+2 = r₂  .............    •  Ree+nn = r_n                        initial approximations of the roots ( not too closed )
Flags: /
Subroutines: /

01 LBL "WDK"  02 STO 01  03 STO 02          04 INT  05 CHS  06 LASTX  07 FRC  08  E3  09 *  10 +  11 1  12 +  13 STO 00

14 .1  15 %  16 +  17 ST+ 01  18 CLX  19 STO 03          20 RCL 01  21 STO 04  22 LBL 01  23 RCL IND 01  24 ENTER  25 ENTER  26 ENTER

27 CLX  28 SIGN  29 LBL 02          30 *  31 RCL IND 02  32 +  33 ISG 02  34 GTO 02  35 LBL 03  36 X<>Y  37 RCL IND 04  38 -  39 X=0?

40 SIGN  41 /  42 ISG 04  43 GTO 03  44 ST- IND 01  45 ABS  46 ST+ 03  47 RCL 00          48 ST- 02  49 ST- 04  50 ISG 01  51 GTO 01  52 ST- 01

53 CLX  54 X<> 03  55 RCL 00  56 /  57 RND  58 VIEW X  59 X#0?  60 GTO 01  61 RCL 02          62 X<> L  63 RCL 01  64 CLD  65 END

 
      ( 94 bytes / SIZE var. )

 

STACK	INPUTS	OUTPUTS
Y	/	epsilon
X	bbb.eee	1+eee.eee+nnn
L	/	bbb.eee

 
Example:      Find all the roots of  p(x) = 2 x⁵ + 3 x⁴ - 35 x³ - 10 x² + 128 x - 74

-First, divide all the coefficients by a₅ = 2 , it yields     x⁵ + 1.5 x⁴ - 17.5 x³ - 5 x² + 64 x - 37

-Then, store these coefficients in - for instance - R11 to R15:  

       1.5   STO 11
    -17.5  STO 12
      -5     STO 13      and the initial approximations of the roots ( say:  -4  -2  0  2  4  )  into   R16-R17-R18-R19-R20   ( just after the coefficients )
      64    STO 14
     -37    STO 15       -4  STO 16   -2  STO 17   0  STO 18  2  STO 19  4  STO 20

-Let's try the program in FIX 8:

     11.015   XEQ "WDK"   >>>>    16.020                                ---Execution time = 74s---
                                           X<>Y   5.7 E-10


-And we have the 5 roots in R16 to R20

    R16 = -4.373739464
    R17 = -2.455070117
    R18 =  0.703611644
    R19 =  1.641131729
    R20 =  2.984066207

Notes:

-The HP41 displays the successive sums of the | differences | between 2 consecutive approximations
-An infinite loop may happen in FIX 9:  FIX 8 is often enough to get the best precision.

-This routine doesn't work if all the roots are not real.
-In this case, employ the program below:


2°)  Complex Polynomials


-Now, we solve:  p(z) =  zⁿ + a_n-1 z^n-1 + ..... + a₁ z + a₀ = 0

  where   a_k = b_k + i c_k

 from n initial approximations   r_k = s_k + i t_k

Data Registers:              R00 to R10: temp                                    ( Registers Rbb thru Ree+2n are to be initialized before executing "WDKZ1" )

                                      •  Rbb = b_n-1   •  Rbb+1 = c_n-1  .........    •  Ree-1 = b₀          •  Ree = c₀                         coefficients of the polynomial        ( bb > 10 )

                                      •  Ree+1 = s₁   •  Ree+2 = t₁  .............    •  Ree+2n-1 = s_n    •  Ree+2n = t_n                    initial approximations of the roots ( not too closed )
Flags: /
Subroutines: /


-Lines 125 & 134 are three-byte GTOs

01 LBL "WDKZ1"  02 STO 01  03 STO 02            04 INT  05 CHS  06 LASTX  07 FRC  08  E3  09 *  10 +  11 1  12 +  13 STO 00  14 .1  15 %  16 +  17 ST+ 02  18 CLX  19 STO 03  20 RCL 02

21 STO 04  22 LBL 01  23 RCL IND 02  24 STO 05            25 ISG 02  26 RCL IND 02  27 STO 06  28 CLX  29 STO 08  30 STO 10  31 SIGN  32 STO 07  33 STO 09  34 LBL 02  35 RCL 05  36 RCL 07  37 *  38 RCL 06  39 RCL 08  40 *

41 -  42 RCL IND 01  43 +  44 X<> 07  45 RCL 06            46 *  47 RCL 05  48 RCL 08  49 *  50 +  51 ISG 01  52 RCL IND 01  53 +  54 STO 08  55 ISG 01  56 GTO 02  57 LBL 03  58 RCL 05  59 RCL IND 04  60 -

61 ISG 04  62 RCL 06  63 RCL IND 04  64 -  65 X=Y?  66 X#0?  67 FS? 30  68 GTO 03  69 STO Z  70 RCL 09            71 *  72 RCL 10  73 X<> Z  74 ST* Z  75 RDN  76 +  77 X<> 10  78 *  79 X<>Y  80 RCL 09

81 *  82 X<>Y  83 -  84 STO 09  85 LBL 03  86 ISG 04  87 GTO 03  88 RCL 08  89 RCL 09            90 *  91 RCL 07  92 RCL 10  93 *  94 -  95 RCL 09  96 X^2  97 RCL 10  98 X^2  99 + 100 STO 06

101 / 102 ST- IND 02 103 ABS 104 ST+ 03 105 RCL 02           106 ENTER 107 SIGN 108 - 109 RCL 07 110 RCL 09 111 * 112 RCL 08 113 RCL 10 114 * 115 + 116 RCL 06 117 / 118 ST- IND Y 119 ABS 120 ST+ 03

121 RCL 00           122 ST- 01 123 ST- 04 124 ISG 02 125 GTO 01 126 ST- 02 127 CLX 128 X<> 03 129 RCL 00 130 / 131 RND 132 VIEW X 133 X#0? 134 GTO 01 135 RCL 01 136 X<> L 137 RCL 02 138 CLD 139 END

       
      ( 187 bytes / SIZE var.)

 

STACK	INPUTS	OUTPUTS
Y	/	epsilon
X	bbb.eee	1+eee.eee+2n
L	/	bbb.eee

 
Example:      Find all the roots of  p(z) = (1+2.i) z⁵ + (3-4.i) z⁴ + (5-i) z³ + (2+3.i) z² + (1-3.i) z + (3+4.i)

-We must divide all the coefficients by 1+2.i , it yields:

    z⁵ + (-1-2.i) z⁴ + (0.6-2.2 i) z³ + (1.6-0.2 i) z² + (-1-i) z + (2.2-0.4 i)
 
-For instance, let's store the coefficients into R11 R12 ...

   -1  STO 11  -2   STO 12
  0.6 STO 13 -2.2  STO 14           and if we choose  1+i   1-i   -1+i   -1-i    0   as initial approximations:
  1.6 STO 15 -0.2  STO 16
   -1  STO 17  -1   STO 18           1   STO 21  STO 22  STO 23   -1  STO 24  STO 25   1  STO 26   -1  STO 27  STO 28  0  STO 29  STO 30
   2.2 STO 19 -0.4 STO 20

-In FIX 8  

    11.020   XEQ "WDKZ1"   >>>>    21.030                                       ---Execution time = 7m22s---
                                              X<>Y   7.7 E-11

-And we find the 5 complex roots in registers R21 to R30:

    1.588278774 + 2.620754096 i  = R21 + i R22
    0.581317453 -  0.537568691 i  = R23 + i R24
    0.220625249 + 0.812343651 i  = R25 + i R26
   -1.209218706 - 0.008454644 i  = R27 + i R28
   -0.181002770 - 0.887074412 i  = R29 + i R30


Notes:

-If the coefficients are real, it's not difficult to divide them by a_n
-But it's less easy when a_n is complex

-In the following variant, we solve directly   p(z) = a_n zⁿ + a_n-1 z^n-1 + ..... + a₁ z + a₀ = 0

  where   a_k = b_k + i c_k

 from n initial approximations   r_k = s_k + i t_k

Data Registers:              R00 to R12: temp                                     ( Registers Rbb thru Ree+2n are to be initialized before executing "WDKZ" )

                                      •  Rbb = b_n      •  Rbb+1 = c_n  ...........    •  Ree-1 = b₀          •  Ree = c₀                         coefficients of the polynomial        ( bb > 12 )

                                      •  Ree+1 = s₁   •  Ree+2 = t₁  .............    •  Ree+2n-1 = s_n    •  Ree+2n = t_n                    initial approximations of the roots ( not too closed )
Flags: /
Subroutines: /


-Lines 149 & 158 are three-byte GTOs

01 LBL "WDKZ"  02 STO 01  03 RCL IND 01  04 STO 11  05 ISG 01  06 RCL IND 01  07 STO 12  08 ISG 01  09 RCL 01  10 STO 02          11 INT  12 CHS  13 LASTX  14 FRC  15  E3  16 *  17 +  18 1  19 +  20 STO 00  21 .1  22 %  23 +

24 ST+ 02  25 CLX  26 STO 03  27 RCL 02  28 STO 04  29 LBL 01  30 RCL IND 02  31 STO 05  32 ISG 02  33 RCL IND 02  34 STO 06          35 CLX  36 STO 10  37 SIGN  38 STO 09  39 RCL 11  40 STO 07  41 RCL 12  42 STO 08  43 LBL 02  44 RCL 05  45 RCL 07  46 *

47 RCL 06  48 RCL 08  49 *  50 -  51 RCL IND 01  52 +  53 X<> 07  54 RCL 06  55 *  56 RCL 05  57 RCL 08          58 *  59 +  60 ISG 01  61 RCL IND 01  62 +  63 STO 08  64 ISG 01  65 GTO 02  66 LBL 03  67 RCL 05  68 RCL IND 04  69 -

70 ISG 04  71 RCL 06  72 RCL IND 04  73 -  74 X=Y?  75 X#0?  76 FS? 30  77 GTO 03  78 STO Z  79 RCL 09          80 *  81 RCL 10  82 X<> Z  83 ST* Z  84 RDN  85 +  86 X<> 10  87 *  88 X<>Y  89 RCL 09  90 *  91 X<>Y  92 -

93 STO 09  94 LBL 03  95 ISG 04  96 GTO 03  97 RCL 09  98 RCL 11  99 * 100 RCL 10 101 RCL 12         102 * 103 - 104 X<> 09 105 RCL 12 106 * 107 RCL 10 108 RCL 11 109 * 110 + 111 STO 10 112 RCL 08 113 RCL 09 114 * 115 X<>Y

116 RCL 07 117 * 118 - 119 RCL 09 120 X^2 121 RCL 10 122 X^2 123 + 124 STO 06         125 / 126 ST- IND 02 127 ABS 128 ST+ 03 129 RCL 02 130 ENTER 131 SIGN 132 - 133 RCL 07 134 RCL 09 135 * 136 RCL 08 137 RCL 10 138 *

139 + 140 RCL 06 141 / 142 ST- IND Y 143 ABS 144 ST+ 03 145 RCL 00         146 ST- 01 147 ST- 04 148 ISG 02 149 GTO 01 150 ST- 02 151 CLX 152 X<> 03 153 RCL 00 154 / 155 RND 156 VIEW X 157 X#0? 158 GTO 01 159 LASTX 160 RCL 02 161 CLD 162 END

   
      ( 213 bytes / SIZE var.)

 

STACK	INPUTS	OUTPUTS
Y	/	epsilon
X	bbb.eee	1+eee.eee+2n
L	/	bbb.eee

 
Example:      Find all the roots of  p(z) = (1+2.i) z⁵ + (3-4.i) z⁴ + (5-i) z³ + (2+3.i) z² + (1-3.i) z + (3+4.i)
 
-For instance, let's store the coefficients into R21 R22 ...

    1  STO 21   2  STO 22
    3  STO 23  -4  STO 24     and if we choose  1+i   1-i   -1+i   -1-i    0   as initial approximations:
    5  STO 25  -1  STO 26
    2  STO 27   3  STO 28      1   STO 33  STO 34  STO 35   -1  STO 36  STO 37   1  STO 38   -1  STO 39  STO 40  0  STO 41  STO 42
    1  STO 29  -3  STO 30
    3  STO 31   4  STO 32

-In FIX 8  

    21.032   XEQ "WDKZ"   >>>>    33.042                                       ---Execution time = 7m43s---
                                            X<>Y   7.4 E-11

-And we find in registers R33 to R42:

    1.588278774 + 2.620754096 i  = R33 + i R34
    0.581317453 -  0.537568691 i  = R35 + i R36
    0.220625249 + 0.812343651 i  = R37 + i R38
   -1.209218706 - 0.008454644 i  = R39 + i R40
   -0.181002770 - 0.887074412 i  = R41 + i R42

Reference:

[1]  https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method