Animate the scattering of a Wave Packet

Question

I know mathematica is probably not the best choice for intense numerical integration, but its the only software I know.

I would like to create an animation (not real-time, but pre-render the frames) of a Gaussian wave-packet scattering against a spherical square-well, whose stationary states are known. Basically, I need to carry out the following integral:

$$\Psi(\vec{x},t)=\int\frac{d^3\mathbf{k}}{\sqrt{(2\pi)^3}}\tilde{f}(\mathbf{k}) \psi_\mathbf{k}(\vec{x})\,e^{-i E(\mathbf{k})t}$$

where $\tilde{f}(\mathbf{k})=\exp(-\frac{(\mathbf{k}-\mathbf{k}_0)^2}{2(\Delta k)^2}-i \mathbf{k}.\vec{x}_0)$ is the Gaussian wavepacket in Fourier space, and $E(\mathbf{k})$ is a very simple function quadratic in $\mathbf{k}$.

The integral needs to be done at several different times, and then compiled together to form an animation (this last step I know how to do). Is there a good reference that I can look up that will accomplish my task?

Answer 1

Not knowing much about quantum mechanics, there's a lot going on in your integral that I know nothing about. The basic ideas are the same for any sequence of objects generated by numerically computed integrals, though. Suppose, for example, you want to estimate $$\iiint_{{\mathbb R}^3} e^{-\langle x,y,z \rangle \cdot \langle x,y,z \rangle / a} dV$$ for a range of $a$ values and the animate the results. We could proceed as follows.

Monitor[
  results = Table[Check[
  {result = NIntegrate[Exp[-{x, y, z}.{x, y, z}/a],
     {x, -Infinity, Infinity},
     {y, -Infinity, Infinity},
     {z, -Infinity, Infinity}]},
  {result, "Uh-oh", a}], {a, 0.1, 10, 0.1}],
a]; // AbsoluteTiming

(* Out: {56.448363, Null} *)

There are a few things going on here. The results might take a while, so I'm using Monitor to watch the progress. Also, I'm using Check to check for errors and returning the parameter that yields the error. Examining the results, it appears that the error couldn't be too bad.

ListPlot[First /@ results]

Using Table, as done here, pre-generates a list of results that can be used to generate a smooth animation with ListAnimate like so:

ListAnimate[First /@ results]

Of course, this just generates the results; you'd want a function that turns those results into pictures, whatever that is.

Also, more parameters get's trickier. Four parameters varying smoothly could take a while but you can check. Consider the following variation.

Monitor[
 results = Table[Check[
  {result = NIntegrate[Exp[-{x, y, z}.{x - x0, y - y0, z - z0}/a],
     {x, -Infinity, Infinity},
     {y, -Infinity, Infinity},
     {z, -Infinity, Infinity}]},
  {result, "Uh-oh", x0, y0, z0, a}],
 {x0, 0, 2, 1}, {y0, 0, 2, 1}, {z0, 0, 2, 1}, {a, 1, 5, 2}],
{x0, y0, z0, a}]; // AbsoluteTiming

(* Out: {164.881392, Null} *)

Runs in 165s on my machine. If I want to take the resolution of all 4 parameters down by the factor 10, it would take 10000 times as long - only three weeks.

I do think that it's hard to beat NIntegrate. In addtion to automatically dealing with 3D improper integrals, here's a couple of integrals that NIntegrate just does:

NIntegrate[x^x^x, {x, 0, 6}]

(* Out: 1.102664999366306*10^36300 *)

NIntegrate[Sin[1/x]/x, {x, 0, 1}]

(* Out: 0.624713 *)

What other systems gets those so easily?