Central Difference Initialization
When we first implemented the central difference method, we pretty much just hand waved away the initialization of buffers for and by using the same values in both buffers.
Let's see what happens if we use a time evolving wave packet to setup and with more realistic values.
The original expression we used for a Gaussian wave packet centered around with momentum was
We can also find several expressions for the time dependent free particle wave function. I chose one from a trusted source. where and .
This expression is perfectly fine for most analysis and calculations you encounter in physics. However, it does not work well with our numerical approach. We need a clear separation between the real and imaginary parts of the wave function.
Let's start with a bit of tuning this expression. He uses , where we have used , we also took and . So then . where and .
We should be able to convince ourselves that this coincides with the original expression for a Gaussian wave packet at .
We need to reorganize the exponent to remove the complex term from the denominator.
Putting this back into our Gaussian wave packet
We want the Gaussian to begin at some arbitrary position so make the transformation
We see some common expressions and
And use Euler's formula to write
Substituting these into the Gaussian wave packet gives an expression quite amenable to a numerical treatment. Similar massaging of mathematical expressions to make them more tractable is not at all uncommon in numerical analysis.
We haven't introduced any new variables in this version of the initializer, so the implementation setup remains the same. There are only some small changes to the computePsi method.
/**
* Compute a time evolving Gaussian wave packet at x for a specific point in time.
*
* globalID An unsigned int giving the thread id for this invocation.
* t A float giving the current time.
*/
fn computePsi(globalID: u32, t: f32) -> vec2f
{
let x = f32(globalID) * parameters.length / f32(parameters.xResolution-1);
let w2 = wavePacketParameters.w * wavePacketParameters.w;
// In our units, hbar = 1
let p0 = wavePacketParameters.k;
let deltaX = x - wavePacketParameters.x0;
let theta = atan(2.0*t/w2)/2.0;
let phi = -theta-p0*p0*t/2.0;
let a = w2*w2 + 4.0*t*t;
let b = (deltaX-p0*t)*(deltaX-p0*t);
let c = phi + p0*deltaX + 2.0*t*b/a;
let phase = vec2(cos(c), sin(c));
let magnitude = pow(2.0*w2/(PI*a), 0.25) * exp(-b*w2/a);
return magnitude*phase;
}
Then when we invoke it, we use and .
waveFunction0[index] = computePsi(index, 0.0);
waveFunction1[index] = computePsi(index, parameters.dt);
We see that the improved initialization has little effect on the overall simulation. Luckily we have grander plans for this code.
