RootFinding Package Updates - New Features in Maple 2019 - Maplesoft

What's New in Maple 2019

RootFinding Package Updates in 2019



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));
 

[x = -1.430647445, x = 3.652088914]  
 

  • 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 {F(x, y) = 0, G(x, y) = 0}:
 

> 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]);
 

Plot_2d
 

  • 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);
 

[-4.795738156, -3.854101966, -1.739664347, 1.250000000, 2.854101966, 3.227646598, 4.307755905, 7.500000000]
[-4.795738156, -3.854101966, -1.739664347, 1.250000000, 2.854101966, 3.227646598, 4.307755905, 7.500000000]
 
 

  • 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 := RootOf(`+`(`*`(`^`(_Z, 4)), `-`(`*`(`^`(_Z, 3))), `-`(`*`(27, `*`(`^`(_Z, 2)))), `*`(28, `*`(_Z)), 116), -4.812 .. -4.758)  
 

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

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

> Isolate(fa);
 

[y = -0.2992556510e-1]  
 

> Isolate(ga);
 

[y = -0.2992556510e-1, y = 0.2992556510e-1]  
 

  • Indeed, there is a common solution close to x = evalf[3](a), y = -0.3e-1:
 

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

true  
 

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

> b := candidates[4];
 

b := `/`(5, 4)  
 

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

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

> Isolate(fb);
 

[y = -5.845766191, y = .8624794396, y = 4.983286752]  
 

> Isolate(gb);
 

[y = -4.630064794, y = 4.630064794]  
 

  • This clearly is a spurious candidate; the roots of F(1.25, y) and G(1.25, y) 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]);
 

[[x = 3.227646598, y = -3.547153427], [x = -3.854101966, y = -2.854101966], [x = -4.795738156, y = -0.2992556510e-1], [x = 4.307755905, y = 2.107899206], [x = 2.854101966, y = 3.854101966], [x = -1.73...
[[x = 3.227646598, y = -3.547153427], [x = -3.854101966, y = -2.854101966], [x = -4.795738156, y = -0.2992556510e-1], [x = 4.307755905, y = 2.107899206], [x = 2.854101966, y = 3.854101966], [x = -1.73...
[[x = 3.227646598, y = -3.547153427], [x = -3.854101966, y = -2.854101966], [x = -4.795738156, y = -0.2992556510e-1], [x = 4.307755905, y = 2.107899206], [x = 2.854101966, y = 3.854101966], [x = -1.73...
 
 

  • 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 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) polynomials 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);
 

[-4.795738156, -3.854101966, -1.739664347, 1.285398164, 2.854101966, 3.227646598, 4.307755905, 7.535398164]
[-4.795738156, -3.854101966, -1.739664347, 1.285398164, 2.854101966, 3.227646598, 4.307755905, 7.535398164]
 
 

> a := candidates[1];
 

a := RootOf(`+`(`*`(`^`(_Z, 4)), `-`(`*`(`^`(_Z, 3))), `-`(`*`(27, `*`(`^`(_Z, 2)))), `*`(28, `*`(_Z)), 116), -4.795738156)  
 

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

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

> Isolate(fa);
 

[y = -0.2992556510e-1]  
 

> Isolate(ga);
 

[y = -0.2992556510e-1, y = 0.2992556510e-1]  
 

> b := candidates[4];
 

b := `+`(`*`(`/`(1, 4), `*`(Pi)), `/`(1, 2))  
 

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

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

> Isolate(fb);
 

[y = -5.827352781, y = .8577160795, y = 4.969636701]  
 

> Isolate(gb);
 

[y = -4.620362709, y = 4.620362709]  
 

  • 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));
 

false, false, false  
 

  • Indeed, gb evaluated over all isolating intervals for 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));
 

true  
 

  • shows that ga, evaluated at the isolating intervals for the root of 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);
 

apx := -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));
 

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.
 

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.15e-1  
 

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

0.27e-1  
 

> time(Isolate(mig(50)));
 

.209  
 

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

.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.56e-1  
 

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

.239  
 

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

.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