The RootFinding:-Isolate command can now be used to isolate the roots of univariate polynomials with arbitrary real coefficients.

Prior to 2019, RootFinding:-Isolate could only determine the roots of polynomials with rational or float coefficients. This restriction is now lifted for univariate polynomials:

with(RootFinding):

Isolate(sqrt(2)*x^2 - Pi*x - exp(2));

$\left[x=\mathrm{-1.430647445}\,x=3.652088914\right]$

In particular, Isolate can now be used to find roots of polynomials with algebraic coefficients. We illustrate this in an example where we manually study the real solutions of a bivariate equation system of the form $\left\{F\left(x\,y\right)=0\,G\left(x\,y\right)=0\right\}$:

F := (2*x^2*y - 2*x^2 - 3*x + y^3 - 33*y + 32) * ((x-2)^2 + y^2 + 3):

G := (x^2 + y^2 - 23) * (x^2 + y^2 + 2):

The common roots of both polynomials are the intersections of their corresponding algebraic curves:

plots[implicitplot]([F = 0, G = 0], x=-16..16, y=-7..6, color=["Teal", "Red"], gridrefine=2, scaling=constrained, size=[0.7,0.35]);

Elimination theory for algebraic equation systems tells us that all x-coordinates of the common solutions are roots of the resultant polynomial of $F$ and $G$ with respect to $y$:

R := resultant(F, G, y):

This resultant is a univariate polynomial with irrational roots, some of which may be complex. The roots are candidate values for the x-coordinate of simultaneous solutions. Note that we store symbolic expressions for the solutions, not just approximations:

candidates := sort([RealDomain[solve](R)], key=evalf):

evalf(candidates);

$\left[\mathrm{-4.795738156}\,\mathrm{-3.854101966}\,\mathrm{-1.739664347}\,1.250000000\,2.854101966\,3.227646598\,4.307755905\,7.500000000\right]$

However, some of the candidates might be spurious. We can use the new interface for RootFinding:-Isolate to determine the roots along the fibers of $F$ and $G$ when we substitute the candidates:

a := candidates[1];

$a≔\mathrm{RootOf}\left({\mathrm{\_Z}}^4-{\mathrm{\_Z}}^3-27{\mathrm{\_Z}}^2+28\mathrm{\_Z}+116\,\mathrm{-4.7957383}..\mathrm{-4.7957372}\right)$

fa := subs(x=a, F):

ga := subs(x=a, G):

Isolate(fa);

$\left[y=\mathrm{-0.02992556510}\right]$

Isolate(ga);

$\left[y=\mathrm{-0.02992556510}\,y=0.02992556510\right]$

Indeed, there is a common solution close to $x={\mathrm{evalf}}_3\left(a\right)$, $y=\mathrm{-0.03}$:

is(RootOf(fa, y, -0.03) = RootOf(ga, y, -0.03));

$\mathrm{true}$

However, let's look at the candidate at $b=1.25$:

b := candidates[4];

$b≔\frac{5}{4}$

fb := subs(x=b, F):

gb := subs(x=b, G):

Isolate(fb);

$\left[y=\mathrm{-5.845766191}\,y=0.8624794396\,y=4.983286752\right]$

Isolate(gb);

$\left[y=\mathrm{-4.630064794}\,y=4.630064794\right]$

This clearly is a spurious candidate; the roots of $F\left(1.25\,y\right)$ and $G\left(1.25\,y\right)$ are distinct.

The example code above can help to filter spurious solutions, but it is not a complete solver for bivariate systems; it allows to filter out suspicious candidates, but does not validate all solutions. Such verification is provided by the multivariate solvers in the RootFinding package:

Isolate([F,G], [x,y]);

$\left[\left[x=\mathrm{-4.795738156}\,y=\mathrm{-0.02992556510}\right]\,\left[x=\mathrm{-3.854101966}\,y=\mathrm{-2.854101966}\right]\,\left[x=\mathrm{-1.739664347}\,y=4.469179786\right]\,\left[x=2.854101966\,y=3.854101966\right]\,\left[x=3.227646598\,y=\mathrm{-3.547153427}\right]\,\left[x=4.307755905\,y=2.107899206\right]\right]$

However, the multivariate polynomial solver requires coefficients of type numeric (that is, rationals or floats). Consider the case where we slightly change $F$ by replacing $3$ with $\mathrm{\pi }$ in the last term:

F := (2*x^2*y - 2*x^2 - 3*x + y^3 - 33*y + 32) * ((x-2)^2 + y^2 + Pi):

Isolate([F,G], [x,y]);

Error, (in RootFinding:-Isolate) polynomial(s) to be solved must have numeric coefficients

The more naive approach above, while uncertified, will work nevertheless:

R := resultant(F, G, y):

candidates := sort([RealDomain[solve](R)], key=evalf):

evalf(candidates);

$\left[\mathrm{-4.795738156}\,\mathrm{-3.854101966}\,\mathrm{-1.739664347}\,1.285398164\,2.854101966\,3.227646598\,4.307755905\,7.535398164\right]$

a := candidates[1];

$a≔\mathrm{RootOf}\left({\mathrm{\_Z}}^4-{\mathrm{\_Z}}^3-27{\mathrm{\_Z}}^2+28\mathrm{\_Z}+116\,\mathrm{-4.795738156}\right)$

fa := subs(x=a, F):

ga := subs(x=a, G):

Isolate(fa);

$\left[y=\mathrm{-0.02992556510}\right]$

Isolate(ga);

$\left[y=\mathrm{-0.02992556510}\,y=0.02992556510\right]$

b := candidates[4];

$b≔\frac{\mathrm{\pi }}{4}+\frac{1}{2}$

fb := subs(x=b, F):

gb := subs(x=b, G):

Isolate(fb);

$\left[y=\mathrm{-5.827352781}\,y=0.8577160795\,y=4.969636701\right]$

Isolate(gb);

$\left[y=\mathrm{-4.620362709}\,y=4.620362709\right]$

Finally, we show how the combination of the constraints and output options of RootFinding:-Isolate can provide certified information, and will allow us to programmatically exclude spurious candidates even with irrational coefficients.

rts_fb, gb_at_rts_fb := Isolate(fb, constraints=[gb], output=interval):

contains_zero := iv -> evalb(iv[1] <= 0 and iv[2] >= 0):

seq(contains_zero(rhs(gb_at_rts_fb[i][1])), i=1..nops(rts_fb));

$\mathrm{false},\mathrm{false},\mathrm{false}$

Indeed, $\mathrm{gb}$ evaluated over all isolating intervals for $\mathrm{fb}$ does not contain zero, which confirms that $F$ and $G$ have no common zero at $x=b$. In contrast,

rts_fa, ga_at_rts_fa := Isolate(fa, constraints=[ga], output=interval):

seq(contains_zero(rhs(ga_at_rts_fa[i][1])), i=1..nops(rts_fa));

$\mathrm{true}$

shows that $\mathrm{ga}$, evaluated at the isolating intervals for the root of $\mathrm{fa}$, contains zero. This still does not validate the simultaneous zero of both systems, but is a strong hint. Techniques along these lines can serve to filter candidates numerically before trying time-consuming symbolic simplification and zero-testing, and can be used as cornerstones for complete solvers.

Note that the aforementioned routines crucially rely on inputs that are served as symbolic expressions rather than approximations:

apx := evalf(a);

$\mathrm{apx}≔\mathrm{-4.795738156}$

fapx := subs(x=apx, F):

gapx := subs(x=apx, G):

rts_fapx, gapx_at_rts_fapx := Isolate(fapx, constraints=[gapx], output=interval):

seq(contains_zero(rhs(gapx_at_rts_fapx[i][1])), i=1..nops(rts_fapx));

$\mathrm{false}$

Even a tiny perturbation of the candidate solution in $x$ will produce distinct roots of $F$ and $G$ in $y$. Thus, the direct handling of arbitrary real coefficients is not only convenient, but required for correctness.

The new default algorithm of Isolate also features vastly improved performance for ill-conditioned polynomials with clustered roots. The root finding method eventually converges quadratically to regions containing roots, rather than just linearly. For example, the following class of polynomials has a cluster of roots extremely close to $2^{\mathrm{-100}}$:

with(RootFinding):

mig := n -> x^n - (nextprime(2^100)*x^2 - 1)^2:

time(Isolate(mig(10)));

$0.015$

time(Isolate(mig(10), method=RS));

$0.027$

time(Isolate(mig(50)));

$0.209$

time(Isolate(mig(50), method=RS));

$0.877$

time(Isolate(mig(100)));

$1.125$

time(Isolate(mig(100), method=RS));

$7.609$

time(Isolate(mig(200)));

$7.469$

timelimit(600, Isolate(mig(200), method=RS));

$''time\; expired''$

The same technique allows even more dramatic improvements for root finding requests with high accuracy even on well-conditioned problems:

f := add(rand(-1. .. 1.)() * x^i, i=0..100):

time(Isolate(f, digits=100));

$0.056$

time(Isolate(f, digits=100, method=RS));

$0.239$

time(Isolate(f, digits=1000));

$0.293$

time(Isolate(f, digits=1000, method=RS));

$107.328$

time(Isolate(f, digits=10000));

$10.824$

timelimit(600, Isolate(f, digits=10000, method=RS));

$''time\; expired''$