3
$\begingroup$

I have a list of polynomials polys.

I want a set of $f(x) \in$ polys which have roots strictly in $\{S^1,0\}$, with $S^1$ being the unit sphere.

They are all degree $n$, so I have crudely picked out those $f(x)$ such that the sum of the norms of the roots of $f$ is $\leq n$.

rts = Table[Values[{ToRules[Roots[Part[polys, i] == 0, x]]}], {i, 1, Length[polys]}] 

rtsAbsSum = Flatten[Table[Sum[Abs[Part[Part[rts, i], j]], {j, 1, Length[Part[rts, i]]}], {i, 1, Length[polys]}]]

parts = Select[Range[Length[polys]], Part[rtsAbsSum, #] <= n &]

polysIWant = Table[Part[polys, i], {i, parts}]

Now the list polysIWant contains the polynomials I want. However, it also contains many false positives.

This code is crude. What its doing:

  1. Takes each polynomial of polys and replaces it with a list of its roots (multiplicity not counted.)
  2. Adds together the Abs of the elements of each list.
  3. Finds the polynomials whose sum of root norms is $\leq n$.

Is there an easier way? ie. Instead of summing roots, checking the inequality, etc., how can I write a program that tosses out the polynomials in polys with a root whose norm is neither $0$ nor $1$, and keeps the rest?

$\endgroup$
  • $\begingroup$ If you have a function roots that gets the roots of a polynomial into a list form of consistent dimension, then why not something like: Select[polys, Not[0!=Norm[roots[#]]!=1] &] $\endgroup$ – flinty Aug 6 at 18:08
  • $\begingroup$ I assume that by $S^1$ you mean that unit circle in the complex plane? So you are allowing complex roots? Just for clarification. $\endgroup$ – march Aug 6 at 20:34
  • $\begingroup$ Apparently there is something called the Jury stability criterion which you can use to determine if the roots of a polynomial are all inside the unit disk. $\endgroup$ – march Aug 6 at 20:42
  • 2
    $\begingroup$ Are "roots of undesirable magnitude" related to "rodents of unusual size"? $\endgroup$ – Daniel Lichtblau Aug 8 at 16:33
  • 2
    $\begingroup$ R.O.U.M.s? I do not think they exist. $\endgroup$ – march Aug 11 at 2:58
6
$\begingroup$

Here is an implementation of the Jury stability criterion (also known as the Jury-Marden criterion; there is also the Schur-Cohn criterion), which determines if the roots of a real-coefficient polynomial are all in the open unit disk. The function returns True if all of the roots of the polynomial are in the open unit disk and returns False otherwise.

juryStabilityCheck[poly_, var_] := juryStabilityCheck[Reverse@CoefficientList[poly, var]]
juryStabilityCheck[coefs_List] := Module[
   {lst = coefs},
   If[First@lst < 0, lst = -lst];
   While[lst != {} && First@lst > 0, lst = lst - Last@lst/First@lst Reverse@lst // Most];
   lst == {}
  ]

Here is a polynomial whose roots are all in the open unit disk:

poly = (z - 0.5) (z - I/2) (z + I/2) (z - 0.3 I) (z + 0.3 I);

If we run the check on it, we get:

juryStabilityCheck[poly, z]
(* True *)

It works for any scalar multiple of that list:

juryStabilityCheck[-10 poly, z]
(* True *)

If some roots are outside the unit disk:

poly = (z - 0.5) (z - 3 I/2) (z + 3 I/2) (z - 0.3 I) (z + 0.3 I);
juryStabilityCheck[poly, z]
(* False *)

If a root is on the unit circle (i.e. the boundary of the unit disk), then the check finds that they are outside the open unit disk:

poly = (z - 1) (z - I/2) (z + I/2) (z - 0.3 I) (z + 0.3 I);
juryStabilityCheck[poly, z]
(* False *)

Now, if you want to include polynomials that have a root on the unit circle, the quick fix I've come up with is to "enlarge" the circle slightly by "shrinking" the z's. You must choose a tolerance, and be okay with polynomials being included that have roots within this tolerance outside the unit circle.

Add to the code the definition

juryStabilityCheck[poly_, var_, tolerance_] := juryStabilityCheck[poly /. var -> var (1 + tolerance), var]

As an example:

poly = (z - 1) (z - I/2) (z + I/2) (z - 0.3 I) (z + 0.3 I) // Expand;
juryStabilityCheck[poly, z]
juryStabilityCheck[poly, z, 10^-9]
(* False *)
(* True *)

Of course, this will also include a polynomial whose roots are within the tolerance but outside the unit circle:

poly = (z - 1 + 10^-10) (z - I/2) (z + I/2) (z - 0.3 I) (z + 0.3 I);
juryStabilityCheck[poly, z, 10^-9]
(* True *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, this is excellent. I've made a one-character change to your code that I believe was a typo (at least, everything worked AFTER I made it.) $\endgroup$ – Descartes Before the Horse Aug 6 at 23:02
  • $\begingroup$ For searching purposes: this is also sometimes referred to as the "Jury-Marden" criterion. There is also the more classical Schur-Cohn criterion for doing the same test. $\endgroup$ – J. M.'s discontentment Aug 7 at 9:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.