POLYROOTS Function for Excel

Work in Progress :: Please contact us to report errors, typos, etc.

=L_POLYROOTS(coeffs, [iterations], [epsilon], [iterations], [xstart], [return_abs])

Argument	Description	Example
coeffs	Coefficients of the polynomial as a row vector in decreasing order of power	{1,7,0,4}
iterations	(default=42) The number of iterations	42
epsilon	(if blank, runs all iterations) Desired value of epsilon where \|p(x)\|<ε	0.000000000001
xstart	(default="0.4+0.9i") A starting value of x for the algorithm	"0.4+0.9i"
return_abs	(default=FALSE) If TRUE, returns IMABS(p(x)) for checking convergence	FALSE

Requires: L_CIRCSHIFT, L_POLYVAL

Description

L_POLYROOTS is the first function in the Lambda library to be able to calculate complex roots of a polynomial (besides the QUADRATIC case). It attempts to solve the equation p(x)=0 for all the values of x (the roots).

For an overview of the Weierstrass method used in the L_POLYROOTS function, see the article on the Weierstrass or Durand-Kerner method. It is an iterative procedure that will require the user to specify the number of iterations (default is 20).

It is important to check your result because if the function reaches the maximum number of iterations without converging, you will not know that unless you check the results. The return_abs option will cause the function to return both the roots and a corresponding vector IMABS(p(x)) evaluated at those roots.

You can also use the L_POLYVAL function to check the value of p(x) at the roots like this:

=LET(
    coeffs, {1,2,3,4,5},
    roots, L_POLYROOTS(coeffs),
    check, L_POLYVAL(coeffs,roots),
    VSTACK("Roots",roots,"Check",check)
)

Result:

Roots
0.287815479557648-1.41609308017191i
0.287815479557648+1.41609308017191i
-1.28781547955765+0.857896758328489i
-1.28781547955765-0.857896758328491i
Check
-9.76996261670138E-15+1.95399252334028E-14i
-9.76996261670138E-15-1.95399252334028E-14i
-2.04281036531029E-14i
-1.95399252334028E-14+9.76996261670138E-15i

Specifying Iterations and/or Epsilon

At each iteration, the value of p(x) is computed for each x. Because the result may be a complex number, IMABS(p(x)) is used to compute the absolute value of z, or SQRT(a^2+b^2) where z=a+bi. If the maximum value of |p(x)| is less than epsilon, the algorithm ends the iteration and returns the result. This allows you to specify a large number of iterations without requiring the function to continue to run when it is no longer necessary.

Notice that the "Check" values in the above example are very close to zero (all less than 1E-12 or 0.000000000001). In this case, IMABS(check) is approximately 2E-14 for each of the roots.

Note

We eventually want to implement the Francis QR algorithm, but our current EIGENVALUE function must first be modified to handle complex numbers. For now, many practical cases can use the Weierstrass method.

Lambda Formula

This code for using L_POLYROOTS in Excel is provided under the License as part of the LAMBDA Library, but to use just this function, you may copy the following code directly into your spreadsheet.

Code for AFE Workbook Module (Excel Labs Add-in)

/**
* Find the real and complex roots of a polynomial using the Weierstrauss method
*/
L_POLYROOTS = LAMBDA(coeffs,[iterations],[epsilon],[xstart],[return_abs],
LET(doc,"https://www.vertex42.com/lambda/polyroots.html",
    iterations,IF(ISBLANK(iterations),10,iterations),
    xstart,IF(ISBLANK(xstart),"0.4+0.9i",xstart),
    return_abs,IF(ISBLANK(return_abs),FALSE,return_abs),
    n,COLUMNS(coeffs)-1,
    xo,TRANSPOSE(IMPOWER(xstart,SEQUENCE(1,n,0))),
    result,REDUCE(xo,SEQUENCE(iterations),LAMBDA(acc,k,
    IF(AND(NOT(ISBLANK(epsilon)),k>1,IFERROR(MAX(INDEX(acc,0,3))<epsilon,FALSE)),acc,
        LET(
            fo,IF(k>1,INDEX(acc,0,2),L_POLYVAL(coeffs,INDEX(acc,0,1))),
            xo,INDEX(acc,0,1),
            xn,MAKEARRAY(n,1,LAMBDA(i,j,
                IMSUB(INDEX(xo,i),IMDIV(INDEX(fo,i),
                  IMPRODUCT(IMSUB(INDEX(xo,i),
                  DROP(L_CIRCSHIFT(xo,1-i,1),1)))
                ))
            )),
            fn,L_POLYVAL(coeffs,xn),
            fabs,BYROW(fn,LAMBDA(value,IMABS(value))),
            HSTACK(xn,fn,fabs)
        )
    ))),
    IF(return_abs=TRUE,CHOOSECOLS(result,{1,3}),INDEX(result,0,1))
));

Named Function for Google Sheets

Name: L_POLYROOTS
Description: Find all real and complex roots of a polynomial using the Weierstrass method
Arguments: coeffs, iterations, epsilon, xstart, return_abs
Function:

=LET(doc,"https://www.vertex42.com/lambda/polyroots.html",
    iterations,IF(ISBLANK(iterations),10,iterations),
    xstart,IF(ISBLANK(xstart),"0.4+0.9i",xstart),
    return_abs,IF(ISBLANK(return_abs),FALSE,return_abs),
    n,COLUMNS(coeffs)-1,
    xo,ARRAYFORMULA(TRANSPOSE(IMPOWER(xstart,SEQUENCE(1,n,0)))),
    result,REDUCE(xo,SEQUENCE(iterations),LAMBDA(acc,k,
    IF(AND(NOT(ISBLANK(epsilon)),k>1,IFERROR(MAX(INDEX(acc,0,3))<epsilon,FALSE)),acc,
        LET(
            fo,IF(k>1,INDEX(acc,0,2),L_POLYVAL(coeffs,INDEX(acc,0,1))),
            xo,INDEX(acc,0,1),
            xn,MAKEARRAY(n,1,LAMBDA(i,j,
                IMSUB(INDEX(xo,i),IMDIV(INDEX(fo,i),IMPRODUCT(ARRAYFORMULA(IMSUB(INDEX(xo,i),L_DROP(L_CIRCSHIFT(xo,1-i,1),1,0)))))
            )),
            fn,L_POLYVAL(coeffs,xn),
            fabs,BYROW(fn,LAMBDA(value,IMABS(value))),
            HSTACK(xn,fn,fabs)
        )
    ))),
    IF(return_abs=TRUE,CHOOSECOLS(result,{1,3}),INDEX(result,0,1))
)

L_POLYROOTS Examples

Example 1: Accuracy Questionable

To demonstrate the significance of the possible errors with this function, consider the polynomial with coefficients {1,0,0,0,0} which has the roots 0,0,0,0 and the polynomial {1,4,6,4,1} which has the roots -1,-1,-1,-1. These are both 4th degree polynomials.

Test: Copy and Paste this LET function into a cell
=LET(
    coeffs, {1,0,0,0,0},
    n, 42,
    epsilon, 1E-12,
    L_POLYROOTS(coeffs,n,epsilon)
)

Result:
0.000607-0.000442i
0.000442+0.000607i
-0.000607+0.000442i
-0.000442-0.000607i

=LET(
    coeffs, {1,4,6,4,1},
    n, 42,
    epsilon, 1E-12,
    L_POLYROOTS(coeffs,n,epsilon)
)

Result:
-1.000668-0.000724i
-1.000726+0.000669i
-0.999330+0.000725i
-0.999275-0.000670i

In these cases, IMABS(p(x)) is less than 1E-12, but keep in mind that the higher degree polynomials can propogate small differences dramatically. In the above two examples, the error in calculation of the roots is approximately 0.001. Note that if x=0.001, then x^4 = 1E-12 or 0.000000000001.

It is important to remember that epsilon in this function is a measure of how close p(x) is to zero at the roots and not a measure of the number of significant digits in the roots themselves.

If you leave epsilon blank, then the function runs through all iterations instead of stopping at |p(x)|<ε.

Example 2

The polynomial x^4+2x^3+3x^2=0 has coefficients {1,2,3,0,0} and the roots (0, 0, -1+SQRT(2)i, -1-SQRT(2)i). Run L_POLYROOTS using 100 iterations for this example.

Test: Copy and Paste this LET function into a cell
=LET(
    coeffs, {1,2,3,0,0},
    n, 100,
    L_POLYROOTS(coeffs,n)
)

Result:
1.73809752232799E-31-1.49432331814218E-30i
-1+1.4142135623731i
-1.738097522328E-31+1.49432331814218E-30i
-1-1.4142135623731i

Check using L_POLYVAL(coeffs,roots):
-6.60837704751668E-60-1.55836779409192E-60i
2.93098878501041E-14+4.08562073062058E-14i
-6.60837704751668E-60-1.55836779409193E-60i
2.93098878501041E-14-4.08562073062058E-14i

In this case, SQRT(2) is accurate to the 15th digit (the maximum precision possible with the text-based complex numbers), and the zeros are indeed very close to zero.

Change History

3/05/2024 - v1.0.6 - Updated to use the new L_POLYVAL instead of IM_POLYVAL.

References & Resources

Disclaimer: This article is meant for educational purposes only. See the License regarding the LAMBDA code, and the site Terms of Use for the documentation.

L_POLYROOTS