Animate the scattering of a Wave Packet

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?

• Could you please post your animation's result as an answer when you're done? – Dr. belisarius Sep 24 '12 at 0:17
• Why do you think M- is "not the best choice for intense numerical integration" - do you have examples where it does not work well? – user21 Sep 24 '12 at 2:15
• @ruebenko Agreed. I tried to address that issue in my response. – Mark McClure Sep 24 '12 at 2:16
• @MarkMcClure, and that's another nice answer from you! – user21 Sep 24 '12 at 2:17

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?

• Wow, I really had no idea I could throw my function right in, and have it integrate in 3 dimensions. I'll try it out! Just a quick question. The integrand is complex-valued. Will that be a problem? – QuantumDot Sep 23 '12 at 23:52
• Complex values should be no problem. – Mark McClure Sep 23 '12 at 23:54