[**Update**: The bug **seems** fixed in the latest version, 10.0.2.]

I am in my third year teaching a course in Quantum Mechanics, and we spend a lot of time working with a very simple system known as the harmonic oscillator — the physics of a pendulum, or a spring. In fact, the simple harmonic oscillator (SHO) is ubiquitous in almost all of physics, because we can often represent the behaviour of some system as approximately the motion of an SHO, with some corrections that we can calculate using a technique called perturbation theory.

It turns out that in order to describe the state of a quantum SHO, we need to work with the Gaussian function, essentially the combination `exp(-y²/2)`

, multiplied by another set of functions called Hermite polynomials. These latter functions are just, as the name says, polynomials, which means that they are just sums of terms like `ayⁿ`

where `a`

is some constant and `n`

is 0, 1, 2, 3, … Now, one of the properties of the Gaussian function is that it dives to zero really fast as `y`

gets far from zero, so fast that multiplying by any polynomial still goes to zero quickly. This, in turn, means that we can integrate polynomials, or the product of polynomials (which are just other, more complicated polynomials) multiplied by our Gaussian, and get nice (not infinite) answers.

Unfortunately, Wolfram Inc.‘s Mathematica (the most recent version 10.0.1) disagrees:

The details depend on exactly which Hermite polynomials I pick — 7 and 16 fail, as shown, but some combinations give the correct answer, which is in fact zero unless the two numbers differ by just one. In fact, if you force Mathematica to split the calculation into separate integrals for each term, and add them up at the end, you get the right answer.

I’ve tried to report this to Wolfram, but haven’t heard back yet. Has anyone else experienced this?

## 3 responses to “Oscillators, Integrals, and Bugs”

The example works with Mathematica 9 for me.

I’ve caught at least a couple of Mathematica bugs over the years — google seems to confirm that they’re not very unusual. I haven’t seen *this* one, but I haven’t dealt with Hermite polynomials in a long time. I have noticed based on my own experience and other’s reports that Wolfram is not particularly good about rapid fixes, which is bothersome given that so many people depend on Mathematica for their research and could be seriously led astray by defects.

Anyway, there’s a Mathematica “Stack Exchange” site, and I suggest asking there about whether others have observed the behavior. They even maintain a category of confirmed bugs…

In case you are interested, I would suggest you to take a look at sagemath.org, a free an open source (as science should be) mathematical environment with a great many features. The following input

var(‘y’)

print integrate(hermite(7,y)*y*hermite(16,y)*exp(-y^2),y,-infinity,infinity)

print integrate(hermite(7,y)*y*hermite(8,y)*exp(-y^2),y,-infinity,infinity)

gives

05160960*sqrt(pi)

Cheers!