Symbolic integration with an open source CAS

When I use Maple with my first year students, and we are experimenting with integration, I challenge them to produce a function which Maple can’t integrate in closed form. Given the huge number of special functions in Maple, and my students’ lack of imagination in creating functions, this is a hard exercise for them (until they get the idea of composition). Mathematica is the same; its vast library of special functions means that many integrals can be expressed in closed form. The open source CAS’s known to me: Sage, Maxima, Axiom, REDUCE, and Yacas, are generally much less able for integration. (Yacas is not being currently developed, but apparently some work is being done on a fork of it.) These systems know about the standard transcendental functions: trigonometric and hyperbolic and their inverses, exponential and logarithmic functions. However, there are a few limits even here.

I started thinking about this with the integral

displaystyle int frac{sqrt{x+sqrt{1+x^2}}}{x},dx

which can be expressed entirely in terms of elementary functions. You can see it done in Wolfram|Alpha and check it out. However, neither Maxima (and hence Sage) can solve this, and REDUCE can only do it with a special package designed for integrals of this type. Interestingly, Axiom does this with no trouble.

Some of the open source systems know a little about other functions: the error function, or elliptic integrals, but often not enough to use them as results of integration.

For example

displaystyle int e^{-x^2},dx

is solved by all of Sage/Maxima, REDUCE and Axiom in terms of the error function, but the similar function

displaystyle int e^{x^2},dx

can’t be solved in Axiom.

Sage/Maxima can solve the integral

displaystyle intsin(x^2),dx

in terms of the error function, but REDUCE and Axiom just return the integral unevaluated. Axiom’s knowledge of special functions is the poorest of all these systems, although it would probably not be hard to build such knowledge in. However, Axiom can evaluate

displaystyle int e^{e^x},dx

in terms of the exponential integral Ei(x), as can Sage/Maxima (but they express the result in terms of the incomplete gamma function), but REDUCE can’t do anything.

None of these systems can evaluate

displaystyle int frac{1}{sqrt{1-x^4}},dx = F(sin^{-1}(x),-1)

which involves an elliptic integral. However, Maxima gets at least half marks for being able to differentiate the result properly:

sage: maxima.diff(elliptic_f(arcsin(x),-1),x)

For one final example:

displaystyle intlog(sin(x)),dx

is solved by Sage/Maxima using the polylogarithm function; Axiom and REDUCE can do nothing.

Maxima and REDUCE come with test suites for integration, but I don’t think there’s such a suite for Axiom. There has not been, as far as I know, a published comparison of the integration capabilities of these systems, but it seems that of the systems mentioned here, Maxima has the best all round strength.

7 thoughts on “Symbolic integration with an open source CAS

  1. I just want to point out that there is also sympy. If you have Sage you can call it via the algorithm parameter.

    sage: ex(x)= sin(x)
    sage: ex.integrate(x,algorithm=’sympy’)
    x |–> -cos(x)

    or directly in sympy of course …

  2. Axiom is usually very good at integrating any function, f(x), for which there exists an elementary function, F(x), such that dF(x)/dx = f(x). In fact, using Axiom it is possible to determine if an elementary function, F(x) exists at all — even if Axiom is unable to determine its form.

    This is because Axiom implements the Risch algorithm, or one of its many revisions. Using the Risch algorithm it is possible to: integrate any function for which the integral can be expressed in terms of elementary functions and if not, prove that no such function, F(x) exists. It is a decision procedure.

    However, it is very, very difficult to implement. Axiom, unlike any other software I know implements the complete negative case: i.e., it knows the difference between “I can’t do it” and “no one can do it.”

    Maxima partially implements the Risch algorithm; as does SymPy (a simplified version at least). Due to its complexity most software will use a combination of look-up tables first before hitting their Risch implementations.

    Now, lets get onto integrals for which no elementary function exists — as is the case for many of your examples. The Risch algorithm is useless here with many pieces of software just falling back to table look-ups and pattern patching. This is why Axiom performs so poorly and Maxima not much better. I suspect Maple and Mathematica use the properties of the Meijer G-function in order to handle such non-elementary integrals.

    But, indefinite integrals are boring; if the function can not be represented in terms of elementary functions chances are I don’t care. Why? As Ei, Erf and friends are not nice to manipulate symbolically and often a real pain to evaluate. Most nasty integrals which I, as a Physics student, stumble upon are definite integrals. The limits are often between -inf and inf. (And, if you take a look at G&R you’ll find 50% of the integrals are definite integrals where one of the limits is +/- inf or +/- pi.) Try getting Axiom, Maxima and friends to evaluate them.

    (You’ll most probably find that Maxima has no real concept of infinity — failing to evaluate integral of exp(-x^2) between -inf and inf while SymPy, thanks to its impressive limit code, performs a lot better. I am unsure about Axiom.)

  3. Very good points – thank you very much! As for Axiom, here’s the integral:

    (1) -> integrate(exp(-x^2),x=%minusInfinity..%plusInfinity)

    (1) |%pi
    Type: Union(f1: OrderedCompletion(Expression(Integer)),…)

    Axiom in fact handles infinity – and other cardinals – very well; the only CAS, to my knowledge, to do so.

  4. Well I have tried exp(-x^2) between -inf and +inf with Maxima and it answers
    sqrt(pi). command:

  5. strange enough with Maxima, if it is entered the command
    integrate( ( x + ( 1 + x^2 )^1/2 )^1/2 * 1/x, x)
    a result should be given

  6. Ernesto: I’ll necro this thread to point of that you’re in fact raising to the 1 power and dividing by 2. The equivalent command would be integrate( ( x + ( 1 + x^2 )^(1/2) )^(1/2) * 1/x, x).

Leave a Reply

Your email address will not be published. Required fields are marked *