RegularChains - New Features in Maple 2020 - Maplesoft

What's New in Maple 2020

RegularChains

The Regular Chains package is a collection of commands for solving systems of algebraic equations, inequations, and inequalities symbolically. It also enables users to manipulate and study the solutions of such systems. Maple's RegularChains package has received many updates since Maple 2019. Below you will find some improvements for Maple 2020.


restart; 1 

with(RegularChains); -1 

with(SemiAlgebraicSetTools); -1 

with(AlgebraicGeometryTools); -1 

with(ChainTools); -1


Quantifier Elimination 

Quantifier elimination is the process of taking a formula that involves statements of the form "there exists a value for variable x; such that..." or "for all values of variable x; , we have..." (represented with the quantifiers "there exists", or , and "for all", or , respectively), where the rest of the statement typically involves logical connectives such as "and" or "not" and polynomial formulas, and expressing this as an equivalent statement that does not include such quantifiers but depends only on the other variables that are not quantified over. These formulas are all interpreted over the real domain. 

An example is probably in order. The following is known as the Davenport-Heintz problem: 

 

This can be encoded as follows. 

input := `&E`([c]), `&A`([b, a]), `&implies`(`&or`(`&and`(a = d, b = c), `&and`(a = c, b = 1)), `*`(`^`(a, 2)) = b);  

Typesetting:-mprintslash([input := `&E`([c]), `&A`([b, a]), `&implies`(`&or`(`&and`(a = d, b = c), `&and`(a = c, b = 1)), `*`(`^`(a, 2)) = b)], [`&E`([c]), `&A`([b, a]), `&implies`(`&or`(`&and`(a = d,...
 

The solution is very simple. The QuantifierElimination command is in the SemiAlgebraicSetTools subpackage of RegularChains. 

QuantifierElimination(input);  

`&or`(`+`(d, `-`(1)) = 0, `+`(d, 1) = 0)
 

It is easy to see by hand that these are indeed solutions to the given problem (with e.g. c = 1; in both cases), but it is not obvious that there are no other solutions. 

 

We can use this functionality to find all cases where the quadratic formula has roots. We get a more legible answer by using some optional arguments to the command: 

R := PolynomialRing([x, c, b, a]); -1 

QuantifierElimination(`&E`([x]), `+`(`*`(a, `*`(`^`(x, 2))), `*`(b, `*`(x)), c) = 0, R, output = rootof);  

`&or`(`&and`(`<`(a, 0), `<=`(`+`(`/`(`*`(`/`(1, 4), `*`(`^`(b, 2))), `*`(a))), c)), `&and`(a = 0, `<>`(b, 0)), `&and`(a = 0, b = 0, c = 0), `&and`(`<`(0, a), `<=`(c, `+`(`/`(`*`(`/`(1, 4), `*`(`^`(b, ...
 

Cylindrical Algebraic Decomposition 

Cylindrical Algebraic Decomposition is a partition of the real, n; -dimensional space into finitely many connected semi-algebraic subsets, called cells, such that any two cells are cylindrically arranged, that is, their projection onto a lower-dimensional space is either identical or disjoint. Here, a semi-algebraic set is a set given by polynomial equations, inequations, and inequalities. 

The CylindricalAlgebraicDecompose command finds such a decomposition where the number of solutions for a given semi-algebraic system stays constant on each cell. A simple example is given as follows: 

R := PolynomialRing([y, x]); -1 

cad := CylindricalAlgebraicDecompose([`<=`(`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))), 1)], R);  

Typesetting:-mprintslash([cad := c_a_d], [c_a_d])
 

Display(cad, R);  

[PIECEWISE([y = y, ``], [`<`(x, -1), ``]), PIECEWISE([`<`(y, 0), ``], [x = -1, ``]), PIECEWISE([y = 0, ``], [x = -1, ``]), PIECEWISE([`<`(0, y), ``], [x = -1, ``]), PIECEWISE([`<`(y, `+`(`-`(`*`(`^`(`...
[PIECEWISE([y = y, ``], [`<`(x, -1), ``]), PIECEWISE([`<`(y, 0), ``], [x = -1, ``]), PIECEWISE([y = 0, ``], [x = -1, ``]), PIECEWISE([`<`(0, y), ``], [x = -1, ``]), PIECEWISE([`<`(y, `+`(`-`(`*`(`^`(`...
[PIECEWISE([y = y, ``], [`<`(x, -1), ``]), PIECEWISE([`<`(y, 0), ``], [x = -1, ``]), PIECEWISE([y = 0, ``], [x = -1, ``]), PIECEWISE([`<`(0, y), ``], [x = -1, ``]), PIECEWISE([`<`(y, `+`(`-`(`*`(`^`(`...
 

A more compact output is given with the option output = rootof; . 

output = rootof
 

CylindricalAlgebraicDecompose([`<=`(`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))), 1)], R, output = rootof);  

[[And(`<=`(-1, x), `<=`(x, 1)), And(`<=`(`+`(`-`(`*`(`^`(`+`(`-`(`*`(`^`(x, 2))), 1), `/`(1, 2))))), y), `<=`(y, `*`(`^`(`+`(`-`(`*`(`^`(x, 2))), 1), `/`(1, 2)))))]]
 

The CylindricalAlgebraicDecompose; command has been available for a long time. In Maple 2020, it has some new options. One of them is the optimization option. In the simplest form, it is simply a true/false option (with default value true). With this option, Maple tries to reduce the number of cells returned. This example uses the Info command, which returns a list with some information about each cell found; the number of elements of this list is the number of cells. 

F := [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = 1, `*`(`^`(y, 2)) = x]; -1 

cad_nonopt := CylindricalAlgebraicDecompose(F, R, optimization = false); -1 

numelems(Info(cad_nonopt, R)); 1 

53
 

cad_opt := CylindricalAlgebraicDecompose(F, R, optimization = true); -1 

numelems(Info(cad_opt, R)); 1 

9
 

Another change is that the default algorithm used under the hood has changed. The previous algorithm can still be selected using the method option. The default value for the method option is incremental; the other possible value is recursive. For large systems, the incremental method is typically substantially faster. In this small example, the same difference in speed occurs. 

CodeTools:-Usage(CylindricalAlgebraicDecompose(F, R, method = recursive)); -1 

memory used=14.99MiB, alloc change=0 bytes, cpu time=130.00ms, real time=129.00ms, gc time=0ns
 

CodeTools:-Usage(CylindricalAlgebraicDecompose(F, R, method = incremental)); -1 

memory used=1.78MiB, alloc change=0 bytes, cpu time=15.00ms, real time=15.00ms, gc time=0ns
 

Computing intersection multiplicities 

The new command IntersectionMultiplicity returns the intersection multiplicity of a space curve at every point of that curve defined by a regular chain. This is shown in the following example. 

R := PolynomialRing([x, y, z]); -1 

F := [`+`(`*`(`^`(x, 2)), y, z, `-`(1)), `+`(`*`(`^`(y, 2)), x, z, `-`(1)), `+`(`*`(`^`(z, 2)), x, y, `-`(1))]; -1 

We consider the points of intersection of all three surfaces. First, we compute these as regular chains. 

dec := Triangularize(F, R); -1 

Display(dec, R); 1 

[PIECEWISE([`+`(x, `-`(z)) = 0, ``], [`+`(y, `-`(z)) = 0, ``], [`+`(`*`(`^`(z, 2)), `*`(2, `*`(z)), `-`(1)) = 0, ``]), PIECEWISE([x = 0, ``], [y = 0, ``], [`+`(z, `-`(1)) = 0, ``]), PIECEWISE([x = 0, ...
 

By inspection, these solutions correspond to the points [x, y, z] = `&+-`([-1, -1, -1], [sqrt(2), sqrt(2), sqrt(2)]); (two points), [0, 0, 1]; , [0, 1, 0]; , and [1, 0, 0]; , respectively. 

The three surfaces look as follows. We show a larger plot that includes all intersection points of all three surfaces, and two close ups: one of the point [`+`(`-`(1), sqrt(2)), `+`(`-`(1), sqrt(2)), `+`(`-`(1), sqrt(2))]; and one of the point [0, 0, 1]; . 

plots:-implicitplot3d(F, x = -3 .. 1, y = -3 .. 1, z = -3 .. 1, grid = [20, 20, 20], color = [red, green, blue], style = surface);

Plot_2d
 

plots:-implicitplot3d(F, x = .35 .. .45, y = .35 .. .45, z = .35 .. .45, grid = [20, 20, 20], color = [red, green, blue], style = surface);

Plot_2d

plots:-implicitplot3d(F, x = -.1 .. .1, y = -.1 .. .1, z = .9 .. 1.1, grid = [20, 20, 20], color = [red, green, blue], style = surface); &

Plot_2d

We see that in the second close up, the pairwise intersection curves appear to be locally tangent, which is not the case in the first close up. This is confirmed by the IntersectionMultiplicity command. 

map(IntersectionMultiplicity, dec, F, R);  

Typesetting:-mprintslash([[[[1, regular_chain]], [[2, regular_chain]], [[2, regular_chain]], [[2, regular_chain]]]], [[[[1, Record(property = isPrime, polynomials = [[x, 1, 1, `+`(, `-`(z)), isPrimeIr...
 

This confirms that the intersection multiplicity of all points in the first regular chain is 1, whereas it is 2 for the others. 

Tangent cones 

The new command TangentCone computes a tangent cone to a given space curve at a given set of points. Consider the following example. 

R := PolynomialRing([x, y, z]); -1 

rc := Chain([`+`(z, `-`(1)), y, x], Empty(R), R); -1 

rc; , defined above, represents the point [x, y, z] = [0, 0, 1]; . 

F := [`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2)), `*`(`^`(z, 2)), `-`(1)), `+`(`*`(`^`(x, 2)), `-`(`*`(`^`(y, 2))), `-`(`*`(z, `*`(`+`(z, `-`(1))))))]; -1 

We will find the tangent cone to the intersection of the surfaces defined in F; at the point defined in rc; . First let us examine these surfaces visually. 

solutions := map(solve, F, z);  

Typesetting:-mprintslash([solutions := [`*`(`^`(`+`(`-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2))), 1), `/`(1, 2))), `+`(`-`(`*`(`^`(`+`(`-`(`*`(`^`(x, 2))), `-`(`*`(`^`(y, 2))), 1), `/`(1, 2))))), `+`(`/`(...
 

We are interested in the first and third solution. We create a plot of these two surfaces, call it p; , and display it. 

p := plot3d(solutions[[1, 3]], x = -.5 .. .5, y = -.5 .. .5, color = [red, blue]); -1; p 

Plot_2d

We now compute the tangent cone. 

tc := TangentCone(rc, F, R, 'equations', [X, Y, Z]);  

Typesetting:-mprintslash([tc := [[[`+`(Z, `-`(1)), `+`(`*`(3, `*`(`^`(X, 2))), `-`(`*`(`^`(Y, 2))))], regular_chain]]], [{[[`+`(Z, `-`(1)), `+`(`*`(3, `*`(`^`(X, 2))), `-`(`*`(`^`(Y, 2))))], Record(pr...
 

We see that it corresponds to the lines given by z = 1; and y = `*`(`&+-`(sqrt(3)), `*`(x)); . We show this in the same plot. 

lines := seq(plottools:-line([`+`(`-`(`/`(`*`(`/`(1, 2)), `*`(sqrt(3))))), `+`(`-`(`*`(`/`(1, 2), `*`(s)))), 1], [`+`(`/`(`*`(`/`(1, 2)), `*`(sqrt(3)))), `+`(`*`(`/`(1, 2), `*`(s))), 1], color = green... 

plots:-display(p, lines);  

Plot_2d

Limit points of implicitly defined curves 

Another new feature of the RegularChains library is its ability to determine the geometry of curves specified as the solution set of some equations. In particular, the library can find where the leading coefficients of these equations vanish, and whether the curve is locally real at those points or not. An example follows. 

Consider the equations `+`(`*`(`^`(z, 5)), `-`(`*`(`^`(y, 3))), `*`(`^`(y, 2))) = 0; and `+`(`*`(x, `*`(`^`(z, 4))), `*`(`^`(y, 3)), `-`(`*`(`^`(y, 2)))) = 0; . The surfaces defined by them look as follows; this is purely for illustration and the method of drawing it won't be explained here. The first equation's surface is rendered in red, the second in blue, and the curve that is their intersection in green. You can also see a transparent black plane at z = 0; , and the two intersection points of the curve with this plane in black; these are explained below. 

Plot_2d

If we view the variables as having the order `and`(`>`(x, y), `>`(y, z)); , then the leading coefficients of these equations are 1 and `*`(`^`(z, 4)); , respectively. The first coefficient doesn't vanish; the second vanishes at z = 0; . The LimitPoints command can find the intersections of the curve, that is, the intersection of the red and blue surfaces. By default, this computation is done over the complex numbers. 

p1 := `+`(`*`(`^`(z, 5)), `-`(`*`(`^`(y, 3))), `*`(`^`(y, 2))); -1 

p2 := `+`(`*`(x, `*`(`^`(z, 4))), `*`(`^`(y, 3)), `-`(`*`(`^`(y, 2)))); -1 

R := PolynomialRing([x, y, z]); -1 

rc := Chain([p1, p2], Empty(R), R); -1 

Display(LimitPoints(rc, R), R);  

[PIECEWISE([x = 0, ``], [y = 0, ``], [z = 0, ``]), PIECEWISE([x = 0, ``], [`+`(y, `-`(1)) = 0, ``], [z = 0, ``])]
 

We see that there are two such limit points. Indeed, if we evaluate the polynomials defining the equations at these points, they vanish. These are the two points shown in the graph above. 

eval([p1, p2], {x = 0, y = 0, z = 0});  

[0, 0]
 

eval([p1, p2], {x = 0, y = 1, z = 0});  

[0, 0]
 

It looks like there is a cusp at the origin: the curve is not smooth there. We can confirm that this is indeed the case: 

Display(LimitPoints(rc, R, coefficient = real), R);  

[PIECEWISE([Typesetting:-mrow(Typesetting:-mi(
 

By using the coefficient = real; option, we include only points where the curve is locally smooth in real space: where it locally has a Puiseux series with real coefficients. This computation uses the underlying command RegularChainBranches, which returns a truncated Puiseux series at each of these points for all branches. 

rbs := RegularChainBranches(rc, R, [z]); -1 

map(print, rbs); -1 

 

 

[z = `*`(`^`(_T, 2)), y = `+`(`-`(`*`(`/`(1, 2), `*`(`^`(_T, 5), `*`(`+`(`*`(`^`(_T, 5)), `*`(2, `*`(RootOf(`+`(`*`(`^`(_Z, 2)), 1)))))))))), x = `+`(`*`(`/`(1, 8), `*`(`^`(_T, 2), `*`(`+`(`*`(`^`(_T,...
[z = `*`(`^`(_T, 2)), y = `+`(`-`(`*`(`/`(1, 2), `*`(`^`(_T, 5), `*`(`+`(`*`(`^`(_T, 5)), `-`(`*`(2, `*`(RootOf(`+`(`*`(`^`(_Z, 2)), 1))))))))))), x = `+`(`*`(`/`(1, 8), `*`(`^`(_T, 2), `*`(`+`(`*`(`^...
[z = _T, y = `+`(`*`(`^`(_T, 5)), 1), x = `+`(`-`(`*`(`^`(_T, 11))), `-`(`*`(2, `*`(`^`(_T, 6)))), `-`(_T))]
 

The first two results each contain a RootOf; function which represents I; (or `+`(`-`(I)); ). We see that they have nonreal coefficients. 

The last result is a truncated Puiseux series with only real coefficients. One might wonder if there are later terms of this series that have nonreal coefficients; but the algorithm used guarantees that this is not the case. The RegularChainBranches command can be instructed to return only truncated Puiseux series for real branches, as follows: 

real_branch := RegularChainBranches(rc, R, [z], coefficient = real);  

Typesetting:-mprintslash([real_branch := [[z = _T0, y = `+`(`*`(`^`(_T0, 5)), 1), x = `+`(`-`(`*`(`^`(_T0, 11))), `-`(`*`(2, `*`(`^`(_T0, 6)))), `-`(_T0))]]], [[[z = _T0, y = `+`(`*`(`^`(_T0, 5)), 1),...
 

When we examine it in the plot, we can see that this is locally a good approximation to the intersection curve for real values of the parameter, but the other two branches are not. Indeed, the z; -component of the curve is nonnegative, whereas the curve is present only for `<=`(z, 0); . When we re-parametrize these nonreal branches by the change of variables _T = `*`(I, `*`(T)); , we get the following. (The call to convert/radical converts the RootOf function into I; .) 

complex_branches := convert(eval(rbs[1 .. 2], _T = `*`(I, `*`(T))), radical); -1 

map(print, complex_branches); -1 

 

[z = `+`(`-`(`*`(`^`(T, 2)))), y = `+`(`-`(`*`(`+`(`*`(`/`(1, 2), `*`(I))), `*`(`^`(T, 5), `*`(`+`(`*`(I, `*`(`^`(T, 5))), `*`(2, `*`(I)))))))), x = `+`(`-`(`*`(`/`(1, 8), `*`(`^`(T, 2), `*`(`+`(`*`(`...
[z = `+`(`-`(`*`(`^`(T, 2)))), y = `+`(`-`(`*`(`+`(`*`(`/`(1, 2), `*`(I))), `*`(`^`(T, 5), `*`(`+`(`*`(I, `*`(`^`(T, 5))), `-`(`*`(2, `*`(I))))))))), x = `+`(`-`(`*`(`/`(1, 8), `*`(`^`(T, 2), `*`(`+`(...

Now the leading coefficient for each of the coordinates is real. Of course the (truncated) series itself can still not be embedded in real space, but we will be able to examine its real part, which should approximate the series well for small values of T; . Below, this real part of the truncated Puiseux series is shown in purple and the real truncated Puiseux series in pink. 

Plot_2d