An integral from quantum mechanics in Sage

Symbolic integration is tricky

Symbolic integration is tricky, both for humans and for computers! Various software systems offer symbolic integration functionality, with varying abilities to deal with all sorts of integrals.

Sage delegates the computation of many symbolic integrals to integration engines provided by other software.

Currently (in Sage 9.2) for symbolic integrals one can call one of FriCAS, Giac, Maxima, or SymPy.

The is a default choice is to call Maxima, but users can also decide which to call, via the optional parameter algorithm=....

Among the available choices, Giac, Maxima and SymPy are "standard packages", meaning they are always installed with Sage. By contrast, FriCAS is an "optional package" in Sage, meaning it does not come installed by default with Sage, but users can decide to install it additionally.

Making the problem look simpler

The integral in the question gets additional complexity from involving a number of constants each with their own quirks:

• $\hbar$ is a weird letter with its extra bar,
• $\pi$ is irrational and even transcendent,
• $2\pi$ should actually be called some name for itself,
• $\sigma_p$ has a greek letter and a superscript,
• there are square roots and $2^{1/4}$ is even a fourth root

There are extra variables t and x, an offset x_0, products, sums and quotients all over the place...

Hard to think clearly with all those distractions.

Why not make life simpler for ourselves and for the software we use by decomposing things a little bit.

The function to integrate is, really, of the form

F = A * exp(-(a*(p-p0)^2 + b*i*(p-p1)^2))

with A, a, b positive real parameters, p0, p1 real parameters of arbitrary sign, and finally our integration variable p, also real.

So we might want to try and compute

J = integral(F, p, -oo, +oo)

and if that succeeds, we can then substitute the values of the parameters to get our answer.

(To avoid any risk of confusion with the imaginary unit $i$, let us call our integral $J$ rather than "I for integral".)

Does this integral even converge?

Before jumping to computing or to asking software to compute for us, we might pause for a second and wonder if the integral actually converges.

We know the integral of exp(-p^2) from -oo to +oo converges, and neither shifting by p0 nor scaling by a can change that.

So G = exp(-a*(p-p0)^2), which by the way takes only real positive values, has a convergent integral from -oo to +oo.

Now exp(-(a*(p-p0)^2 + b*i*(p-p1)^2)) = G * H where H = exp(-b*i*(p-p1)^2), and since the integral of the real positive valued function G converges, winding it by H (which takes values on the unit circle) is not going to change that: the integral of G * H converges and its absolute value is less than the integral of G.

Good, so F converges. Let's compute it!

How various integrators deal with our integral

We want to try

J = integral(F, p, -oo, +oo, algorithm='thing')

For all the available choices of 'thing', i.e., taking them in reverse alphabetic order (why not), 'sympy', 'maxima', 'giac', 'fricas'.

Let's try that!

SymPy

Let's tell Sage to ask SymPy to compute:

sage: J = integral(F, p, -oo, +oo, algorithm='sympy')

Disappointingly, this causes an attribute error:

Traceback (most recent call last)
...
AttributeError: 'And' object has no attribute '_sage_'

Looking more carefully we could figure out what is going on: Sage is calling SymPy, which returns an answer, but converting this answer back to Sage fails.

If we know a little bit of SymPy, we can ask SymPy directly. Maybe someone will cover that in a different answer.

Maxima

Let's tell Sage to ask Maxima to compute:

sage: J = integral(F, p, -oo, +oo, algorithm='maxima')

Disappointingly, we get an error similar to the one in the question:

Traceback (most recent call last)
...
AttributeError: 'And' object has no attribute '_sage_'

Giac

Let's tell Sage to ask Giac to compute:

sage: J = integral(F, p, -oo, +oo, algorithm='giac')

This succeeds and we get the following answer:

sage: J
+Infinity

However we know the integral should have been finite. Either Giac is confused, or there was some communication problem between Sage and Giac, either when asking the question or when interpreting the answer. This would deserve to be investigated more thoroughly.

FriCAS

Let's tell Sage to ask FriCAS to compute:

sage: J = integral(F, p, -oo, +oo, algorithm='fricas')

Fricas not installed?

If FriCAS is not installed we get a long (very long) error message ending with

TypeError: unable to start fricas because the command
'fricas -nosman' failed: The command was not found
or was not executable: fricas.

In order to use the FriCAS interface you need to have FriCAS installed.
You can either run 'sage -i fricas' to install FriCAS as an optional
package within SageMath, or install FriCAS separately, see
http://fricas.sourceforge.net.

Hopefully in some future version of Sage, the message can be made a lot shorter: when running integrate(..., algorithm='fricas'), Sage could maybe check whether FriCAS is available and, if not, raise an error with a short informative error message (something like the fragment quoted above for instance).

Fricas installed?

If FriCAS is installed, the computation succeeds and we get:

sage: J
sqrt(pi)*A*e^((-I*a*b*p0^2 + 2*I*a*b*p0*p1 - I*a*b*p1^2)/(I*a + b))/sqrt(I*a + b)

which seems plausible. The sqrt(pi) reminds us of the integral of $\int_{-\infty}^{+\infty}e^{-t^2}dt$, and there are some multiplicative constants involving a, b, p0 and p1, as one would expect.

Factoring some parts of this expression, this can be rewritten as:

K = sqrt(pi)/sqrt(I*a + b)*A*e^((-I*a*b/(I*a + b)*(p0 - p1)^2))

and we can substitute all the messy constants if we wish.

Note that the Sage + FriCAS combination is installed both at SageCell and at CoCalc, so one can try it online without having to install FriCAS oneself.