Real Root Finding - Maple Help

Home : Support : Online Help : System : Information : Updates : Maple 2019 : Real Root Finding

Description

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

Isolating roots of univariate polynomials with 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}{=}{-1.430647445}{,}{x}{=}{3.652088914}\right]$ (1)
 • 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[{-4.795738156}{,}{-3.854101966}{,}{-1.739664347}{,}{1.250000000}{,}{2.854101966}{,}{3.227646598}{,}{4.307755905}{,}{7.500000000}\right]$ (2)
 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}{,}{-4.7957383}{..}{-4.7957372}\right)$ (3)
 > fa := subs(x=a, F):
 > ga := subs(x=a, G):
 > Isolate(fa);
 $\left[{y}{=}{-0.02992556510}\right]$ (4)
 > Isolate(ga);
 $\left[{y}{=}{-0.02992556510}{,}{y}{=}{0.02992556510}\right]$ (5)
 Indeed, there is a common solution close to $x={\mathrm{evalf}}_{3}\left(a\right)$, $y=-0.03$:
 > is(RootOf(fa, y, -0.03) = RootOf(ga, y, -0.03));
 ${\mathrm{true}}$ (6)
 However, let's look at the candidate at $b=1.25$:
 > b := candidates[4];
 ${b}{≔}\frac{{5}}{{4}}$ (7)
 > fb := subs(x=b, F):
 > gb := subs(x=b, G):
 > Isolate(fb);
 $\left[{y}{=}{-5.845766191}{,}{y}{=}{0.8624794396}{,}{y}{=}{4.983286752}\right]$ (8)
 > Isolate(gb);
 $\left[{y}{=}{-4.630064794}{,}{y}{=}{4.630064794}\right]$ (9)
 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}{=}{3.227646598}{,}{y}{=}{-3.547153427}\right]{,}\left[{x}{=}{-3.854101966}{,}{y}{=}{-2.854101966}\right]{,}\left[{x}{=}{-4.795738156}{,}{y}{=}{-0.02992556510}\right]{,}\left[{x}{=}{4.307755905}{,}{y}{=}{2.107899206}\right]{,}\left[{x}{=}{2.854101966}{,}{y}{=}{3.854101966}\right]{,}\left[{x}{=}{-1.739664347}{,}{y}{=}{4.469179786}\right]\right]$ (10)
 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]);
 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[{-4.795738156}{,}{-3.854101966}{,}{-1.739664347}{,}{1.285398164}{,}{2.854101966}{,}{3.227646598}{,}{4.307755905}{,}{7.535398164}\right]$ (11)
 > a := candidates[1];
 ${a}{≔}{\mathrm{RootOf}}{}\left({{\mathrm{_Z}}}^{{4}}{-}{{\mathrm{_Z}}}^{{3}}{-}{27}{}{{\mathrm{_Z}}}^{{2}}{+}{28}{}{\mathrm{_Z}}{+}{116}{,}{-4.795738156}\right)$ (12)
 > fa := subs(x=a, F):
 > ga := subs(x=a, G):
 > Isolate(fa);
 $\left[{y}{=}{-0.02992556510}\right]$ (13)
 > Isolate(ga);
 $\left[{y}{=}{-0.02992556510}{,}{y}{=}{0.02992556510}\right]$ (14)
 > b := candidates[4];
 ${b}{≔}\frac{{\mathrm{\pi }}}{{4}}{+}\frac{{1}}{{2}}$ (15)
 > fb := subs(x=b, F):
 > gb := subs(x=b, G):
 > Isolate(fb);
 $\left[{y}{=}{-5.827352781}{,}{y}{=}{0.8577160795}{,}{y}{=}{4.969636701}\right]$ (16)
 > Isolate(gb);
 $\left[{y}{=}{-4.620362709}{,}{y}{=}{4.620362709}\right]$ (17)
 • 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}}$ (18)
 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}}$ (19)
 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}}{≔}{-4.795738156}$ (20)
 > 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}}$ (21)
 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.

Performance improvements for root isolation and root refinement

 • 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}^{-100}$:
 > with(RootFinding):
 > mig := n -> x^n - (nextprime(2^100)*x^2 - 1)^2:
 > time(Isolate(mig(10)));
 ${0.015}$ (22)
 > time(Isolate(mig(10), method=RS));
 ${0.027}$ (23)
 > time(Isolate(mig(50)));
 ${0.209}$ (24)
 > time(Isolate(mig(50), method=RS));
 ${0.877}$ (25)
 > time(Isolate(mig(100)));
 ${1.125}$ (26)
 > time(Isolate(mig(100), method=RS));
 ${7.609}$ (27)
 > time(Isolate(mig(200)));
 ${7.469}$ (28)
 > timelimit(600, Isolate(mig(200), method=RS));
 ${"time expired"}$ (29)
 • 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}$ (30)
 > time(Isolate(f, digits=100, method=RS));
 ${0.239}$ (31)
 > time(Isolate(f, digits=1000));
 ${0.293}$ (32)
 > time(Isolate(f, digits=1000, method=RS));
 ${107.328}$ (33)
 > time(Isolate(f, digits=10000));
 ${10.824}$ (34)
 > timelimit(600, Isolate(f, digits=10000, method=RS));
 ${"time expired"}$ (35)