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.
Welcome to Ask Sage! And to Sage!
Thank you for your question!
f(p) looks like a normal law. So why the exponent 1/4 on the constant.
@Cyrille I think the constants are generally required in these integrals for normalization purposes. For e.g., in this case, see here: https://www.youtube.com/watch?v=wcrXc...
Which version of sage are you using? Sage 9.2 is using an updated version of Maxima. I get the same kind of errors in 9.0. Sometimes Maxima can't even show that an exponential is positive, probably because symbolic variables are complex by default in sage, but using assumptions gives a good old segfault.
I am indeed using Sage 9.0. Surely there can't be a huge difference between 9.0 and 9.2? I'll definitely consider updating as a last option if I have no other solutions since I have sage working for everything else I need.