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, 2020 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$ Feb 28, 2020 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, 2020 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, 2020 at 22:46

4 Answers 4


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$ Feb 29, 2020 at 8:46
  • $\begingroup$ Sorry for coming back, but I am wondering if similar solution exists for inhomogeneous dirichlet boundary condition? $\endgroup$ Dec 31, 2022 at 3:11
  • $\begingroup$ meaning non-zero dirichlet boundary $\endgroup$ Dec 31, 2022 at 8:47
  • $\begingroup$ @WhatsupAndThanks You can homogenize your Dirichlet BCs to use the same solution above. Note that if $\psi$ is the solution of the Schrodinger equation, then any function in the form of $\phi = \psi + ax + by + cz + d$ is also a solution of the Schrodinger equation, due to that $\nabla^{2} \phi = \nabla^{2} \psi$ and $\partial_{t} \phi = \partial_{t} \psi$. $\endgroup$ Jan 3, 2023 at 18:33
  • $\begingroup$ Since the homogenous case is basicly the particle in a box model, is it right that even if we adopt the transformation 𝜙=𝜓+𝑎𝑥+𝑏𝑦+𝑐𝑧+𝑑, an initial moving gaussian wave will not behave correctly at the boundary? Im asking because I have used the transformation, but there are numerical oscillation/reflection when the moving gaussian wave touched the boundary. however, we expect the wave simply disappears as if there is no boundary. $\endgroup$ Jan 4, 2023 at 8:12

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 and acknowledge you have read our privacy policy.

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