Friday, July 31, 2015


Maxima versus Mathematica


Should I go for Maxima instead of Mathematica? NO! DON'T!

    Prolog  [top]

    General purpose Computer Algebra Systems (CAS) have become a part of many scientist's toolboxes. In natural sciences and in part in engineering, by far the most wide-spread CAS are  Mathematica and Maple. These are proprietary software systems which rest heavily on the budget of any university department. Yet, these products are so forcefully present on any campus, that it is practically unavoidable to get dependent on them early on, not only after diving into actual research, but also on the level of course homework and in graduate projects.

    A natural question to ask is, whether there are open-source alternative to these CAS. The answer may depend on your needs - see http:// en.wikipedia.org/ wiki/ List of computer algebra systems. There are some highly specialized systems out there which are of little help to natural sciences and/or engineering. Pure mathematician may find their poison there. They are no part of this post. I have no clue what these systems are good for.

    I.m.h.o. most CAS users in science and engineering do things related to expression evaluation and simplification in its wides sense, including symbolic high-school math, linear algebra, all forms and usage of calculus in any dimension, differential equations, differential geometry, and analysis in the complex plane. A smaller fraction of  CAS codes involves home-grown symbolic pattern manipulation, using pattern matching, substitution, and rule-based algorithms (these can be pretty cool). Based on the preceding my very personal view is, that there are essentially three open source alternatives one might consider: (i) Maxima, (ii) Sage, (iii) Axiom. I have no experience with the latter. It seems less oriented towards natural sciences. There may be also be a new kid in town, i.e. (iv) SymPy, but it is much too early to say. Sage is not really a CAS only. It is an attempt at a common Python-istic interface to as many open source mathematics packages as possible. The main point is, that if it comes to science and engineering CAS type-of-work, it is primarily Maxima under the hood of Sage which is doing all the work. In turn, for me, there is really only one open source competitor, namely Maxima.

    This blog tells you about some of my Maxima experience. Don't get me wrong, this is not some geeky, full fledged test of Maxima internals, rather it is like always in science: first I'd like to see if Maxima can do basic stuff right, before I'll trust it to solve a quantum mechanics problem. What I learned is: Maxima gets elementary math wrong, it knows very little about integration, it is weak on simplifying expressions, it fails on trivial equations, it is slow, Maxima usage is discouraged by CAS power users. In summary: do not go and waist your time or wreck your nerves using Maxima. JUST DON'T DO IT. So, let me start

    Maxima gets elementary math wrong  [top]

    One of my very first experiences with Maxima involved the branch cut behavior of an elementary function in the complex plane, namely the natural logarithm. We all know, that in principle one can choose whichever crazy angle one likes to set the natural log's branch cut, but in practice most people attach it to the negative real axis. This seems also true for Maxima
    (%i) log(-2 + 0.001*%i);
    (%o) 0.6931473055599296+3.14109265363146*%i
    (%i) log(-2 - 0.001*%i);
    (%o) 0.6931473055599296-3.14109265363146*%i
    Mathematica follows the same convention. Now let's look a simple symbolic calculation using this property$$a,b\in\mathrm{reals},\hphantom{aa} a>b>0$$ $$\lim_{\eta\rightarrow 0^+} \ln (b-a+i\eta) = \ln (a-b) + i \pi \hphantom{aaaa} \lim_{\eta\rightarrow 0^-} \ln (b-a+i\eta) = \ln (a-b) - i \pi$$Watch the different order of $a,b$ on both sides of the equations. Ok., so here is Maxima's result 
    (%i) declare(a,real);
    (%o) done
    (%i) declare(b,real);
    (%o) done
    (%i) assume(a>b,b>0);
    (%o) [a > b,b > 0]
    (%i) limit(log(b-a+%i*eta),eta,0,plus);
    (%o) log(b-a)+%i*%pi
    (%i) limit(log(b-a+%i*eta),eta,0,minus);
    (%o) log(b-a)-%i*%pi
    Both of the last outputs are simply wrong. Now, here is Mathematica's result
    In[]:= Limit[Log[b - a + I eta], eta -> 0, Direction -> -1,Assumptions -> {a > 0, b > 0, a > b}]
    Out[]= I Pi + Log[a - b]
    In[]:= Limit[Log[b - a + I eta], eta -> 0, Direction -> 1,Assumptions -> {a > 0, b > 0, a > b}]
    Out[]= -I Pi + Log[a - b]
    Obviously, both outputs are correct.
    In many areas of natural sciences, employing the limiting properties of much more complicated functions in the complex plane is of tantamount importance.  If you want to base your next research paper on 'results' from Maxima - I hope you sleep well. To me things like this are really frightening.  How can you ever trust a CAS that doesn't even get such simple things straight.

    Maxima knows very little about integration  [top]

    Certainly calculus is more than just differentiation and integration. Yet these are among the most basic ingredients of calculus. So a CAS better be good at them. The famous mathematician Hilbert once said that "differentiation is a plumbers job, while integration is an art". So let us look at integration only.
    Many online demonstrations of Maxima will show you trivial integration of mostly rational functions for which I really don't know why I would need a CAS to begin with.
    The very moment one departs from such simple stuff Maxima routinely fails. Here are some examples.

    For starters, let's try this one$$\int_0^1 \frac{1}{(x^2-x^3)^{1/3}}\, dx$$Maxima provides no solution
    (%i) integrate(1/(x^2-x^3)^(1/3),x,0,1);
    (%o) 'integrate(1/(x^2-x^3)^(1/3),x,0,1)
    Mathematica does the job
    In[]:= Integrate[1/(x^2-x^3)^(1/3),{x,0,1}]
    Out[]= 2 Pi/Sqrt[3]

    Next is a so-called n-th order Bose integral$$I_n = \int_0^\infty \frac{x^n}{e^x-1}\,dx$$Let's look at $n=1$. Maxima provides no solution
    (%i) integrate(x/(%e^x-1),x,0,inf);
    (%o) 'limit(li[2](%e^x)+x*log(1-%e^x)-x^2/2,x,inf,minus)-%pi^2/6
    Mathematica does the job
    In[]:= Integrate[x/(Exp[x]-1),{x,0,Infinity}]
    Out[]= Pi^2/6
    Needless to say that Mathematica will properly use recursion relations to produce the result for any $n$, say for $n=10$ we get $I_n=3628800\, \zeta(11)$, while Maxima will get lost in ever longer expressions involving polylogs.

    Ok., one more$$\int_{-\infty}^\infty \frac{1}{x^2+1}\frac{1}{e^x+1}\,dx$$Maxima provides no solution
    (%i) integrate(1/(x^2+1)/(exp(x)+1),x,minf,inf);
    (%o) 'integrate(1/((x^2+1)*(%e^x+1)),x,minf,inf)
    Mathematica does the job
    In[]:= Integrate[1/(x^2+1)/(Exp[x]+1),{x,-Infinity,Infinity}]
    Out[]= Pi/2
    Again, there are more complicated versions of this type of integral, Mathematica will continue to solve, while Maxima just fails.

    PS.: just for the weary, all of the preceding integral are a little more complicated than trivial rational integrals in the sense that they involve fun-games in the complex plane to actually be evaluated.

    Finally, and even though CAS designers may hate this, one of the prime reasons for natural scientists and engineers to use CAS is to have voluminous handbooks of mathematical functions and integral tables just at the tip of their fingers (=keyboards).
    Therefore many researchers tend to assess the quality of a CAS by checking at how good it is to recognize all of the integrals in the two seminal handbooks AS (for imperialists ;)) and GR (for commies ;)).
    To make a long story short, Maxima fails mega-miserably at this, to say the very least. It has some special functions and their derivatives and series representations implemented, but that does not seem to help. Just one trivial example: elliptic integrals$$K(x)=\int_0^{\pi/2}(1-x\sin(y)^2)^{-1/2}\,dy$$Maxima provides no solution
    (%i) integrate(1/sqrt(1 - m*sin(x)^2), x, 0,%pi);
    (%o) 'integrate(1/sqrt(1-m*sin(x)^2),x,0,%pi)
    Mathematica does the job
    In[]:= Integrate[1/Sqrt[1 - x Sin[y]^2], {y, 0, Pi/2}]
    Out[]= ConditionalExpression[EllipticK[x/(-1 + x)]/Sqrt[1 - x],
           Im[1/(-1 + x)] != 0 || Re[1/(-1 + x)] <= 0]
    You can spend a lot of time in going through AS and GR - with essentially the same picture for any special function you look at.

    These where just four personally biased examples, but in a nutshell, when you start really using Maxima, you will soon be very frustrated to find that is extremely weak at integration and will very often return ... just nothing.

    Maxima is weak on simplifying expressionn  [top]

    Simplification of expressions is an 'evergreen' in the CAS community. This is because  - like it or not - simplification of expressions remains one of the prime tasks for which CAS are used in natural sciences and engineering. However, what looks 'simple' to many scientists and engineers may not conform with CAS' designers puristic definitions of such.
    Therefore it comes as no surprise that claims of a failed simplification, which are directed towards the maintainers of a CAS by a frustrated user, tend to be refuted immediately.
    Dismissing such dogmatic discussions, I like to follow a practitioner's approach to see how well a particular CAS is at what I believe are trivial simplifications. For this I give it a totally meaningless, sufficiently large, expression, which only contains elementary operations and elementary functions, for which I know, that if properly arranged, using only basic math, it will look MUCH simpler beyond any doubt.
    So, let's chose this one:
    (%i) B:1/4*(%e^(sqrt(-2*cos(l) + 2)) + 1)*(%i*sin(l) - cos(l))*%e^(-1/2*sqrt(-2*cos(l) + 2)) - 1/4*sqrt(-2*cos(l) + 2)*(%e^(sqrt(-2*cos(l) + 2))*sin(l)^2/sqrt(-2*cos(l) + 2) + %i*%e^(sqrt(-2*cos(l) + 2))*sin(l)*cos(l)/sqrt(-2*cos(l) + 2) - (%i*%e^(sqrt(-2*cos(l) + 2)) - %i)*sin(l) + (%e^(sqrt(-2*cos(l) + 2)) - 1)*cos(l) - %i*%e^(sqrt(-2*cos(l) + 2))*sin(l)/sqrt(-2*cos(l) + 2))/(%e^(1/2*sqrt(-2*cos(l) + 2))*cos(l) - %e^(1/2*sqrt(-2*cos(l) + 2))) + 1/8*sqrt(-2*cos(l) + 2)*((%i*%e^(sqrt(-2*cos(l) + 2)) - %i)*cos(l) + (%e^(sqrt(-2*cos(l) + 2)) - 1)*sin(l) - %i*%e^(sqrt(-2*cos(l) + 2)) + %i)*(%e^(1/2*sqrt(-2*cos(l) + 2))*sin(l)*cos(l)/sqrt(-2*cos(l) + 2) - 2*%e^(1/2*sqrt(-2*cos(l) + 2))*sin(l) - %e^(1/2*sqrt(-2*cos(l) + 2))*sin(l)/sqrt(-2*cos(l) + 2))/(%e^(1/2*sqrt(-2*cos(l) + 2))*cos(l) - %e^(1/2*sqrt(-2*cos(l) + 2)))^2 - 1/4*((%i*%e^(sqrt(-2*cos(l) + 2)) - %i)*cos(l) + (%e^(sqrt(-2*cos(l) + 2)) - 1)*sin(l) - %i*%e^(sqrt(-2*cos(l) + 2)) + %i)*sin(l)/(sqrt(-2*cos(l) + 2)*(%e^(1/2*sqrt(-2*cos(l) + 2))*cos(l) - %e^(1/2*sqrt(-2*cos(l) + 2))));

    Feeding this into Mathematica, including the proper syntax adjustments, i.e. cos()->Cos[] a.s.o., it does a terrific job
    In[]:= FullSimplify[B]
    Out[]:= (I*Sin[l/2]^3*(Sqrt[1 - Cos[l]]*Cosh[Sin[l/2]] - 
            Sqrt[2]*Sinh[Sqrt[Sin[l/2]^2]]))/(E^((I/2)*l)*(1 - Cos[l])^(3/2))

    Regarding Maxima, I'm not going to waste blog-space for the many amusing 'results' you can get, if you apply any combination of its built-in 'simplification' algorithms. You will realize that Maxima does not have a unified simplification interface like Mathematica, but a plethora of different routines, involving modules like ratsimp, fullratsimp, radcan, logcontract, rootscontract, radexpand, trigsimp, trigexpand, trigreduce, and trigrat and many more of that. For the user it remains a frustrating trial-and-error game to find out if any of them might be of help. The bottom line is: Maxima fails to produce anything simpler from expression B - even if assessed by very mild standard. It fails similarly for many other simplification benchmarks of this type.

    Maxima cannot solve simple equations  [top]

    Many people use CAS to free them from the tedious task of finding all solutions to sets of equations. For linear systems of equations there are standard algorithms, any CAS will have implemented more or less safely. However, if it comes to non-linear equations the wheat gets separated from the chaff. So, by example, lets look at two very simple non-linear sets of equations.

    Number 1: $\hspace1cm \sin(x-y)-\sin(x)=0,\hspace1cm\sin(y)=0$

    You know the maifold of solutions ... right ;-) ? Now, let's try Maxima. There are various types of 'solvers'. Let's first try its  basic solver, which is available immediately on start up
    (%i) %solve([sin(y) = 0,sin(x-y)-sin(y) = 0],[x,y])
    (%o)[]
    so nothing. Can you believe this ... its not rocket science ... its just basic high-school math ... nothing!!!. Ok. then, let's try their next built-in solver
    (%i) load(to_poly_solve)$
    (%i) to_poly_solve([sin(x-y)-sin(y)=0,sin(y)=0],[x,y]);
    Unable to solve
    Unable to solve
    Unable to solve
    (%o) %solve([sin(y) = 0, (- sin(y - x)) - sin(y) = 0], [x, y])
    wow, great, still nothing ... but there is yet another 'super'-solver which is available
    (%i7) load(solver);
    (%i11) Solver([sin(x-y)-sin(y)=0,sin(y)=0],[x,y]);
    solve: using arc-trig functions to get a solution.
    Some solutions will be lost.
    solve: using arc-trig functions to get a solution.
    Some solutions will be lost.
    (%o11) [[x = 0, y = 0]] 
    so, depressingly, even this last resort returns only a partial solution. Btw., this solver routine is an extension of Maxima written for a Diploma thesis in 1994, in Germany, and completely unmaintained ever since then ... this gives me the shivers. Now, lets try Mathematica
    In[]:= Solve[Sin[y]==0 && Sin[x-y]-Sin[y]==0,{x,y}]
    Out[]= {{x -> ConditionalExpression[2 Pi C[1], C[1] ∈ Integers],
         y -> ConditionalExpression[-2 Pi C[2], C[2] ∈ Integers]}, 
        {x -> ConditionalExpression[2 Pi C[1], C[1] ∈ Integers], 
         y -> ConditionalExpression[-Pi - 2 Pi C[2], C[2] ∈ Integers]}, 
        {x -> ConditionalExpression[Pi + 2 Pi C[1], C[1] ∈ Integers], 
         y -> ConditionalExpression[-2 Pi C[2], C[2] ∈ Integers]}, 
        {x -> ConditionalExpression[Pi + 2 Pi C[1], C[1] ∈ Integers], 
         y -> ConditionalExpression[-Pi - 2 Pi C[2], C[2] ∈ Integers]}}
    This what I'd call a complete solution. You may want to see solutions, removing some trivial multiples due to trigonometric function's periodicity
    In[]:= Solve[Sin[y]==0 && Sin[x-y]-Sin[y]==0,{x, y}]/.{C[1]->0,C[2]->0}
    Out[]= {{x->0,y->0},{x->0,y->-Pi},{x->Pi,y->0},{x->Pi,y->-Pi}}
    which clearly shows you that Maxima truly misses solutions.

    Depending on the problem, there may be ways of rewriting sets of equations, to help Maxima to eventually find all solutions. E.g. in the present case we could use that $\sin(x)=(z-1/z)/(2i)$ with $z=e^{ix}$ and turn the set of equations into a set of polynomial equations, for which Maxima will hopefully return all solutions. But that is not at all the issue: a CAS should give the solution to you - not the other way around.

    Number 2:$\hspace1cm\sin(x/2)\cos(\sqrt{3} y/2)+\sin(x) = 0,\hspace1cm \sqrt{3}\cos(x/2)\sin(\sqrt{3} y/2) = 0$

    Maxima finds this
    (%i)a:[sin(x/2)*cos(sqrt(3)*y/2)+sin(x)=0,sqrt(3)*cos(x/2)*sin(sqrt(3)*y/2)=0]$
    (%i)Solver(a,[x,y]);
    solve: using arc-trig functions to get a solution.
    Some solutions will be lost.
    solve: using arc-trig functions to get a solution.
    Some solutions will be lost.
    (%o)[[x = %pi,y = %pi/sqrt(3)]]
    yep(!) true indeed :( , at least it warns you that again it finds, only a partial solution. Now, lets try Mathematica, already suppressing some periodical solutions
    In[]:=Solve[Sin[x/2]Cos[Sqrt[3]y/2]+Sin[x]==0&&Sqrt[3]Cos[x/2]Sin[Sqrt[3]y/2]==0,{x,y}]/.{C[1]->0,C[2]->0}
    Out[]:={{x->0,y->0},{x->0,y->(2*Pi)/Sqrt[3]},
    {x->(-4*Pi)/3,y->0},{x->-Pi,y->-(Pi/Sqrt[3])},
    {x->-Pi,y->Pi/Sqrt[3]},{x->(-2*Pi)/3,
    y->(2*Pi)/Sqrt[3]},{x->(2*Pi)/3,y->(2*Pi)/Sqrt[3]},
    {x->Pi,y->-(Pi/Sqrt[3])},{x->Pi,y->Pi/Sqrt[3]},
    {x->(4*Pi)/3,y->0},{x->2*Pi,y->0},
    {x->2*Pi,y->(2*Pi)/Sqrt[3]}}
    there are still more solutions returned then necessary if one invokes symmetry - but better more than less. So, again Mathematica returns all solutions, several of which Maxima misses, e.g. the one at $(x,y)=(0,0)$ a.s.o..

    In conclusion: at its very best, Maxima will find 'some' of the solutions of even very simple non-linear sets of equations. What if your science project is in need of exactly one of those solutions it did not find?

    Maxima is slow  [top]

    There are rather diverse dogmas on the performance of the ancient language underlying Maxima, i.e. LISP. I'd like to stay away from these and rather select only one big issue. For many practitioners in natural sciences and engineering one of the very reasons for using computers is to set up many nested loops of action, accessing and working in some way on big sets of data arranged in various types of containers. Now, Mathematica uses packed arrays for container management, while Maxima is based on linked lists (... LISP ;-) ). Therefore access/operations to/on anything ranging from vectors, to matrices, to multidimensional arrays can be painfully slow in Maxima, meaning that Maxima tends to be slower by a factor of $\sim N$, or $\sim\log(N)$ at best, as compared to Mathematica. Here $N$ roughly refers to the number of 'objects' stored in your container. (Having said that, we all know that there are some rare algorithmic cases where linked lists are the way go, but that's not my point.)

    Maxima usage is discouraged by CAS power users  [top]

    There exist some huge  CAS-based projects which try to solve tomorrows natural science problems with todays hard- and software. Their developers surely have a much more profound insight into the strength and weaknesses of various CAS systems, than I will ever have.
    In this context one remarkably open-hearted statement is due to the core developer of  FeynCalc, a powerful package for algebraic calculations in elementary particle physics. In a nutshell: FeynCalc had actually begun as a MacSyma code - a closed source version of Maxima. But soon after that the very core developers ditched the MacSyma code in favor of  a complete Mathematica rewrite - primarily because of the exceedingly better performance of the Mathematica version - and this even with some of the very early versions of Mathematica.
    Later attempts of the Maxima community to simulate a roll back of the code have been given a rather disillusioning and well-funded 'don't even think about it' by the FeynCalc developer.
    I view such actual power users experiences far more informative then the many introductory "comparisons" on the web, like e.g. The Playsome Threesome: Maxima, Maple, and Mathematica, which try to may make you believe that on the average there is little difference between all of these CAS. Which is just not true.

    Epilog  [top]

    All of the issues reported here (and more) have been posted by me to official user forums and mailing lists of Maxima and also SageMath. Some of them more then three years ago. Since then nothing has changed about them. Arguably, it might take more time to strengthen things like Maxima's integration capabilities dramatically. But other problems, like accepting Maxima to fool its users with a broken behavior of elementary functions in the complex plane is remarkable to me, to say the least.

    Companies like Wolfram Alpha employ hundreds of professional natural scientists, mathematicians, and programmers full time for nice salaries, to write, expand, and maintain a CAS like Mathematica. The 20-40 open source developers who work on Maxima in their spare time, or even as retirees, are certainly heroes. Some of them are even very friendly and helpful heroes, as you can experience on maxima-discuss@lists.sourceforge.net. Yet, given this kind of a dramatic difference in man and financial powers, Maxima will continue to fall behind Mathematica and, most likely, as of today, has lost any chance of catching up again.

    PS.: For a kind-of similar point of view, see this blog-post by the founding father of Sage on the issue of "Sage is Failing".
    Posted by at
    Labels: , , , , , , , , , , ,

    No comments:

    Post a Comment

    Subscribe to: Post Comments (Atom)