# Wavepackets as solutions of PDEs

I'm currently doing a mathematica for physicists course and am struggling to solve a problem we were given!

I'm supposed to define an initial wavefunction in a harmonic oscillator (displaced, stretched or compressed ground state for example). Then I should propagate the wave function in time and use NDSolve to solve the resulting PDE. To keep the expressions simple, constants are set to 1.

Now I really don't know how to use mathematica well yet so I'm struggling to start this exercise. I first define an initial wavefunction.

ψinit[x_] := Exp[-x^2/2]


Then to propagate in time I can multiply it by the time evolution:

ψinit[x_, t_] := ψinit[x] Exp[-I*E_n*t/ħ]


But I do not understand how I'm supposed to get a PDE from this that I can solve.

Any help appreciated :) Thanks

• You need to use the time-dependent Schrödinger equation (which is a PDE), not the time-independent one (which is an eigenvalue equation). What you wrote is the time propagation of an eigenstate of the time-independent Schrödinger equation, which will never get you any interesting dynamics. Jan 26, 2019 at 21:37

A time-dependent wavefunction ψ[x,t] must satisfy the time-dependent Schrödinger equation $$i\hbar\frac{\partial\psi(x,t)}{\partial t}=H\psi(x,t)$$. For a harmonic oscillator with Hamiltonian $$H=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+\frac12m\omega^2x^2$$ this would be

I*ħ*D[ψ[x,t],t] == -ħ^2/(2m)*D[ψ[x,t],{x,2}] + 1/2*m*ω^2*x^2*ψ[x,t]


Here's an example, using a finite calculation box and setting $$\hbar=m=\omega=1$$:

a = 5;
NDSolve[{I*D[ψ[x,t],t] == -1/2*D[ψ[x,t],{x,2}] + 1/2*x^2*ψ[x,t],
ψ[x,0]==Exp[-x^2/2], ψ[-a,t]==ψ[a,t]==0}, ψ, {x,-a,a}, {t,0,10}]


(not being very careful with boundary conditions here).

To be more careful with the boundary conditions, let's use a beta distribution as a starting wavefunction: it is almost indistinguishable from a normal distribution, with the difference that it reaches zero precisely at the boundaries of the simulation box. For a Beta-Gaussian centered at $$x=0$$ with standard deviation $$\sigma^2$$ this would be

$$\psi_0(x) \propto \left[1-\left(\frac{x}{a}\right)^2\right]^{\frac{1}{4} \left[\left(\frac{a}{\sigma}\right)^2-3\right]}$$

a = 5;
σ = 0.5;
ψ0[x_] = (1 - (x/a)^2)^(((a/σ)^2 - 3)/4);
NDSolve[{I*D[ψ[x, t], t] == -1/2*D[ψ[x, t], {x, 2}] + 1/2*x^2*ψ[x, t],
ψ[x, 0] == ψ0[x], ψ[-a, t] == ψ[a, t] == 0}, ψ, {x, -a, a}, {t, 0, 10}]