Problem 492979 · hard · Level 04 Non-Linear Data Structures

Finding the Poles: a Root Finder in Plain Python

polynomial roots · Durand-Kerner method · iteration · complex numbers · poles

Whether a filter or a control loop is stable depends on its poles, the values of z where A(z) = 0. Multiply A(z) = a[0] + a[1] * z^-1 + ... + a[N] * z^-N by z^N and you get the ordinary polynomial a[0] * z^N + a[1] * z^(N-1) + ... + a[N]: the same list, read as powers of z with the highest power first. So the poles are the roots of the list a. For degree 1 or 2 there are formulas, but a fourth-order drone filter needs a numerical method, and you will write one now.

Find the roots of p[0] * z^n + p[1] * z^(n-1) + ... + p[n] with the Durand-Kerner method. It improves guesses for all n roots together:

  1. Divide every coefficient by p[0], so the leading coefficient is 1 (same roots). Call the polynomial P.
  2. Start from the guesses z_i = (0.4 + 0.9j) ** (i + 1) for i = 0, ..., n-1 (complex, all different, and none of them on the real axis).
  3. One sweep: for each i in turn, replace z_i by z_i - P(z_i) / ((z_i - z_0) * (z_i - z_1) * ...), the product running over every j != i, using the newest values of the other guesses. (If two guesses ever become exactly equal, the product is 0: use 1e-12 in its place.)
  4. Repeat sweeps until the largest change in a sweep is below 1e-12, or after 500 sweeps.

Why it works: if the other guesses were exactly the other roots, P(z) / (product of (z - z_j)) would be just z - root_i, so one step lands on the root; with all guesses improving together, the method closes in on all the roots at once.

Write poly_roots(p) that returns the n roots as a list of pairs (re, im), in any order. The answer is checked by multiplying (z - r) over your roots and comparing the result with p / p[0], coefficient by coefficient, with a relative tolerance of 1e-6. Real roots may come out with a tiny imaginary part such as 1e-17; that is fine. Press Run with plot([r[0] for r in roots], [r[1] for r in roots], kind="scatter") to see where they lie.

Examples

Input:  p = [1, -3, 2]
Output: [(1.0, 0.0), (2.0, 0.0)]          (any order)
Explanation: z^2 - 3z + 2 = (z - 1)(z - 2).

Input:  p = [1, 0, 0.81]
Output: [(0.0, 0.9), (0.0, -0.9)]         (any order)
Explanation: the notch filter's poles: z^2 = -0.81.

Input:  p = [5]
Output: []
Explanation: a constant has no roots.

Constraints

  • p[0] != 0, len(p) <= 11 (degree at most 10), real coefficients
  • the roots are all different, though some lie close together (as little as 0.04 apart)

Goals

  • Find all the roots of a polynomial, real and complex, by iteration
  • Use the Durand-Kerner update, which improves every root guess at once
  • Read the denominator list of H(z) as the polynomial whose roots are the poles
Starting Python…