@Carl Love If you give it all the equations at once, it computes a Groebner basis for the side relations using a lexicographical ordering.  The basis captures all of the algebraic consequences of the relations.  In this case the computation becomes:

F := [phi-s[1], -eta[2]*w[2]+eta[3]*w[1]+eta[3]*w[2]-mu-eta[3]-s[2],
 2*mu*eta[2]*w[2]-2*mu*eta[3]*w[1]-2*mu*eta[3]*w[2]+eta[2]^2*w[2]
 -eta[3]^2*w[1]-eta[3]^2*w[2]+mu^2+2*mu*eta[3]+eta[3]^2-s[3],
 3*mu^2*eta[2]*w[2]-3*mu^2*eta[3]*w[1]-3*mu^2*eta[3]*w[2]
 +3*mu*eta[2]^2*w[2]-3*mu*eta[3]^2*w[1]-3*mu*eta[3]^2*w[2]
 +eta[2]^3*w[2]-eta[3]^3*w[1]-eta[3]^3*w[2]+mu^3+3*mu^2*eta[3]
 +3*mu*eta[3]^2+eta[3]^3-s[4], 4*mu^3*eta[2]*w[2]-4*mu^3*eta[3]*w[1]
 -4*mu^3*eta[3]*w[2]+6*mu^2*eta[2]^2*w[2]-6*mu^2*eta[3]^2*w[1]
 -6*mu^2*eta[3]^2*w[2]+4*mu*eta[2]^3*w[2]-4*mu*eta[3]^3*w[1]
 -4*mu*eta[3]^3*w[2]+eta[2]^4*w[2]-eta[3]^4*w[1]-eta[3]^4*w[2]
 +mu^4+4*mu^3*eta[3]+6*mu^2*eta[3]^2+4*mu*eta[3]^3+eta[3]^4-s[5]];

tord := plex(phi,w[1],w[2],mu,eta[2],eta[3],s[5],s[4],s[3],s[2],s[1]);

infolevel[Groebner] := 5:

# G := CodeTools:-Usage(Groebner:-Basis(F, tord));  # takes a long time

Maple's default strategy for lexicographical bases is to compute a total degree basis and convert, because direct elimination is often infeasible.  But in this case, direct elimination is much faster because the polynomials are close enough to a triangular form, and the total degree basis is large.  It took me about 2 minutes to compute directly.

G := CodeTools:-Usage(Groebner:-Basis(F, tord, method=direct));

It would be nice if Maple could detect this situation, but how to do that is unclear.  If you have computed G via the direct method, you can give it to simplify/siderels, i.e.:

result := simplify(kappa, G, tord):

I got:

[s[1]*s[2]+1, -s[2]*s[1], s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^2*s[2]-s[1]^2*s[3], -s[1]^2*s[3]-s[1]*s[2], 

s[3]*s[1]^2, 

s[1]^3*s[2]+2*s[1]^3*s[3]-s[1]^3*s[4]+s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^3*s[2]-2*s[1]^3*s[3]+s[1]^3*s[4], 

-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[2]-s[1]^2*s[3], s[1]^3*s[3]-s[1]^3*s[4], 

-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[3]-s[1]*s[2], s[1]^3*s[3]-s[1]^3*s[4], 

-s[1]^3*s[4]+s[1]^2*s[3], s[1]^3*s[4], 

s[1]^4*s[2]+3*s[1]^4*s[3]-3*s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[2]

+2*s[1]^3*s[3]-s[1]^3*s[4]+s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^4*s[2]-3*s[1]^4*s[3]+3*s[1]^4*s[4]-s[1]^4*s[5], 

-s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[2]-2*s[1]^3*s[3]+s[1]^3*s[4], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], 

-s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[2]-s[1]^2*s[3], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[3]-s[1]^3*s[4], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[3]+s[1]^3*s[4]

-s[1]^2*s[3]-s[1]*s[2], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[3]-s[1]^3*s[4], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]-s[1]^3*s[4]+s[1]^2*s[3], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[5]+s[1]^3*s[4], s[5]*s[1]^4]

@Carl Love If you give it all the equations at once, it computes a Groebner basis for the side relations using a lexicographical ordering.  The basis captures all of the algebraic consequences of the relations.  In this case the computation becomes:

F := [phi-s[1], -eta[2]*w[2]+eta[3]*w[1]+eta[3]*w[2]-mu-eta[3]-s[2],
 2*mu*eta[2]*w[2]-2*mu*eta[3]*w[1]-2*mu*eta[3]*w[2]+eta[2]^2*w[2]
 -eta[3]^2*w[1]-eta[3]^2*w[2]+mu^2+2*mu*eta[3]+eta[3]^2-s[3],
 3*mu^2*eta[2]*w[2]-3*mu^2*eta[3]*w[1]-3*mu^2*eta[3]*w[2]
 +3*mu*eta[2]^2*w[2]-3*mu*eta[3]^2*w[1]-3*mu*eta[3]^2*w[2]
 +eta[2]^3*w[2]-eta[3]^3*w[1]-eta[3]^3*w[2]+mu^3+3*mu^2*eta[3]
 +3*mu*eta[3]^2+eta[3]^3-s[4], 4*mu^3*eta[2]*w[2]-4*mu^3*eta[3]*w[1]
 -4*mu^3*eta[3]*w[2]+6*mu^2*eta[2]^2*w[2]-6*mu^2*eta[3]^2*w[1]
 -6*mu^2*eta[3]^2*w[2]+4*mu*eta[2]^3*w[2]-4*mu*eta[3]^3*w[1]
 -4*mu*eta[3]^3*w[2]+eta[2]^4*w[2]-eta[3]^4*w[1]-eta[3]^4*w[2]
 +mu^4+4*mu^3*eta[3]+6*mu^2*eta[3]^2+4*mu*eta[3]^3+eta[3]^4-s[5]];

tord := plex(phi,w[1],w[2],mu,eta[2],eta[3],s[5],s[4],s[3],s[2],s[1]);

infolevel[Groebner] := 5:

# G := CodeTools:-Usage(Groebner:-Basis(F, tord));  # takes a long time

Maple's default strategy for lexicographical bases is to compute a total degree basis and convert, because direct elimination is often infeasible.  But in this case, direct elimination is much faster because the polynomials are close enough to a triangular form, and the total degree basis is large.  It took me about 2 minutes to compute directly.

G := CodeTools:-Usage(Groebner:-Basis(F, tord, method=direct));

It would be nice if Maple could detect this situation, but how to do that is unclear.  If you have computed G via the direct method, you can give it to simplify/siderels, i.e.:

result := simplify(kappa, G, tord):

I got:

[s[1]*s[2]+1, -s[2]*s[1], s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^2*s[2]-s[1]^2*s[3], -s[1]^2*s[3]-s[1]*s[2], 

s[3]*s[1]^2, 

s[1]^3*s[2]+2*s[1]^3*s[3]-s[1]^3*s[4]+s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^3*s[2]-2*s[1]^3*s[3]+s[1]^3*s[4], 

-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[2]-s[1]^2*s[3], s[1]^3*s[3]-s[1]^3*s[4], 

-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[3]-s[1]*s[2], s[1]^3*s[3]-s[1]^3*s[4], 

-s[1]^3*s[4]+s[1]^2*s[3], s[1]^3*s[4], 

s[1]^4*s[2]+3*s[1]^4*s[3]-3*s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[2]

+2*s[1]^3*s[3]-s[1]^3*s[4]+s[1]^2*s[2]+s[1]^2*s[3]+s[1]*s[2]+1, 

-s[1]^4*s[2]-3*s[1]^4*s[3]+3*s[1]^4*s[4]-s[1]^4*s[5], 

-s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[2]-2*s[1]^3*s[3]+s[1]^3*s[4], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], 

-s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[3]+s[1]^3*s[4]-s[1]^2*s[2]-s[1]^2*s[3], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[3]-s[1]^3*s[4], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[3]+2*s[1]^4*s[4]-s[1]^4*s[5]-s[1]^3*s[3]+s[1]^3*s[4]

-s[1]^2*s[3]-s[1]*s[2], 

s[1]^4*s[3]-2*s[1]^4*s[4]+s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]+s[1]^3*s[3]-s[1]^3*s[4], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[4]+s[1]^4*s[5]-s[1]^3*s[4]+s[1]^2*s[3], 

s[1]^4*s[4]-s[1]^4*s[5], -s[1]^4*s[5]+s[1]^3*s[4], s[5]*s[1]^4]

I tried mapping this to an equivalent polynomial system with this program:

polsys := proc(sys,var)
 local F,R,M,V,P,C,S,i,j:
 F,R,M,V := Groebner:-SubstituteRootOfs(sys);
 P := subs(F, sys);
 C := map(a->1/a, indets(P,anything^negative));
 C := {seq(seq(i[1],i=factors(j)[2]), j=C)};
 P := [op(map(numer,P)),op(M)];
 V := [op(var),op(V)];
 P, C, V, R;
end proc:

This function returns four outputs: a set of polynomials, a set of constraints (polynomials not equal to zero), a set of variables, and a set of substitutions to map the problem back to the original inputs.

@Carl Love 200 x 200 is a tiny grid.  If you have a 3 GHz processor running for 2380 seconds, it's consuming over 1.78 * 10^8 cycles per point.  Hopefully it's simulating the weather.  I could understand slow code on a 2000 x 2000 grid (1.7 million cycles per point), or pretty good code on a 20000 x 20000 grid (17850 cycles per point) taking that long.  But 200 x 200 is really tiny.

@J4James Can you post your code?  Maybe we can optimize it.

Software must be written to take advantage of multiple cores.  Some important routines in Maple do this, e.g. numerical linear algebra:

with(LinearAlgebra):

A := RandomMatrix(5000,5000,datatype=float[8]):

B := CodeTools:-Usage( A.A ):

And polynomial arithmetic:

f := expand((1+x+y+z+t)^25):

p := CodeTools:-Usage( expand((f+1)*(f+2)) ):

CodeTools:-Usage( divide(p,f+1) ):

In the case of higher level algorithms such as polynomial factorization below, there is parallel speedup from lower level routines but the overall speedup is limited by sequential operations and Amdahl's law.

CodeTools:-Usage( factor(p) ):

The standard interface should be better for programming tasks, not worse.  A GUI could give you a lot more information about what you're doing. That is the whole purpose of an interactive environment.  I filed a feature request. 

The standard interface should be better for programming tasks, not worse.  A GUI could give you a lot more information about what you're doing. That is the whole purpose of an interactive environment.  I filed a feature request. 

Depending on the equations, there may be different methods that can be used.  Can you post your system, and if applicable, the variables that you want to solve for (leaving other variables as unknown).

I'd be very curious to see a comparison of Maple 17 with some of the older versions, e.g. Maple V r5, Maple 8, and Maple 12, on modern hardware.  I suppose it would have to be 32-bit Linux, but you could show 64-bit Linux as well for the newer versions.

@bryon Please do something about the speed of the site.  I came back to this site after a year or so.  The response time is very annoying.  It often takes several seconds to load a page, and some queries such as "recent posts" can take 10 seconds or more.  This makes the site tedious to use and too cumbersome to check.  I also dislike the fixed width.  My browser window right now is 2000 pixels wide.  On other machines it is 600 pixels wide.  Fixed width doesn't work, and it especially doesn't work for questions.

It looks like a bug in the heuristic splitting routine.  This means there is also a problem in PrimeDecomposition.  Here is a workaround:

restart:
kernelopts(opaquemodules=false):
unprotect(PolynomialIdeals:-Radical):

PolynomialIdeals:-Radical := proc(J::PolynomialIdeal, $)
local vars, char, R, i;
option cache;
    vars := PolynomialIdeals:-IdealInfo:-Variables(J);
    char := PolynomialIdeals:-IdealInfo:-Characteristic(J);
    if char <> 0 and 0 < nops(PolynomialIdeals:-IdealInfo:-Parameters(J)) then error "non-perfect coefficient fields are not supported" end if;
    if PolynomialIdeals:-IsZeroDimensional(J) then R := PolynomialIdeals:-ideal_algorithms:-zradical(J)
    elif char = 0 then
        R := {J};
#        R := {PolynomialIdeals:-ideal_algorithms:-heusplit(J)};
        R := map(PolynomialIdeals:-ideal_algorithms:-gsplit, R);
        R := map(PolynomialIdeals:-ZeroDimensionalDecomposition, R, ':-split' = true);
        R := map(PolynomialIdeals:-ideal_algorithms:-zradical, R);
        R := PolynomialIdeals:-ideal_algorithms:-reduce(R);
        R := select(PolynomialIdeals:-IsProper, R);
        R := map(PolynomialIdeals:-Contract, R, vars);
        R := PolynomialIdeals:-Intersect(op(R))
    else error "positive-dimensional ideals in positive characteristic are not supported"
    end if;
    R
end proc:

with(PolynomialIdeals);
J := :
J2 := :
R := Radical(J): R2 := Radical(J2):
IdealContainment(J2, J, J2);
IdealContainment(R, J);
IdealContainment(R2, J2);

p := b*A*(b+B):
IdealMembership(p, R);
IdealMembership(p, R2);

You can also revert to the algorithm of Maple 15.

restart:
kernelopts(opaquemodules=false):
unprotect(PolynomialIdeals:-Radical):

PolynomialIdeals:-Radical := proc(J::PolynomialIdeal, $)
local vars, char, R, i;
option remember, system;
    vars := PolynomialIdeals:-IdealInfo:-Variables(J);
    char := PolynomialIdeals:-IdealInfo:-Characteristic(J);
    if char <> 0 and 0 < nops(PolynomialIdeals:-IdealInfo:-Parameters(J)) then error "non-perfect coefficient fields are not supported" end if;
    if PolynomialIdeals:-IsZeroDimensional(J) then PolynomialIdeals:-ideal_algorithms:-zradical(J)
    elif char = 0 then
        R := {PolynomialIdeals:-ZeroDimensionalDecomposition(J)};
        R := map(PolynomialIdeals:-ideal_algorithms:-zradical, R);
        PolynomialIdeals:-Intersect(seq(PolynomialIdeals:-Contract(i, vars), i = R))
    else error "positive-dimensional ideals in positive characteristic are not supported"
    end if
end proc:

with(PolynomialIdeals);
J := :
J2 := :
R := Radical(J): R2 := Radical(J2):
IdealContainment(J2, J, J2);
IdealContainment(R, J);
IdealContainment(R2, J2);

p := b*A*(b+B):
IdealMembership(p, R);
IdealMembership(p, R2);

@Alejandro Jakubi "Frequently, along a computation, a lot of shift may occur between the builtin expand routine and the`expand/...` routines in the library which, I guess, do not parallelize. So, isn't there a lot of overhead involved in the process? (pressumably creating and destroyng threads)."

The `expand/bigprod` and `expand/bigdiv` call into a C library (sdmp), which is parallelized, so large problems run in parallel.  The result is converted back to Maple's representation, which is a big source of overhead.  We hope to get rid of this overhead to improve parallel speedup.

@Alejandro Jakubi "Frequently, along a computation, a lot of shift may occur between the builtin expand routine and the`expand/...` routines in the library which, I guess, do not parallelize. So, isn't there a lot of overhead involved in the process? (pressumably creating and destroyng threads)."

The `expand/bigprod` and `expand/bigdiv` call into a C library (sdmp), which is parallelized, so large problems run in parallel.  The result is converted back to Maple's representation, which is a big source of overhead.  We hope to get rid of this overhead to improve parallel speedup.

Is it possible to send me these examples?  I would like to test them for our research project:

http://www.cecm.sfu.ca/~mmonagan/papers/issac12.pdf

My email address is rpearcea at cecm dot sfu dot ca

Any large Maple examples are greatly appreciated.  Worksheets or text files are ok.

Is it possible to send me these examples?  I would like to test them for our research project:

http://www.cecm.sfu.ca/~mmonagan/papers/issac12.pdf

My email address is rpearcea at cecm dot sfu dot ca

Any large Maple examples are greatly appreciated.  Worksheets or text files are ok.