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

  • 1
    $\begingroup$ You can manufacture a solution. Assume some non-trivial solution and find the corresponding potential by putting your assumed solution into the Schrodinger equation. $\endgroup$ – cfdlab Feb 28 '20 at 3:07
  • $\begingroup$ @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. $\endgroup$ – Alone Programmer Feb 28 '20 at 13:31
  • 1
    $\begingroup$ @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. $\endgroup$ – cfdlab Feb 28 '20 at 16:08
  • $\begingroup$ 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. $\endgroup$ – user14717 Mar 1 '20 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)$$


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


$$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$$


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


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


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

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

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


$$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}$$


$$\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)$$


$$\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}$$


$$\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)$$


$$\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})$.


$$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$$


$$\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}})))$$


$$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}$$

  • $\begingroup$ 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! $\endgroup$ – jchan192 Feb 29 '20 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!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.