New Features in Maple 16: Polynomial System Solving

Next Feature

Computing and manipulating the real solutions of a polynomial system is a requirement for many application areas, such as biological modeling, robotics, program verification, and control design, to name just a few. For example, an important problem in computational biology is to study the stability of the equilibria (or steady states) of biological systems. This question can often be reduced to solving a parametric system of polynomial equations and inequalities. 

The RegularChains package in Maple 16 provides a collection of tools for studying systems of polynomial equations, inequations and inequalities. It is particularly useful for solving and working with the real solutions of polynomial systems, such as the steady-state problem, amongst many others.  The new RegularChains features in Maple 16 include: 

  • Set theoretical operations for semi-algebraic sets (Complement,  Difference , Intersection, IsContained,  IsEmpty, Projection)
  • New solvers for popular types of systems (LinearSolve, RealComprehensiveTriangularize)
  • A new command, SuggestVariableOrder, for heuristically selecting a good variable order for computing a triangular decomposition of a polynomial system
  • Significant enhancements and performance improvements for the Triangularize,  RealTriangularize, and CylindricalAlgebraicDecompose commands

Polynomial System Solving

The two following applications highlight some of these new features.  

with(RegularChains); -1 

with(SemiAlgebraicSetTools); -1 

with(ParametricSystemTools); -1 


Application 1: Stability analysis of a parametric dynamical system
 

Consider a biological system described by the nonlinear multiple switch model 

S := Vector([diff(x(t), t) = `+`(`-`(x(t)), `/`(`*`(s), `*`(`+`(1, `*`(`^`(y(t), 2)))))), diff(y(t), t) = `+`(`-`(y(t)), `/`(`*`(s), `*`(`+`(1, `*`(`^`(x(t), 2))))))]) 

Vector[column](%id = 18446744078139347078)
 

where the unknowns x, y represent two protein concentrations and s represents the strength of unrepressed protein expression. This latter quantity is regarded as a time-constant parameter. 

The equilibria of S correspond to `and`(diff(x(t), t) = diff(y(t), t), diff(y(t), t) = 0), or equivalently, F = 0, where 

f := eval(`~`[rhs](S), [x(t) = x, y(t) = y]) 

Vector[column](%id = 18446744078140405750)
 

The following two Hurwitz determinants determine the stability of the hyperbolic equilibria of S: 

d[1] := `+`(`-`(diff(f[1], x)), `-`(diff(f[2], y))) 

2
 

d[2] := `+`(`*`(diff(f[1], x), `*`(diff(f[2], y))), `-`(`*`(diff(f[1], y), `*`(diff(f[2], x))))) 

`+`(1, `-`(`/`(`*`(4, `*`(`^`(s, 2), `*`(y, `*`(x)))), `*`(`^`(`+`(1, `*`(`^`(y, 2))), 2), `*`(`^`(`+`(1, `*`(`^`(x, 2))), 2))))))
 

The semi-algebraic system below encodes the asymptotically stable hyperbolic equilibria: 

P := [numer(f[1]) = 0, numer(f[2]) = 0, `>`(x, 0), `>`(y, 0), `>`(s, 0), `>`(numer(d[2]), 0)]; -1; `~`[print](P); -1 

        
`+`(`-`(x), `-`(`*`(x, `*`(`^`(y, 2)))), s) = 0

`+`(`-`(y), `-`(`*`(y, `*`(`^`(x, 2)))), s) = 0

`<`(0, x)

`<`(0, y)

`<`(0, s)

`<`(0, `+`(1, `*`(2, `*`(`^`(x, 2))), `*`(`^`(x, 4)), `*`(2, `*`(`^`(y, 2))), `*`(4, `*`(`^`(y, 2), `*`(`^`(x, 2)))), `*`(2, `*`(`^`(y, 2), `*`(`^`(x, 4)))), `*`(`^`(y, 4)), `*`(2, `*`(`^`(y, 4), `*`(...

 

We solve this problem by first computing a real comprehensive triangular decomposition of P, with respect to the parameter s, using the new command RealComprehensiveTriangularize:  

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

C := RealComprehensiveTriangularize(P, 1, R) 

[[[1, squarefree_semi_algebraic_system], [2, squarefree_semi_algebraic_system]], [[semi_algebraic_set, []], [semi_algebraic_set, [1]], [semi_algebraic_set, [2]]]]
[[[1, squarefree_semi_algebraic_system], [2, squarefree_semi_algebraic_system]], [[semi_algebraic_set, []], [semi_algebraic_set, [1]], [semi_algebraic_set, [2]]]]
 

This is a decomposition of the original system into several simpler triangular systems and some additional conditions on the parameter s: 

Display(C[1], R) 

 

Display(C[2], R) 

[[[`<`(s, -2), s = -2, And(`<`(-2, s), `<`(s, 0)), s = 0, s = 2], []], [[And(`<`(0, s), `<`(s, 2))], [1]], [[`<`(2, s)], [2]]]
 

Suppose we are interested in those values of s for which the biological system is bistable (that is, there are at least two stable equilibria). This means that we are looking for values of s such that P has at least 2 positive real solutions. We use the RealComprehensiveTriangularize command again and apply it to our previous result C: 

result := RealComprehensiveTriangularize(C, R, 2) 

[[[1, squarefree_semi_algebraic_system]], [[semi_algebraic_set, [1]]]]
 

The condition on s that we are looking for is given by the second entry: 

Display(result[2], R) 

[[[`<`(2, s)], [1]]]
 

The locations of the stable equilibria are described by the following triangular system from the first entry: 

Display(result[1], R) 

 

We now illustrate the result by discussing the special case s = 3. From the equations above, there are 2 stable equilibria given by: 

equilibria := [allvalues(solve(eval([`+`(`*`(x, `*`(y)), `-`(1)) = 0, `+`(`*`(`^`(x, 2)), `-`(`*`(3, `*`(x))), 1) = 0], s = 3)))] 

[{x = `+`(`/`(3, 2), `-`(`*`(`/`(1, 2), `*`(`^`(5, `/`(1, 2)))))), y = `+`(`/`(3, 2), `*`(`/`(1, 2), `*`(`^`(5, `/`(1, 2)))))}, {x = `+`(`/`(3, 2), `*`(`/`(1, 2), `*`(`^`(5, `/`(1, 2))))), y = `+`(`/`...
 

evalf(equilibria) 

[{x = .381966012, y = 2.618033988}, {x = 2.618033988, y = .381966012}]
 

We verify that the last inequality is satisfied: 

condition := eval(`>`(`+`(`*`(8, `*`(x, `*`(`^`(s, 3)))), `-`(`*`(6, `*`(x, `*`(`^`(s, 5))))), `-`(`*`(4, `*`(`^`(s, 2)))), `*`(5, `*`(`^`(s, 4))), `-`(`*`(`^`(s, 6))), `*`(x, `*`(`^`(s, 7)))), 0), s ... 

`<`(0, `+`(`*`(945, `*`(x)), `-`(360)))
 

evalf(eval(condition, equilibria[1])) 

`<`(0., .957881)
 

evalf(eval(condition, equilibria[2])) 

`<`(0., 2114.042119)
 

Below, we graphically display trajectories of the dynamical system in the special case. The first plot is a 2-D animation, and the second one is a 3-D static plot with time as the z-axis. We can see the two stable equilibria well in both plots, but there is also a third, unstable, equilibrium that can be seen best in the 3-D plot. 

special := eval(convert(S, list), s = 3) 

[diff(x(t), t) = `+`(`-`(x(t)), `/`(`*`(3), `*`(`+`(1, `*`(`^`(y(t), 2)))))), diff(y(t), t) = `+`(`-`(y(t)), `/`(`*`(3), `*`(`+`(1, `*`(`^`(x(t), 2))))))]
 

with(DEtools); -1 

IC[2] := [[x(0) = .5, y(0) = 0], seq([x(0) = i, y(0) = 0], i = 1 .. 5), [x(0) = 0, y(0) = .5], seq([x(0) = 0, y(0) = i], i = 1 .. 5), [x(0) = 4.5, y(0) = 5], seq([x(0) = i, y(0) = 5], i = 1 .. 4), [x(...
IC[2] := [[x(0) = .5, y(0) = 0], seq([x(0) = i, y(0) = 0], i = 1 .. 5), [x(0) = 0, y(0) = .5], seq([x(0) = 0, y(0) = i], i = 1 .. 5), [x(0) = 4.5, y(0) = 5], seq([x(0) = i, y(0) = 5], i = 1 .. 4), [x(...
IC[2] := [[x(0) = .5, y(0) = 0], seq([x(0) = i, y(0) = 0], i = 1 .. 5), [x(0) = 0, y(0) = .5], seq([x(0) = 0, y(0) = i], i = 1 .. 5), [x(0) = 4.5, y(0) = 5], seq([x(0) = i, y(0) = 5], i = 1 .. 4), [x(...
IC[2] := [[x(0) = .5, y(0) = 0], seq([x(0) = i, y(0) = 0], i = 1 .. 5), [x(0) = 0, y(0) = .5], seq([x(0) = 0, y(0) = i], i = 1 .. 5), [x(0) = 4.5, y(0) = 5], seq([x(0) = i, y(0) = 5], i = 1 .. 4), [x(...
IC[2] := [[x(0) = .5, y(0) = 0], seq([x(0) = i, y(0) = 0], i = 1 .. 5), [x(0) = 0, y(0) = .5], seq([x(0) = 0, y(0) = i], i = 1 .. 5), [x(0) = 4.5, y(0) = 5], seq([x(0) = i, y(0) = 5], i = 1 .. 4), [x(...
 

DEplot(special, [x(t), y(t)], t = 0 .. 20, x = 0 .. 5, y = 0 .. 5, IC[2], arrows = none, animate = true, numframes = 40, linecolor = sqrt(t))
DEplot(special, [x(t), y(t)], t = 0 .. 20, x = 0 .. 5, y = 0 .. 5, IC[2], arrows = none, animate = true, numframes = 40, linecolor = sqrt(t))
 

Plot_2d


IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
IC[3] := [seq([x(0) = .1, y(0) = `+`(`*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(1, `*`(0.4e-1, `*`(j)))], j = 0 .. 25), seq([x(0) = .1, y(0) = `+`(4, `*`(.1, `*`(j)))], j = 0 .. 1...
 

DEplot3d(special, {x(t), y(t)}, t = 0 .. 50, x = 0 .. 5, y = 0 .. 5, IC[3], scene = [x(t), y(t), t], thickness = 1, linecolor = sqrt(t), axes = frame)
DEplot3d(special, {x(t), y(t)}, t = 0 .. 50, x = 0 .. 5, y = 0 .. 5, IC[3], scene = [x(t), y(t), t], thickness = 1, linecolor = sqrt(t), axes = frame)
 

Plot_2d

 
 
Application 2: Verifying mathematical identities by branch cut analysis 

Let z be a complex variable. Which of the following two identities holds for all `in`(z, complex)? 

identity1 := sqrt(`+`(`*`(`^`(z, 2)), `-`(1))) = `*`(sqrt(`+`(z, `-`(1))), `*`(sqrt(`+`(z, 1)))) 

`*`(`^`(`+`(`*`(`^`(z, 2)), `-`(1)), `/`(1, 2))) = `*`(`^`(`+`(z, `-`(1)), `/`(1, 2)), `*`(`^`(`+`(z, 1), `/`(1, 2))))
 

identity2 := sqrt(`+`(1, `-`(`*`(`^`(z, 2))))) = `*`(sqrt(`+`(1, `-`(z))), `*`(sqrt(`+`(z, 1)))) 

`*`(`^`(`+`(1, `-`(`*`(`^`(z, 2)))), `/`(1, 2))) = `*`(`^`(`+`(1, `-`(z)), `/`(1, 2)), `*`(`^`(`+`(z, 1), `/`(1, 2))))
 

In Maple, the branch cut for the square root function sqrt(z) is the negative real axis `<`(z, 0). If we rewrite z = `+`(x, `*`(I, `*`(y))) in Cartesian coordinates, we can express that condition by the simple semi-algebraic system `and`(Im(z) = y, y = 0), `and`(Re(z) = x, `<`(x, 0)). For example, in identity1 above, the branch cuts of the three square root functions are, from left to right: 

 

assume(x, real) 

assume(y, real) 

bc1 := [Im(`+`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2)), `-`(1))) = 0, `<`(Re(`+`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2)), `-`(1))), 0)] 

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

bc2 := [Im(`+`(x, `*`(I, `*`(y)), `-`(1))) = 0, `<`(Re(`+`(x, `*`(I, `*`(y)), `-`(1))), 0)] 

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

bc3 := [Im(`+`(x, `*`(I, `*`(y)), 1)) = 0, `<`(Re(`+`(x, `*`(I, `*`(y)), 1)), 0)] 

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

The following plot depicts these three branch cuts. 

with(plottools); -1 

with(plots); -1 

branchcuts := display(line([0, -10], [0, 10], color = red), line([-1, -.12], [1, -.12], color = red), line([-10, 0.5e-1], [1, 0.5e-1], color = green), line([-10, -.12], [-1, -.12], color = blue), thic...
branchcuts := display(line([0, -10], [0, 10], color = red), line([-1, -.12], [1, -.12], color = red), line([-10, 0.5e-1], [1, 0.5e-1], color = green), line([-10, -.12], [-1, -.12], color = blue), thic...
branchcuts := display(line([0, -10], [0, 10], color = red), line([-1, -.12], [1, -.12], color = red), line([-10, 0.5e-1], [1, 0.5e-1], color = green), line([-10, -.12], [-1, -.12], color = blue), thic...
branchcuts := display(line([0, -10], [0, 10], color = red), line([-1, -.12], [1, -.12], color = red), line([-10, 0.5e-1], [1, 0.5e-1], color = green), line([-10, -.12], [-1, -.12], color = blue), thic...
 

branchcuts 

Plot_2d

By continuity arguments, it is now sufficient to check the proposed identity at finitely many points covering all possible combinations of branch cuts. Such a set of points can be computed using the command CylindricalAlgebraicDecompose. This example takes advantage of one of the newly supported output types for this command, added in Maple 16. 

b[1] := [bc1, bc2, bc3] 

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

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

S := CylindricalAlgebraicDecompose(b[1], R, output = rootof) 

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

The union of all branch cuts was decomposed into the 7 disjoint intervals above. We pick one sample point for each interval: 

S := [[x = -2, y = 0], S[2], [x = -`/`(1, 2), y = 0], [x = 0, y = -1], S[5], [x = 0, y = 1], [x = `/`(1, 2), y = 0]] 

[[x = -2, y = 0], [x = -1, y = 0], [x = -`/`(1, 2), y = 0], [x = 0, y = -1], [x = 0, y = 0], [x = 0, y = 1], [x = `/`(1, 2), y = 0]]
[[x = -2, y = 0], [x = -1, y = 0], [x = -`/`(1, 2), y = 0], [x = 0, y = -1], [x = 0, y = 0], [x = 0, y = 1], [x = `/`(1, 2), y = 0]]
 

We also add one sample point 1, 1 not contained on any branch cut: 

S := [S[], [x = 1, y = 1]] 

[[x = -2, y = 0], [x = -1, y = 0], [x = -`/`(1, 2), y = 0], [x = 0, y = -1], [x = 0, y = 0], [x = 0, y = 1], [x = `/`(1, 2), y = 0], [x = 1, y = 1]]
[[x = -2, y = 0], [x = -1, y = 0], [x = -`/`(1, 2), y = 0], [x = 0, y = -1], [x = 0, y = 0], [x = 0, y = 1], [x = `/`(1, 2), y = 0], [x = 1, y = 1]]
 

The following graph displays the sample points in relation to the branch cuts. 

display(branchcuts, seq(point(`~`[rhs](P)), `in`(P, S)), symbol = solidcircle, symbolsize = 20) 

Plot_2d

Now we check whether the identity  identity1 holds for all of these sample points: 

f[1] := eval((`+`(lhs, `-`(rhs)))(identity1), z = `+`(x, `*`(I, `*`(y)))) 

`+`(`*`(`^`(`+`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2)), `-`(1)), `/`(1, 2))), `-`(`*`(`^`(`+`(x, `*`(I, `*`(y)), `-`(1)), `/`(1, 2)), `*`(`^`(`+`(x, `*`(I, `*`(y)), 1), `/`(1, 2))))))
 

seq(radnormal(eval(f[1], P)), `in`(P, S)) 

`+`(`*`(2, `*`(`^`(3, `/`(1, 2))))), 0, 0, `*`(`*`(2, `*`(I)), `*`(`^`(2, `/`(1, 2)))), 0, 0, 0, 0
 

We conclude that identity1 is not an identity, and we can even say precisely when it fails: on the first interval `and`(`<`(x, -1), y = 0), and on the fourth interval `and`(x = 0, `<`(y, 0)). 

Now we proceed in the same way for the second proposed identity identity2. The branch cuts are: 

bc1 := [Im(`+`(1, `-`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2))))) = 0, `<`(Re(`+`(1, `-`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2))))), 0)] 

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

bc2 := [Im(`+`(1, `+`(`-`(x), `-`(`*`(`+`(I), `*`(y)))))) = 0, `<`(Re(`+`(1, `+`(`-`(x), `-`(`*`(`+`(I), `*`(y)))))), 0)] 

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

bc3 := [Im(`+`(x, `*`(I, `*`(y)), 1)) = 0, `<`(Re(`+`(x, `*`(I, `*`(y)), 1)), 0)] 

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

branchcuts := display(line([-10, 0], [-1, 0], color = red), line([1, 0], [10, 0], color = red), line([1, .12], [10, .12], color = green), line([-10, -.12], [-1, -.12], color = blue), thickness = 4, vi...
branchcuts := display(line([-10, 0], [-1, 0], color = red), line([1, 0], [10, 0], color = red), line([1, .12], [10, .12], color = green), line([-10, -.12], [-1, -.12], color = blue), thickness = 4, vi...
branchcuts := display(line([-10, 0], [-1, 0], color = red), line([1, 0], [10, 0], color = red), line([1, .12], [10, .12], color = green), line([-10, -.12], [-1, -.12], color = blue), thickness = 4, vi...
branchcuts := display(line([-10, 0], [-1, 0], color = red), line([1, 0], [10, 0], color = red), line([1, .12], [10, .12], color = green), line([-10, -.12], [-1, -.12], color = blue), thickness = 4, vi...
 

branchcuts 


Plot_2d

b[2] := [bc1, bc2, bc3] 

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

S := CylindricalAlgebraicDecompose(b[2], R, output = rootof) 

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

We need only 2+1 sample points in this case: 

S := [[x = -2, y = 0], [x = 0, y = 0], [x = 2, y = 0]] 

[[x = -2, y = 0], [x = 0, y = 0], [x = 2, y = 0]]
 

display(branchcuts, seq(point(`~`[rhs](P)), `in`(P, S)), symbol = solidcircle, symbolsize = 20) 

Plot_2d
 

f[2] := eval((`+`(lhs, `-`(rhs)))(identity2), z = `+`(x, `*`(I, `*`(y)))) 

`+`(`*`(`^`(`+`(1, `-`(`*`(`^`(`+`(x, `*`(I, `*`(y))), 2)))), `/`(1, 2))), `-`(`*`(`^`(`+`(1, `-`(x), `-`(`*`(`+`(I), `*`(y)))), `/`(1, 2)), `*`(`^`(`+`(x, `*`(I, `*`(y)), 1), `/`(1, 2))))))
 

seq(radnormal(eval(f[2], P)), `in`(P, S)) 

0, 0, 0
 

Thus we have proved that identity2 is indeed an identity. 

Next Feature

How to Proceed:   Pricing & Purchase Evaluate Upgrade Get Price Quote Buy Online