# Testing the time dependent Schrodinger Equation with an analytical solution?

I am numerically solving the Schrodinger Equation in 1D first and in higher dimension later, but I want to know the convergence rate of my numerical solver in different grid size and numerical methods.

To calculate the error (or convergence rate), I will need an analytical solution to compare with, but I don't seem to find one for TDSE, even in 1D when potential is 0.

For example, numerically solved Poisson equation can be tested with an analytical solution sin(x). Is there a simple analytical solution for 1D TDSE, or in higher dimension?

Thank you

• You can manufacture a solution. Assume some non-trivial solution and find the corresponding potential by putting your assumed solution into the Schrodinger equation. – cfdlab Feb 28 at 3:07
• @cfdlab generally using manufactured solutions is a bad idea I believe, specially you will stuck at finding proper boundary conditions corresponding to your manufactured solution. – Alone Programmer Feb 28 at 13:31
• @AloneProgrammer The manufactured solution already gives you the boundary condition. Combine that with the knowledge of what kind of bc lead to a well-posed problem for your particular problem. – cfdlab Feb 28 at 16:08
• Schr\"odinger himself found an analytic solutions to the time dependent Schr\"odinger equation which was not separable. This is given in David J. Griffiths book, which sadly I do not have on hand. – user14717 Mar 1 at 22:46

It depends what you consider as an analytical solution. If you consider a Fourier series as an analytical solution, you have this:

$$i \hbar \frac{\partial \Psi (\mathbf{r},t)}{\partial t} = -\frac{\hbar^{2}}{2m} \nabla^{2} \Psi(\mathbf{r},t)$$

Let's say you can use separation of variables and assume this for your wave function as:

$$\Psi(\mathbf{r},t) = F(\mathbf{r}) T(t)$$

So:

$$i \hbar F T^{'} = -\frac{\hbar^{2}}{2m} T \nabla^{2} F$$

or:

$$i \frac{2 m}{\hbar} \frac{T^{'}}{T} = -\frac{1}{F} \nabla^{2} F = \lambda^{2}$$

You have this for temporal part:

$$T(t) = \exp(-i \frac{\hbar}{2m} \lambda^{2} t)$$

For spatial part:

$$\nabla^{2}F + \lambda^{2} F = 0$$

The spatial part equation is just a Helmholtz equation. Now, it depends what geometry you have here. If it is just a cube in three-dimensional space:

$$F(\mathbf{r}) = X(x)Y(y)Z(z)$$

You have:

$$X^{''}YZ + XY^{''}Z + XYZ^{''} + \lambda^{2} XYZ = 0$$

or:

$$\frac{X^{''}}{X} + \frac{Y^{''}}{Y} + \frac{Z^{''}}{Z} + \lambda^{2} = 0$$

or:

$$\frac{X^{''}}{X} = - \frac{Y^{''}}{Y} - \frac{Z^{''}}{Z} - \lambda^{2} = -k_{x}^{2}$$

So:

$$X^{''} + k_{x}^{2} X = 0$$

$$Y^{''} + k_{y}^{2} Y = 0$$

$$Z^{''} + k_{z}^{2} Z = 0$$

and:

$$k_{x}^{2} + k_{y}^{2} + k_{z}^{2} = \lambda^{2}$$

$$X(x) = A_{x} \cos{(k_{x}x)} + B_{x} \sin{(k_{x}x)}$$

$$Y(y) = A_{y} \cos{(k_{y}y)} + B_{y} \sin{(k_{y}y)}$$

$$Z(z) = A_{z} \cos{(k_{z}z)} + B_{z} \sin{(k_{z}z)}$$

If you have Dirichlet boundary condition of zero at all boundaries of cube, you have:

$$X(x) = B_{x} \sin{(\frac{n\pi}{L_{x}}x)}$$

$$Y(y) = B_{y} \sin{(\frac{m\pi}{L_{y}}y)}$$

$$Z(z) = B_{z} \sin{(\frac{l\pi}{L_{z}}z)}$$

Where $$L_{x}$$, $$L_{y}$$, and $$L_{z}$$ are the length of the cube, $$n$$, $$m$$, and $$l$$ are integers, so:

$$k_{x} = \frac{n \pi}{L_{x}}$$

$$k_{y} = \frac{m \pi}{L_{y}}$$

$$k_{z} = \frac{l \pi}{L_{z}}$$

$$\lambda^{2}_{n,m,l} = \Bigg(\frac{n \pi}{L_{x}}\Bigg)^{2} + \Bigg(\frac{m \pi}{L_{y}}\Bigg)^{2} + \Bigg(\frac{l \pi}{L_{z}}\Bigg)^{2}$$

Finally:

$$\Psi_{n,m,l}(x,y,z,t) = B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i \frac{\hbar}{2m} \lambda^{2}_{n,m,l} t)$$

and

$$\Psi(x,y,z,t) = \sum_{n}\sum_{m}\sum_{l} B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i \frac{\hbar}{2m} \lambda^{2}_{n,m,l} t)$$

If for initial condition: $$\Psi(x,y,z,0) = \Psi_{0}(x,y,z)$$

$$B_{n,m,l} = \frac{1}{L_{x}L_{y}L_{z}}\int_{0}^{L_{x}} \int_{0}^{L_{y}} \int_{0}^{L_{z}} \Psi_{0}(x,y,z) \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} dx dy dz$$

Note that the frequency of your system is defined as:

$$\omega = \frac{\hbar}{2m} \mathbf{k} \cdot \mathbf{k}$$

Where:

$$\mathbf{k} = (k_{x},k_{y},k_{z})$$

The energy of the system is defined as:

$$E = \hbar \omega = \frac{\hbar^{2}}{2m} \mathbf{k} \cdot \mathbf{k}$$

So rewriting the solution:

$$\Psi(x,y,z,t) = \sum_{n}\sum_{m}\sum_{l} B_{n,m,l} \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} \exp(-i\omega_{n,m,l} t)$$

Where:

$$\omega_{n,m,l} = \frac{\hbar}{2m} \Bigg( \Big ( \frac{n \pi}{L_{x}} \Big )^{2} + \Big ( \frac{m \pi}{L_{y}} \Big )^{2} + \Big ( \frac{l \pi}{L_{z}} \Big )^{2} \Bigg)$$

So, now you have an analytical solution by means of Fourier series to examine the error and convergence rate your numerical implementation.

Update: For an initial Gaussian wave packet of the form:

$$\Psi(x,y,z,0) = \Psi_{0}(x,y,z) = \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{\sigma_{0}^{2} \pi}} \exp{(i \mathbf{k} \cdot \mathbf{r} - \frac{|\mathbf{r} - \mathbf{r}_{0}|^{2}}{2 \sigma_{0}^{2}})}$$

Where: $$\mathbf{k} = (k_{x}, k_{y}, k_{z})$$, $$\mathbf{r} = (x, y, z)$$, $$\mathbf{r}_{0} = (\frac{L_{x}}{2}, \frac{L_{y}}{2}, \frac{L_{z}}{2})$$.

So:

$$B_{n,m,l} = \frac{1}{L_{x}L_{y}L_{z}}\int_{0}^{L_{x}} \int_{0}^{L_{y}} \int_{0}^{L_{z}} \Psi_{0}(x,y,z) \sin{(\frac{n\pi}{L_{x}}x)} \sin{(\frac{m\pi}{L_{y}}y)} \sin{(\frac{l\pi}{L_{z}}z)} dx dy dz$$

But this 3D integral could be broken down to three 1D integrals of:

$$\mathcal{I}_{x} = \frac{1}{L_{x}} \int_{0}^{L_{x}} \exp{(i k_{x}x - \frac{(x - \frac{L_{x}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{x} x)} dx$$

$$\mathcal{I}_{y} = \frac{1}{L_{y}} \int_{0}^{L_{y}} \exp{(i k_{y}y - \frac{(y - \frac{L_{y}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{y} y)} dy$$

$$\mathcal{I}_{z} = \frac{1}{L_{z}} \int_{0}^{L_{z}} \exp{(i k_{z}z - \frac{(z - \frac{L_{z}}{2})^{2}}{2 \sigma_{0}^{2}})} \sin{(k_{z} z)} dz$$

But:

$$\mathcal{I}_{x} = i\frac{\sigma_{0}}{2L_{x}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{x}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{x}^{2} + i k_{x} L_{x})} (\mathrm{erf} (\frac{- \frac{L_{x}}{2} + i 2 \sigma_{0}^{2} k_{x}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{x}}{2} + i 2 \sigma_{0}^{2} k_{x}}{\sqrt{2} \sigma_{0}})))$$

$$\mathcal{I}_{y} = i\frac{\sigma_{0}}{2L_{y}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{y}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{y}^{2} + i k_{y} L_{y})} (\mathrm{erf} (\frac{- \frac{L_{y}}{2} + i 2 \sigma_{0}^{2} k_{y}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{y}}{2} + i 2 \sigma_{0}^{2} k_{y}}{\sqrt{2} \sigma_{0}})))$$

$$\mathcal{I}_{z} = i\frac{\sigma_{0}}{2L_{z}} \sqrt{\frac{\pi}{2}} (2 \mathrm{erf} (\frac{L_{z}}{2 \sqrt{2} \sigma_{0}}) + \exp{(-2 \sigma_{0}^{2} k_{z}^{2} + i k_{z} L_{z})} (\mathrm{erf} (\frac{- \frac{L_{z}}{2} + i 2 \sigma_{0}^{2} k_{z}}{\sqrt{2} \sigma_{0}}) - \mathrm{erf} (\frac{\frac{L_{z}}{2} + i 2 \sigma_{0}^{2} k_{z}}{\sqrt{2} \sigma_{0}})))$$

So:

$$B(k_{x},k_{y},k_{z}) = \sqrt[\leftroot{-2}\uproot{2}4]{\frac{1}{\sigma_{0}^{2} \pi}} \mathcal{I}_{x} \mathcal{I}_{y} \mathcal{I}_{z}$$

• thank you for the comment. I am now working on heat equation instead, which is simpler but similar to the Schrodinger Eq. Following your appraoch, I could also obtain a solution of fourier series for heat equation, and I think I successfully got it. I would make a detailed answer after I apply your suggested method to Schrodinger Eq. Anyway Thank you! – jchan192 Feb 29 at 8:46

If you have a solution $$\psi$$ to the stationary Schroedinger equation $$H\psi(x) = E \psi(x)$$ then the time dependent Schroedinger equation $$i\hbar \frac{\partial}{\partial t}\Psi(x,t) = H\Psi(x,t)$$ has the solution $$\Psi(x,t) = e^{-iE t/\hbar}\psi(x)$$.

(In other words, every pure eigenstate remains a pure eigenstate -- only the phase rotates with a frequency proportional to $$E$$.)

Take the one-dimensional harmonic oscillator and its (analytical) ground state function $$\Psi_0(x)$$ (which is just a Gaussian) and displace it by a constant length $$a$$. Use the resulting state $$\Psi(x) = \Psi_0(x-a)$$ as initial state to the time-dependent Schrödinger equation. One knows that it should oscillate with a frequency of $$1/2\pi$$ (in atomic units). Evolve the wavepacket for a chosen number of periods, and compare the final state with the initial state (which for an exact propagation method should be identical).

Here is a picture I produced more than a decade ago ($$M$$ in this case is the order of the generalized Crank-Nicolson method). I hope this gets the idea.

AloneProgrammer and Wolfgang provided a general answer, which are very useful for a general initial condition. Here I want to provide a more explicit answer, which is simply the analytical solution of a 1D box.

Let's assume $$V=0$$, and the TDSE has the following form $$i\dfrac{\partial \psi}{\partial t} = -\dfrac{1}{2}\dfrac{\partial^2\psi}{\partial x^2}$$

If the boundary condition is $$\psi(x=0) = \psi(x=1) = 0$$. The time evolution of the real part follows

$$Re(\psi(x,t)) = \sqrt{2}\cos(\frac{n^2 \pi^2}{2}t)\sin(n\pi x)$$

where n is the number of node that you want. You can now compare it with the numerical solution of the TDSE!