New Features in Maple 17: Advanced Mathematics
 

Next Feature

Maple 17 is the result of many thousands of hours of work by world-class mathematicians at Maplesoft as well as experts in research labs and universities all around the world. Maple 17 offers numerous advancements in a variety of branches of mathematics that push the frontiers of mathematical knowledge and Maple’s capabilities.

In addition to the work done in group theory, statistics, physics, differential geometry, and more, Maple 17 also provides:

  • Ground-breaking achievements in solving a whole new class of differential equations
  • Major advancements in solving systems of equations, introducing new solution methods for systems of linear equations and  inequalities and for systems involving nonlinear polynomial inequalities
  • New algorithms for finding limits of bivariate rational functions
  • Advancements in working with algebraic curves
  • Significant improvements in graph theory, including new algorithms and scalability improvements
  • More functionality for handling and exploring branch cuts

Details and Examples

Bivariate Limits

The limit command has been enhanced for the case of limits of bivariate rational functions. Many such limits that could not be determined previously are now computable. The new algorithm specifically handles the case where the function has an isolated singularity at the limit point.

In Maple 16, the following limit calls would return unevaluated, but they can be computed in Maple 17:

Example 1 

Example 2 

Example 3 

> f := `/`(`*`(x, `*`(y)), `*`(`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))))); -1
 

> limit(f, {x = 0, y = 0})
 

undefined

> g := `/`(`*`(x, `*`(`^`(y, 2))), `*`(`+`(`*`(`^`(x, 4)), `*`(`^`(y, 2))))); -1
 

> limit(g, {x = 0, y = 0})
 

0

> h := `/`(`*`(`+`(`*`(`^`(x, 4)), `*`(x, `*`(`^`(y, 2))), `*`(`^`(x, 2)), `*`(`^`(y, 2)))), `*`(`+`(`-`(`*`(`^`(x, 2), `*`(y))), `*`(`^`(x, 2)), `*`(`^`(y, 2))))); -1
 

> limit(h, {x = 0, y = 0})
 

1

Let us plot these three functions in the neighborhood of the origin:

> pf := plot3d(f, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; pf
pf := plot3d(f, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; pf
 

Plot_2d

> pg := plot3d(g, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; pg
pg := plot3d(g, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; pg
 

Plot_2d

> ph := plot3d(h, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; ph
ph := plot3d(h, x = -.1 .. .1, y = -.1 .. .1, axes = boxed); -1; ph
 

Plot_2d

In the last two examples, we can visually verify the existence of the limit.

In the first example, by inspecting the graph we can identify two different directions with different limits, namely y = x and y = `+`(`-`(x)). Indeed:

> eval(f, y = x)
 

`/`(1, 2)

> eval(f, y = `+`(`-`(x)))
 

-`/`(1, 2)

> with(plots); -1
 

> c1 := spacecurve([x, x, `/`(1, 2)], x = -.1 .. .1, color = red, thickness = 3); -1
 

> c2 := spacecurve([x, `+`(`-`(x)), -`/`(1, 2)], x = -.1 .. .1, color = blue, thickness = 3); -1
c2 := spacecurve([x, `+`(`-`(x)), -`/`(1, 2)], x = -.1 .. .1, color = blue, thickness = 3); -1
 

> display(pf, c1, c2)
 

Plot_2d

How does Maple determine these limits? Let us consider a circle C given by `+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))) = `*`(`^`(r, 2)), where we will later let proc (r) options operator, arrow; 0 end proc. Using the theory of Lagrange multipliers, the extremal values (maxima and minima) of the function f  on the circle, for a fixed radius r, satisfy the condition that the gradient of the function and the gradient of the constraint equation of the circle are parallel:

> C := `+`(`-`(`*`(`^`(r, 2))), `*`(`^`(x, 2)), `*`(`^`(y, 2))); -1
 

> with(VectorCalculus); -1
 

> df := `~`[normal](Jacobian([f], [x, y]))
 

df := Matrix(%id = 18446744078171316462)

> dC := Jacobian([C], [x, y])
 

dC := Matrix(%id = 18446744078171317302)

> eq := numer(normal(`+`(`-`(`*`(dC[1, 1], `*`(df[1, 2]))), `*`(dC[1, 2], `*`(df[1, 1])))))
 

`+`(`-`(`*`(2, `*`(`^`(x, 2)))), `*`(2, `*`(`^`(y, 2))))

Thus, the maximal and minimal values of f  on C occur when both C = 0 and eq = 0 For the bivariate limit, this means that it is sufficient to consider only those critical paths satisfying eq = 0. In the example, there are two such paths, namely, y = x and y = `+`(`-`(x)), and indeed they are exactly the same two paths from above leading to the extremal values `/`(1, 2) and -`/`(1, 2). Since we have found two directions with different limits, we can conclude that the bivariate limit does not exist.

The following graph depicts both critical paths in the x, y-plane:

> implicitplot(eq, x = -.1 .. .1, y = -.1 .. .1, numpoints = 10000, thickness = 2)
implicitplot(eq, x = -.1 .. .1, y = -.1 .. .1, numpoints = 10000, thickness = 2)
 

Plot_2d

Note that the critical paths are not necessarily always lines. For example, in the second example above:

> dg := `~`[normal](Jacobian([g], [x, y]))
 

dg := Matrix(%id = 18446744078171312726)

> eq := numer(normal(`+`(`-`(`*`(dC[1, 1], `*`(dg[1, 2]))), `*`(dC[1, 2], `*`(dg[1, 1])))))
 

`+`(`-`(`*`(2, `*`(y, `*`(`+`(`*`(2, `*`(`^`(x, 6))), `*`(3, `*`(`^`(x, 4), `*`(`^`(y, 2)))), `-`(`*`(`^`(y, 4)))))))))

> implicitplot(eq, x = -.1 .. .1, y = -.1 .. .1, numpoints = 100000, thickness = 2)
implicitplot(eq, x = -.1 .. .1, y = -.1 .. .1, numpoints = 100000, thickness = 2)
 

Plot_2d

In fact, we can give nonlinear closed form expressions for the critical paths:

> paths := [solve(eq, y)]
 

[0, `+`(`*`(`/`(1, 2), `*`(`^`(`+`(`*`(2, `*`(x, `*`(`^`(`+`(`*`(9, `*`(`^`(x, 2))), 8), `/`(1, 2))))), `*`(6, `*`(`^`(x, 2)))), `/`(1, 2)), `*`(x)))), `+`(`-`(`*`(`/`(1, 2), `*`(`^`(`+`(`*`(2, `*`(x,...

In our example, the limits along all of the critical paths are identical:

> eval_paths := map(proc (z) options operator, arrow; evala(Normal(eval(g, y = z))) end proc, paths)
 

[0, `+`(`/`(`*`(`/`(1, 4), `*`(`+`(`*`(x, `*`(`^`(`+`(`*`(9, `*`(`^`(x, 2))), 8), `/`(1, 2)))), `*`(3, `*`(`^`(x, 2))), `-`(4)), `*`(x))), `*`(`+`(`*`(2, `*`(`^`(x, 2))), `-`(1))))), `+`(`/`(`*`(`/`(1...

> `~`[limit](eval_path, x = 0)
 

[0, 0, 0, 0, 0]

This is a proof that the bivariate limit of g at the origin exists and is equal to 0.

The following graph depicts the function g as well as all of the critical paths:

> curves := map(proc (z) options operator, arrow; plots[spacecurve]([x, z, eval(g, y = z)], x = -.1 .. .1, color = red, thickness = 3) end proc, paths); -1
curves := map(proc (z) options operator, arrow; plots[spacecurve]([x, z, eval(g, y = z)], x = -.1 .. .1, color = red, thickness = 3) end proc, paths); -1
 

> display(pg, op(curves))
 

Plot_2d

Branch Cuts of Mathematical Expressions

The FunctionAdvisor can now compute, algebraically, branch cuts (discontinuities) of mathematical expressions involving compositions of +, *, ^, and possibly fractional powers, with every mathematical function call for which FunctionAdvisor knows the branch cuts, with no restrictions to the level of nesting in the expression. This new feature includes the ability to generate 2-D and 3-D plots of the cuts, so as to view where the cuts are located and how the function evolves from one cut to another one.

Example

> expression := ln(sqrt(`*`(`^`(z, 3))))
 

> FunctionAdvisor(branch_cuts, expression, plot = 3.)
 

Plot_2d Plot_2d

[`+`(`*`(`/`(1, 2), `*`(ln(`*`(`^`(z, 3)))))), `<`(z, 0), And(Re(z) = `+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2)), `*`(Im(z)))))), `<`(Im(z), 0)), And(Re(z) = `+`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))...

The segments of discontinuity are therefore located in the imaginary part of expression, Im(`+`(`*`(`/`(1, 2), `*`(ln(`*`(`^`(z, 3))))))).

You can rotate these plots with the mouse and use any of the options of plot, plot3d and plotcompare to observe the cuts in different ways.

A 2-D more abstract visualization of the same cuts:

> FunctionAdvisor(branch_cuts, expression, plot = `2D`)
 

Plot_2d
[`+`(`*`(`/`(1, 2), `*`(ln(`*`(`^`(z, 3)))))), `<`(z, 0), And(Re(z) = `+`(`-`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2)), `*`(Im(z)))))), `<`(Im(z), 0)), And(Re(z) = `+`(`*`(`/`(1, 3), `*`(`^`(3, `/`(1, 2))...

Example

A classic example in the theory of branch cut calculation: Kahan's teardrop

> expression := `+`(`*`(2, `*`(arccosh(`+`(1, `*`(`+`(`*`(2, `*`(z))), `/`(1, 3)))))), `-`(`*`(2, `*`(arccosh(`/`(`*`(`+`(12, `*`(5, `*`(z)))), `*`(`+`(`*`(3, `*`(z)), 12))))))), `-`(`*`(2, `*`(arccosh(...
 

This expression is equal to zero for all values of z in the complex plane except for a small teardrop shaped region over the negative real axis, as shown in the following 3-D plot on the right (you can rotate it with the mouse). To see just the plot, without displaying the complicated algebraic description of the cuts, end the input with colon ':'

> FunctionAdvisor(branch_cuts, expression, plot = 3., shift_range = -3, scale_range = 2, orientation = [-124, 20, 14])
 

Plot_2d Plot_2d

Differential Equations

A new class of ordinary differential equations of Abel type, with non-constant invariants and depending on two arbitrary parameters, is now fully solvable in Maple 17.  Additionally, a large number of changes were added to the routines for tackling ordinary and partial differential equations.

Ordinary differential equation

For 1st order ODEs, the simplest problem beyond the reach of complete solving algorithms is known as Abel equations. These are equations of the form

where `≡`(y, y(x)) is the unknown and the `≡`(f[i], f[i](x)) and `≡`(g[j], g[j](x)) are arbitrary functions of x. The biggest subclass of Abel equations known to be solvable was discovered by our research team and is the AIR 4-parameter class. New in Maple 17, another 2-parameter class of Abel equations, beyond the AIR class, is now known and solvable.

Example

>
 

> sol[1] := dsolve(ode[1])
 

Typesetting:-mprintslash([sol[1] := `+`(_C1, `/`(`*`(`^`(`+`(`-`(`/`(`*`(y(x)), `*`(beta))), `-`(1)), `+`(beta, 1)), `*`(`^`(`+`(`-`(`/`(`*`(y(x)), `*`(alpha))), 1), `+`(alpha, 1)))), `*`(`^`(`+`(`*`(...

> odetest(sol[1], ode[1])
 

0

The related class of Abel equations that is now entirely solvable consists of the set of equations that can be obtained from equation above by changing variables

proc (x) options operator, arrow; F(x) end proc, proc (y) options operator, arrow; `/`(`*`(`+`(`*`(P[1](x), `*`(y)), Q[1](x))), `*`(`+`(`*`(P[2](x), `*`(y)), Q[2](x)))) end proc

where and the four are arbitrary rational functions of x; this is the most general rational transformation that preserves the form of Abel equations and thus generates Abel ODE classes.

Several changes were performed in the existing algorithms to optimize their functioning, resulting in new solutions for problems that were out of reach in previous releases.

Example: Second order linear equations

> ode[2] := `+`(diff(y(x), x, x), `-`(`*`(`+`(l, `/`(`*`(`*`(2, `+`(`*`(`^`(a, 2)), `*`(`^`(b, 2))))), `*`(`^`(`+`(`*`(a, `*`(cos(x))), `*`(b, `*`(sin(x)))), 2)))), `*`(y(x))))) = 0; -1
 

> sol[2] := dsolve(ode[2], useInt)
 

Typesetting:-mprintslash([sol[2] := y(x) = `+`(`*`(_C1, `*`(`^`(`+`(`*`(`+`(`*`(`^`(b, 2), `*`(l)), `-`(`*`(`^`(a, 2)))), `*`(`^`(tan(x), 2))), `*`(2, `*`(a, `*`(b, `*`(`+`(l, 1), `*`(tan(x)))))), `*`...
Typesetting:-mprintslash([sol[2] := y(x) = `+`(`*`(_C1, `*`(`^`(`+`(`*`(`+`(`*`(`^`(b, 2), `*`(l)), `-`(`*`(`^`(a, 2)))), `*`(`^`(tan(x), 2))), `*`(2, `*`(a, `*`(b, `*`(`+`(l, 1), `*`(tan(x)))))), `*`...
Typesetting:-mprintslash([sol[2] := y(x) = `+`(`*`(_C1, `*`(`^`(`+`(`*`(`+`(`*`(`^`(b, 2), `*`(l)), `-`(`*`(`^`(a, 2)))), `*`(`^`(tan(x), 2))), `*`(2, `*`(a, `*`(b, `*`(`+`(l, 1), `*`(tan(x)))))), `*`...
Typesetting:-mprintslash([sol[2] := y(x) = `+`(`*`(_C1, `*`(`^`(`+`(`*`(`+`(`*`(`^`(b, 2), `*`(l)), `-`(`*`(`^`(a, 2)))), `*`(`^`(tan(x), 2))), `*`(2, `*`(a, `*`(b, `*`(`+`(l, 1), `*`(tan(x)))))), `*`...
Typesetting:-mprintslash([sol[2] := y(x) = `+`(`*`(_C1, `*`(`^`(`+`(`*`(`+`(`*`(`^`(b, 2), `*`(l)), `-`(`*`(`^`(a, 2)))), `*`(`^`(tan(x), 2))), `*`(2, `*`(a, `*`(b, `*`(`+`(l, 1), `*`(tan(x)))))), `*`...

> odetest(sol[2], ode[2])
 

0

Example: Nonlinear ODE problems with initial values

> ode[3] := {diff(y(x), x) = `+`(`-`(`*`(`+`(1, `-`(y(x))), `*`(sqrt(abs(y(x))))))), y(0) = `/`(1, 2)}; -1
 

> sol[3] := dsolve(ode[3], implicit)
 

Typesetting:-mprintslash([sol[3] := `+`(x, `-`(PIECEWISE([`+`(`*`(2, `*`(arctan(`*`(`^`(`+`(`-`(y(x))), `/`(1, 2))))))), `<=`(y(x), 0)], [`+`(`-`(`*`(2, `*`(arctanh(`*`(`^`(y(x), `/`(1, 2)))))))), `<`...

> odetest(sol[3], ode[3])
 

{0}

Example: Third order nonlinear ODE depending on arbitrary function, a reduction to a Riccati equation

>
 

> sol[4] := dsolve(ode[4], y(x))
 


> odetest(sol[4], ode[4])
 

0

Partial differential equations

Changes were done in the existing algorithms for linear and nonlinear PDEs and systems of them, resulting in new solutions for more PDE problems.

Examples: Linear PDEs

>
 

> sol[1] := pdsolve(pde[1])
 

Results can be verified with pdetest

>
 

0

> pde[2] := `+`(`*`(x, `*`(diff(u(x, y), x, x))), `-`(`*`(2, `*`(diff(u(x, y), x)))), `-`(`*`(x, `*`(diff(u(x, y), y, y))))) = 0; -1
 

> sol[2] := pdsolve(pde[2])
 

Typesetting:-mprintslash([sol[2] := u(x, y) = `+`(`-`(`*`((D(_F1))(`+`(`-`(x), y)), `*`(x))), `-`(_F1(`+`(`-`(x), y))), `*`((D(_F2))(`+`(x, y)), `*`(x)), `-`(_F2(`+`(x, y))))], [u(x, y) = `+`(`-`(`*`(...

> pdetest(sol[2], pde[2])
 

0

> pde[3] := `+`(`*`(2, `*`(diff(u(x, y), y, x))), diff(u(x, y), y, y)) = f(x, y); -1
 

> sol[3] := pdsolve(pde[3])
 

Typesetting:-mprintslash([sol[3] := u(x, y) = `+`(_F1(x), _F2(`+`(`*`(2, `*`(y)), `-`(x))), `*`(`/`(1, 2), `*`(Int(Int(f(_b, _a), _a = `` .. `+`(`*`(`/`(1, 2), `*`(_b)), `-`(`*`(`/`(1, 2), `*`(x))), y...

> pdetest(sol[3], pde[3])
 

0

Examples: PDE problems with boundary conditions

> pde[4] := diff(u(x, t), t) = `+`(`*`(5, `*`(diff(u(x, t), x, x)))); -1
 

> bc[4] := u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = `+`(`*`(3, `*`(sin(`+`(`*`(2, `*`(x)))))), `*`(6, `*`(sin(`+`(`*`(5, `*`(x)))))))
 

Typesetting:-mprintslash([bc[4] := u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = `+`(`*`(3, `*`(sin(`+`(`*`(2, `*`(x)))))), `*`(6, `*`(sin(`+`(`*`(5, `*`(x)))))))], [u(0, t) = 0, u(Pi, t) = 0, u(x, 0) = `+`(`*...

> sol[4] := pdsolve([pde[4], bc[4]])
 

Typesetting:-mprintslash([sol[4] := u(x, t) = `+`(`*`(3, `*`(exp(`+`(`-`(`*`(20, `*`(t))))), `*`(sin(`+`(`*`(2, `*`(x))))))), `*`(6, `*`(exp(`+`(`-`(`*`(125, `*`(t))))), `*`(sin(`+`(`*`(5, `*`(x))))))...

The verification with pdetest also verifies that the solution satisfies the boundary conditions given

> pdetest(sol[4], [pde[4], bc[4]])
 

[0, 0, 0, 0]

Nonlinear PDE systems

> pde[5] := {`+`(diff(u(x, y), x, x), diff(u(x, y), y, y)) = 0, `+`(`*`(`^`(diff(u(x, y), x), 2)), `*`(`^`(diff(u(x, y), y), 2)), `-`(exp(`+`(`*`(2, `*`(x)))))) = 0}; -1
 

> sol[5] := pdsolve(pde[5])
 

Typesetting:-mprintslash([sol[5] := {u(x, y) = `+`(`/`(`*`(`/`(1, 4), `*`(`+`(exp(`+`(`-`(`*`(`+`(I), `*`(y))), x)), `*`(4, `*`(`^`(_C2, 2), `*`(exp(`+`(x, `*`(I, `*`(y))))))), `*`(`+`(`*`(4, `*`(_C1)...

> pdetest(sol[5], pde[5])
 

{0}

Graph Theory

Maple 17 has significant improvements to several GraphTheory commands, including:

In addition, two new commands were added to the package:

Improvements to TuttePolynomial, ChromaticPolynomial, FlowPolynomial, AcyclicPolynomial and RankPolynomial commands

The TuttePolynomial command uses a new algorithm that is faster and consumes less memory.

For example, computation of the Tutte polynomial of the 12 vertex complete graph (CompleteGraph(12)) now requires 1/3 the memory, and completes in 1/3 the time.

For the 9,3 generalized Petersen graph (PetersenGraph(9,3)) the memory usage has dropped to 1/6, and the time to 1/8 of the time needed in Maple 16.

> with(GraphTheory); -1
 

> with(SpecialGraphs); -1
 

Example

> G := CompleteGraph(12); -1
 

> CodeTools:-Usage(TuttePolynomial(G, 'x', 'y')); -1
 

memory used=106.98MiB, alloc change=0 bytes, cpu time=2.67s, real time=2.83s

> G := GeneralizedPetersenGraph(9, 3); -1
 

> CodeTools:-Usage(TuttePolynomial(G, 'x', 'y')); -1
 

memory used=3.21MiB, alloc change=0 bytes, cpu time=125.00ms, real time=149.00ms

These improvements propagate to the commands ChromaticPolynomial, FlowPolynomial, AcyclicPolynomial and RankPolynomial as well.

Improvements to IsIsomorphic command

The IsIsomorphic command now utilizes a new algorithm that scales better to larger problems. As an example, two isomorphic soccer ball graphs (SpecialGraphs[SoccerBallGraph]) with 60 vertices and 90 edges can be tested for isomorphism in less than one second in Maple 17, while in Maple 16 and earlier, this same computation did not complete in under 10 minutes.

Example

> S := SoccerBallGraph(); -1
 

> G := IsomorphicCopy(S); -1
 

> H := IsomorphicCopy(S); -1
 

> CodeTools:-Usage(IsIsomorphic(G, H, 'phi')); -1
 

memory used=448.29KiB, alloc change=0 bytes, cpu time=47.00ms, real time=50.00ms

The LaplacianMatrix command

The new LaplacianMatrix command can be used to compute the Laplacian matrix representation of a graph:

Example

> G := CompleteGraph(5); -1
 

> LaplacianMatrix(G)
 

Matrix(%id = 18446744078171284654)

The ReliabilityPolynomial command

The new ReliabilityPolynomial command can be used to compute the reliability polynomial of a graph:

Example

> ReliabilityPolynomial(G, 'x')
 

`*`(`+`(`*`(24, `*`(`^`(x, 6))), `*`(36, `*`(`^`(x, 5))), `*`(30, `*`(`^`(x, 4))), `*`(20, `*`(`^`(x, 3))), `*`(10, `*`(`^`(x, 2))), `*`(4, `*`(x)), 1), `*`(`^`(`+`(1, `-`(x)), 4)))

Intersecting Plane Curves

For Maple17, a new command, intersectcurves, has been added to the algcurves package.

> with(algcurves); -1
 

Example

Given two plane curves expressed as polynomials in two variables:

> c1 := `+`(`-`(`*`(`^`(x, 3))), `*`(2, `*`(`^`(y, 2))), x); -1
 

> c2 := `+`(`*`(`^`(x, 2)), `-`(`*`(`^`(y, 2))), `-`(1)); -1
 

The intersectcurves command computes the points of intersection along with their multiplicities. Each set of intersection points is expressed as a list containing three values. The first two are irreducible polynomials, the second of which is in only one of the variables. The third value is either 0 or 1. When 1, the point is affine; when 0 it is a point in the projective plane (at infinity).

> intersectcurves(c1, c2, x, y)
 

[[2, [`+`(1, x), y, 1]], [2, [`+`(`-`(1), x), y, 1]], [1, [`+`(`-`(2), x), `+`(`*`(`^`(y, 2)), `-`(3)), 1]]]

This is to be interpreted as follows. There are two intersection points of multiplicity 2, namely, (-1,0) and (1,0), indicating that they are points where the curves intersect tangentially. In addition, there are two intersection points of multiplicity 1, namely, 2, sqrt(3) and 2, `+`(`-`(sqrt(3))). We can see all intersection points graphically with a simple implicit plot.

> with(plots); -1
 

> implicitplot([c1, c2], x = -3 .. 3, y = -3 .. 3, gridrefine = 3, scaling = constrained, color = [
implicitplot([c1, c2], x = -3 .. 3, y = -3 .. 3, gridrefine = 3, scaling = constrained, color = [
 

Plot_2d

The following example illustrates the case when there are intersection points at infinity.

> d1 := `+`(`*`(`^`(x, 2), `*`(y)), `-`(y), `-`(1)); -1
 

> d2 := `+`(`*`(`^`(x, 2)), `-`(x)); -1
 

> intersectcurves(d1, d2, x, y)
 

[[1, [x, `+`(1, y), 1]], [5, [x, 1, 0]]]

> implicitplot([d1, d2], x = -3 .. 3, y = -3 .. 3, gridrefine = 3, scaling = constrained, color = [
implicitplot([d1, d2], x = -3 .. 3, y = -3 .. 3, gridrefine = 3, scaling = constrained, color = [
 

Plot_2d

There is one affine intersection point (0,-1) of multiplicity 1, plus one intersection point at infinity, of multiplicity 5. The tangent at this point is the line x = 0. Each of the two blue vertical lines intersects each of the two asymptotes of the red curve at the infinite point All but one of these intersections have multiplicity 1, and the exceptional one, between the blue line and the asymptote to the red curve at x = 1, comes from a tangential intersection, of multiplicity 2.

We can perform a change of coordinates in order to plot this in the affine plane. First, we homogenize the equations, introducing a third variable z, and then we substitute y = 1 in order to move the intersection point that was previously at infinity to the origin.

> h1 := eval(`+`(`*`(`^`(x, 2), `*`(y)), `-`(`*`(y, `*`(`^`(z, 2)))), `-`(`*`(`^`(z, 3)))), y = 1)
 

`+`(`-`(`*`(`^`(z, 3))), `*`(`^`(x, 2)), `-`(`*`(`^`(z, 2))))

> h2 := eval(`+`(`*`(`^`(x, 2)), `-`(`*`(x, `*`(z)))), y = 1)
 

`+`(`*`(`^`(x, 2)), `-`(`*`(x, `*`(z))))

> implicitplot([h1, h2], x = -1.5 .. 1.5, z = -1.5 .. 1.5, gridrefine = 3, scaling = constrained, color = [
implicitplot([h1, h2], x = -1.5 .. 1.5, z = -1.5 .. 1.5, gridrefine = 3, scaling = constrained, color = [
 

> intersectcurves(h1, h2, x, z)
 

[[5, [x, z, 1]], [1, [x, `+`(z, 1), 1]]]

After the coordinate transformation, all intersection points (of h1 and h2) are visible in the plot and there are no intersection points at infinity.

Plot_2d

Linear Inequality Solving

The LinearSolve command in the RegularChains package, whose purpose is to solve systems of linear equations and inequalities, has been enhanced in several ways:

  • A more efficient algorithm based on Fourier-Motzkin elimination, has been implemented.
  • Two new options, canonical and projection, for controlling the form of the output have been added.

The examples below illustrate the use of this command and the two new options for computations with polyhedra.

Change of representation for polyhedra

There are at least two ways to mathematically describe a polyhedron, either by specifying its corners, or by specifying its bounding hyperplanes. We start by illustrating this using a cube of side length 1 with one corner at the origin and its three adjacent edges parallel to the three coordinate axes.

Example

Plot the cube:

> with(plots); -1with(plottools); -1
 

> display(cuboid([0, 0, 0], [1, 1, 1]), axes = boxed, labels = [x, y, z])
 

Plot_2d

Define this cube by bounding planes:

> c1 := [`<=`(0, x), `<=`(x, 1), `<=`(0, y), `<=`(y, 1), `<=`(0, z), `<=`(z, 1)]
 

[`<=`(0, x), `<=`(x, 1), `<=`(0, y), `<=`(y, 1), `<=`(0, z), `<=`(z, 1)]

Alternatively, define the vertices of the same cube:

> V := Matrix([[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]])
 

V := Matrix(%id = 18446744078171276942)

Suppose now that we are given the vertex representation and want to automatically derive the bounding planes representation. The cube is defined as the convex hull of all its vertices:

> op(convert(`~`[`=`](Vector[row]([x, y, z]), add(`*`(lambda[i], `*`(V[i])), i = 1 .. 8)), list))
 

x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8])

> add(lambda[i], i = 1 .. 8) = 1
 

`+`(lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6], lambda[7], lambda[8]) = 1

> seq(`<=`(0, lambda[i]), i = 1 .. 8)
 

`<=`(0, lambda[1]), `<=`(0, lambda[2]), `<=`(0, lambda[3]), `<=`(0, lambda[4]), `<=`(0, lambda[5]), `<=`(0, lambda[6]), `<=`(0, lambda[7]), `<=`(0, lambda[8])

> seq(`<=`(lambda[i], 1), i = 1 .. 8)
 

`<=`(lambda[1], 1), `<=`(lambda[2], 1), `<=`(lambda[3], 1), `<=`(lambda[4], 1), `<=`(lambda[5], 1), `<=`(lambda[6], 1), `<=`(lambda[7], 1), `<=`(lambda[8], 1)

> sys := [x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8]), `+`(lambda[1], lambda[2], lambda[...
sys := [x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8]), `+`(lambda[1], lambda[2], lambda[...
sys := [x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8]), `+`(lambda[1], lambda[2], lambda[...
 

[x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8]), `+`(lambda[1], lambda[2], lambda[3], lam...
[x = `+`(lambda[5], lambda[6], lambda[7], lambda[8]), y = `+`(lambda[3], lambda[4], lambda[7], lambda[8]), z = `+`(lambda[2], lambda[4], lambda[6], lambda[8]), `+`(lambda[1], lambda[2], lambda[3], lam...

We can now use the LinearSolve command to eliminate all lambda[i] from this system of equations and inequalities. We use the projection option to specify that elimination is required, and the canonical option to eliminate redundant equations and inequalities.

> with(RegularChains); -1with(SemiAlgebraicSetTools); -1
 

> vars := [seq(lambda[i], i = 1 .. 8), z, y, x]
 

[lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6], lambda[7], lambda[8], z, y, x]

> LinearSolve(sys, PolynomialRing(vars), canonical, projection = 3)
 

[`<=`(0, x), `<=`(x, 1), `<=`(0, y), `<=`(y, 1), `<=`(0, z), `<=`(z, 1)]

This agrees with the representation c1 above.

Dual polyhedra

Example

For every convex 3-D polyhedron, there is a dual polyhedron with the roles of corners and faces swapped. Let us compute the dual polyhedron of the cube. We start with a representation of the cube with the origin as the center.

> c2 := [`<=`(x, 1), `<=`(`+`(`-`(x)), 1), `<=`(y, 1), `<=`(`+`(`-`(y)), 1), `<=`(z, 1), `<=`(`+`(`-`(z)), 1)]
 

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

Now we rewrite this in matrix form.

> B := Matrix([[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]])
 

B := Matrix(%id = 18446744078171273326)

> v := Vector([x, y, z])
 

v := Vector[column](%id = 18446744078171273806)

> c := Vector([1, 1, 1, 1, 1, 1])
 

c := Vector[column](%id = 18446744078171274286)

> `~`[`<=`](Typesetting:-delayDotProduct(B, v), c)
 

Vector[column](%id = 18446744078171267070)

In general, if a bounded convex polyhedron is given by a set of inequalities of the form `<=`(Typesetting:-delayDotProduct(B, v), 1), then the vertices of its dual are exactly the columns of the matrix B. The convex hull of these six vertices defines an octahedron. Using the same method as before, we compute the bounding planes for this octahedron.

> bound1 := op(convert(`~`[`=`](Vector[row]([x, y, z]), add(`*`(lambda[i], `*`(B[i])), i = 1 .. 6)), list))
 

x = `+`(lambda[1], `-`(lambda[2])), y = `+`(lambda[3], `-`(lambda[4])), z = `+`(lambda[5], `-`(lambda[6]))

> bound2 := add(lambda[i], i = 1 .. 6) = 1
 

`+`(lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6]) = 1

> bound3 := seq(`<=`(0, lambda[i]), i = 1 .. 6)
 

`<=`(0, lambda[1]), `<=`(0, lambda[2]), `<=`(0, lambda[3]), `<=`(0, lambda[4]), `<=`(0, lambda[5]), `<=`(0, lambda[6])

> bound4 := seq(`<=`(lambda[i], 1), i = 1 .. 6)
 

`<=`(lambda[1], 1), `<=`(lambda[2], 1), `<=`(lambda[3], 1), `<=`(lambda[4], 1), `<=`(lambda[5], 1), `<=`(lambda[6], 1)

> sys := [bound1, bound2, bound3, bound4]
 

[x = `+`(lambda[1], `-`(lambda[2])), y = `+`(lambda[3], `-`(lambda[4])), z = `+`(lambda[5], `-`(lambda[6])), `+`(lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6]) = 1, `<=`(0, lambda[1...
[x = `+`(lambda[1], `-`(lambda[2])), y = `+`(lambda[3], `-`(lambda[4])), z = `+`(lambda[5], `-`(lambda[6])), `+`(lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6]) = 1, `<=`(0, lambda[1...

> vars := [seq(lambda[i], i = 1 .. 6), z, y, x]
 

[lambda[1], lambda[2], lambda[3], lambda[4], lambda[5], lambda[6], z, y, x]

> S := LinearSolve(sys, PolynomialRing(vars), canonical, projection = 3)
 

[`<=`(-1, x), `<=`(x, 1), `<=`(`+`(`-`(x), `-`(1)), y), `<=`(`+`(`-`(1), x), y), `<=`(y, `+`(1, x)), `<=`(y, `+`(1, `-`(x))), `<=`(`+`(`-`(y), `-`(x), `-`(1)), z), `<=`(`+`(`-`(y), x, `-`(1)), z), `<=...
[`<=`(-1, x), `<=`(x, 1), `<=`(`+`(`-`(x), `-`(1)), y), `<=`(`+`(`-`(1), x), y), `<=`(y, `+`(1, x)), `<=`(y, `+`(1, `-`(x))), `<=`(`+`(`-`(y), `-`(x), `-`(1)), z), `<=`(`+`(`-`(y), x, `-`(1)), z), `<=...

The first two inequalities represent that projection of the octahedron onto the x axis:

> S[1 .. 2]
 

[`<=`(-1, x), `<=`(x, 1)]

Together with the following four inequalities, which involve only x and y but not z, this defines the projection of the octahedron onto the x, y-plane:

> S[1 .. 2]; 1S[3 .. 6]
 

[`<=`(-1, x), `<=`(x, 1)]
[`<=`(`+`(`-`(x), `-`(1)), y), `<=`(`+`(`-`(1), x), y), `<=`(y, `+`(1, x)), `<=`(y, `+`(1, `-`(x)))]

The following plot shows this projection (blue) together with the four boundary lines.

> inequal(S[3 .. 6], x = -2 .. 2, y = -2 .. 2)
 

Plot_2d

Finally, the last 8 inequalities, i.e., the ones containing z, define the bounding planes of the octahedron.

> S[-8 .. -1]
 

[`<=`(`+`(`-`(y), `-`(x), `-`(1)), z), `<=`(`+`(`-`(y), x, `-`(1)), z), `<=`(`+`(`-`(1), `-`(x), y), z), `<=`(`+`(`-`(1), x, y), z), `<=`(z, `+`(y, x, 1)), `<=`(z, `+`(y, `-`(x), 1)), `<=`(z, `+`(1, x...
[`<=`(`+`(`-`(y), `-`(x), `-`(1)), z), `<=`(`+`(`-`(y), x, `-`(1)), z), `<=`(`+`(`-`(1), `-`(x), y), z), `<=`(`+`(`-`(1), x, y), z), `<=`(z, `+`(y, x, 1)), `<=`(z, `+`(y, `-`(x), 1)), `<=`(z, `+`(1, x...
[`<=`(`+`(`-`(y), `-`(x), `-`(1)), z), `<=`(`+`(`-`(y), x, `-`(1)), z), `<=`(`+`(`-`(1), `-`(x), y), z), `<=`(`+`(`-`(1), x, y), z), `<=`(z, `+`(y, x, 1)), `<=`(z, `+`(y, `-`(x), 1)), `<=`(z, `+`(1, x...

The following 3-D plot displays the tetrahedron as well as the last of these 8 bounding planes. Rotate the plot with the mouse in order to see that the plane actually touches one of the faces.

> plane := Student:-LinearAlgebra:-PlanePlot((`+`(lhs, `-`(rhs)))(S[-1]) = 0, [x, y, z], caption =
plane := Student:-LinearAlgebra:-PlanePlot((`+`(lhs, `-`(rhs)))(S[-1]) = 0, [x, y, z], caption =
plane := Student:-LinearAlgebra:-PlanePlot((`+`(lhs, `-`(rhs)))(S[-1]) = 0, [x, y, z], caption =
 

> display(octahedron([0, 0, 0]), plane, axes = boxed, labels = [x, y, z])
 

Plot_2d

Semialgebraic System Solving

The SolveTools[SemiAlgebraic] command has been integrated directly into the solve command, such that many systems involving non-linear polynomial inequalities that could not be solved previously, are solved.

In Maple 16, no solutions were found for the following system, but in Maple 17 it is easily solved:

Example

> sol := solve({`<=`(0, `+`(`-`(`*`(`^`(y, 2))), x)), `<`(0, `+`(`-`(y), `*`(3, `*`(x)), `-`(2))), `<`(`+`(`*`(`^`(x, 2)), `*`(`^`(y, 2))), 9)}, [x, y])
 

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

A solution to such a system of inequalities is a decomposition of the feasible region of the inequalities into bands as can be seen in the inequality plot:

> plots[inequal](sol, x = 0 .. 3, y = -2 .. 2)
 

Plot_2d

Additionally, the SemiAlgebraic command now can build case discussions for systems with real-valued parameters.

> SolveTools[SemiAlgebraic]({`<`(a, x), `<`(`*`(a, `*`(`^`(x, 2))), b)}, [x])
 

piecewise(And(a = 0, `<`(0, b)), [[`<`(0, x)]], And(`<`(0, a), `<`(`*`(`^`(a, 3)), b)), [[`<`(a, x), `<`(x, `/`(`*`(sqrt(`*`(a, `*`(b)))), `*`(a)))]], And(`<`(a, 0), b = 0), [[`<`(0, x)], [`<`(a, x), ...

Solving Equations with Branch Cuts

New functionality has been added to the solve command to account for branch cuts in the input equations and build piecewise expressions that are correct for substitution of parameters. This functionality is controlled with the new option symbolic = false;. The default behavior in Maple 17 and previous versions of Maple, is the same as specifying symbolic = true;.

Example

> expr := a = `+`(sqrt(`+`(a, y)), 1); -1
 

> sol1 := solve(expr, [y], symbolic = true)
 

[[y = `+`(`*`(`^`(a, 2)), `-`(`*`(3, `*`(a))), 1)]]

> sol2 := solve(expr, [y], symbolic = false)
 

piecewise(And(`<=`(`+`(`*`(2, `*`(argument(`+`(a, `-`(1)))))), Pi), `<`(`+`(`-`(Pi)), `+`(`*`(2, `*`(argument(`+`(a, `-`(1)))))))), [[y = `+`(`*`(`^`(a, 2)), `-`(`*`(3, `*`(a))), 1)]], [])

While both sol1 and sol2 are valid for some values of the parameter a:

> eval(sol1, a = 1), eval(sol2, a = 1), solve(eval(expr, a = 1), [y])
 

[[y = -1]], [[y = -1]], [[y = -1]]

Only the solution using the symbolic = false option are valid for other values of the parameter:

> eval(sol1, a = -1), eval(sol2, a = -1), solve(eval(expr, a = -1), [y])
 

[[y = 5]], [], []