John May

Dr. John May

2351 Reputation

17 Badges

12 years, 338 days
Pasadena, California, United States

Social Networks and Content at

Maple Application Center

I am a Senior Developer in the Mathematical Software Group and have been with Maplesoft since 2007. I am also an Adjunct Assistant Professor in the School of Computer Science at the University of Waterloo.

I have a Ph.D in Mathematics from North Carolina State University as well as Masters and Bachelors degrees from the University of Oregon. I have been working on research in computational mathematics since 1997.

My main research interests in are computational linear and polynomial algebra, especially numerical polynomial algebra. I currently work on the exact algebraic solvers as well as other subsystems of Maple.

MaplePrimes Activity

These are answers submitted by John May

 Take a look at the ImageTools package, in particular, the Read function.  For an example of using ImageTools, you can also take a look at this blog post.


 Looks like these messages are being printed on the 'hints' infolevel. 

infolevel[hints] := 0;

should turn them off.


Since you are looking for a floating point answer, and all your expressions have floating point numbers, you should really use  fsolve instead of solve.  That is, replace T= solve(f=0,T); with T = fsolve(f,T);

If you use solve, your f is converted to (notice no floating point values):


                                  /                   173063     \
                                  |40356350000 - ----------------|
                                  |                  /    116713\|
                                  |              100 |T + ------||
                                  \                  \     500  //
                      657894737 10
            f := -1 + --------------------------------------------

which Maple tries to solve exactly (it then calls evalf on the answer in the end so you don't notice), and that appears to cause big expression blow up.


You can of course do this with ListTools as well:



FindRepetitions removes one copy of each unique element which is why you get the ordering of the elements you see here.

Lots of useful stuff in ListTools.


Strangely, even the output of Groebner[Solve] does not get solve() to return a solution even though Groebner clearly says the system is two dimensional.

I was able to use triangular decomposition to get an answer however:

sol:=Triangularize(sys, R);
solve(Equations(sol[1], R));

This gave me the following solution with x8 and y2 free:

      15751       7813
{x3 = -----, x2 = ----, x8 = x8, y2 = y2,
       125        125

         -7812               15876
    x1 = -----, y1 = 0, x4 = -----, x9 = RootOf(
          125                 125

               2             2
    15625 y2 _Z  + (-15625 y2  - 61027344) _Z

     - 15625 y2), x7 = -----, x5 =

             2                               2
    125 y2 x8  - 15876 y2 + 15876 x8 - 125 y2  x8
                       125 x8

         125 x8  - 15876
    x6 = ---------------}
             125 x8

John May
Developer, Maplesoft

Setting the environment variable _MaxSols can force solve to return only one solution:

This can also make solve run faster internally in some, but certainly not all, cases.


One posibility is to avoid simplifying entirely.  In this case, you can compute the symbolic matrix "Du" by turning off normalization with:

Normalizer:=x->x; # Tells LinearAlgebra not to call normal() on intermediate results
Testzero:=testeq; # Tells LinearAlgebra to use the testeq(foo) instead of evalb(Normalizer(foo)=0)

Du := -(G21.(G11-Pt)^(-1).G12-G22)^(-1): # <- do not try to print this result

But now, Du is extremely big:

> map(length, Du)[1,1];


It would be a very bad idea to try to print or simplify Du, but it is diagonal:

> map(testeq, Du);

	[ false  true   true  ]
	[ true   false  true  ]
        [ true   true   false ]

If you are careful, you might be able to do something with this result, but since it is so large, it is going to be tough.


Assuming eq1 is expanded out into sums of the form c*y^n you could do

pp, c := selectremove(has, [op(eq1)], y); # c is a list of constant terms
pw, pd := selectremove(type, pp, `^`); # pd is products, pw is powers
pd, lin := selectremove(type, map2(op, -1, pd), `^`); # lin is the linear terms in y
[op(map(0,c)), op(map(1, lin)), op(map2(op, 2, pw)), op(map2(op, 2, pd))];

Of course, unlike Robert's solution, this won't preserve the order of the exponents.

I haven't really dug into the source of eliminate, but it looks like it is trying to preserve the variables in the lhs of the original equations as much as possible in the output equations.

It looks like you can trick it into given the form you want by having it eliminate the squares of the original equations:

elm:=eliminate({x^2=(cos(v) *cos(u))^2, y^2=(cos(v)*sin(u))^2, z^2=sin(v)^2},{u,v});

The downside to this, is that it gives 4 answers with different substitions (but with the same output "eliminated" equation) which makes the output slightly annoying to process. In this case, this works:

op(map2(op,2, {elm}));


John May
Mathematical Software

Actually, the solutions returned in Maple 9.5 are not complex, they just happen to be written in terms of complex numbers.  If you want to get the solutions written without complex numbers you can do something like:

evalf[50](sol); # Imaginary parts show up, but they are numerically 0.
evalf(sol); # no imaginary parts

In Maple 11, the fact that RealDomains package has a hard time figuring out that all the roots are real and it only returns the ones that it can tell for sure are actually real. This is a weakness in the RealDomains package (it doesn't use evalc).

John May
Mathematical Software

If you change your '.'s to '*', it should work fine (as long as you change eq1 to eqn1 in the call to solve).

eqn1 := x + y + z = 0;
eqn2:= a*x + b*y + c*z=0;
eqn3:= b*c*x+c*a*y+a*b*z=1;

solve( {eqn1,eqn2,eqn3},{x,y,z} );

The easy answer to the question:

So the question is, how do I get ReducedRowEchelonForm to stop doing this?

is probably this:


I don't think you can avoid at least some discussion of the difference between exact numbers and floating point numbers, however.  That is, to Maple (and most computer software), 27.6 is very different from 276/10.

John May
Mathematical Software

You can convert the piecewise function into a list.  Look at ?convert,pwlist

In this list every other entry will be a function, and the others will be the bounds.  Depending on whether the domain of the piecewise function is bounded, the first entry of the list may be a function or a bound.  Assuming the first entry is a polynomial, you can do the following:

#f is your original piecewise function
flist:=[seq](ftmp[i], i=1..nops(ftmp), 2); #may need to start at i=3 instead
# flist is a list of all the polynomials.

John May
Mathematical Software

solve() isn't really the right tool for this, but it should be able to find the same solution that fsolve does: P = 0.5772706212 (instead of finding nothing).

Solve works on equations with floats by running convert/rational on them, then running evalf on the output.  Part of what trips this process up is the numerical error in the argument to tan() in the input expression.    If the input expression is written only in terms of cosine:

.2666845389e-11*Pi^2/(.1450695744e-10*Pi+.1281254481e-13*Pi*cos(.5637758661e-1 *P^(1/2)))=P

then solve works.

John May
Mathematical Software

If you are only interested in computing over the real numbers, you can use: with(RealDomain): at the beginning of your Maple session.  This package redefines a number of Maple commands, including solve, to ignore complex answers.

John May
Mathematical Software

First 7 8 9 10 Page 9 of 10