PENG9570: The Schrödinger Equation
Sølve Selstø
March 12, 2026
In this session we will take the heat equation as our starting point and, by slight adjustments, introduce the Schrödinger equations. “Equations” in plural as the Schrödinger equation comes in both a time-dependent and a timeindependent version.
1 The Heat Equation
The heat equation reads
where is the temperature at position ( ) at time , and is the thermal diffusivity, which is a measure on how fast the temperature tends to spread out.
In one dimension it is simpler:
where .
- Download the Python script SolveHeatEq.py and run it - preferibly using Spyder.
- In this script, an initial condition is already implemented. What is this condition; what is the initial temperature distribution ?
- Run your code again for a few choices of inputs - both for and for the numerical ones, i.e., the box size and the temporal and spatial step sizes and . For a given , how small must and be for your solution to make sense?
In our implementation, the temperature distribution at time is approximated by a vector which, ideally, interpolates the actual solution,
where and . The double spatial derivative is implemented by means of the three point finite difference stencil
which allows us to approximate the double derivative as a matrix multiplication on ,
- Now, replace the somewhat primitive forward Euler solution,
with the Crank-Nicolson scheme:
Does the performance improve? 5) What is it about the “big O-terms” in Eqs. (6) and (7) that suggest that the Crank-Nicolson scheme is superior?
2 The Schrödinger Equation for a Particle Moving Freely
In quantum physics, all information there is to obtain for a quantum particle is contained in what is called the wave function, . Much can be said about this function; here we will leave it to this: is the probability density that a position measurement at time on the quantum particle provides the position .
The wave function will in general be time-dependent, and when left undisturbed, it will evolve according to the time-dependent Schrödinger equation:
Here is the potential that the quantum particle is exposed to. In classical physics, this would correspond to the force in Newton’s second law. (We do not use the “F-word” in quantum physics.) is a very small natural constant called (the reduced) Planck constant, and is the mass of the quantum particle.
In case there is no potential, , and we restrict our quantum particle to one single dimension, Eq. (8) becomes much simpler:
16) The wave function is not determined by the Schrödinger equation alone. There are two other circumstances that affects it. Which ones? 7) By comparing Eq. (9) with Eq. (2) we may see that these equations are very similar. How so? And what is the big difference?
Despite this big difference, implementing the solution of Eq. (9) is very similar to that of Sec. 1; we can still solve the partial differential equation using our matrix approximation for the double spatial derivative combined with the Crank-Nicolson method.
Again, we discretize the wave function in the same manner as before:
and write the Schrödinger equation as a coupled ordinary differential equation:
where is a matrix. 8) How must we modify Eq. (7) in order to solve the Schrödinger equation, Eq. (9), instead of the heat equation, Eq. (2)? 9) Now, solve the Schrödinger equation numerically by implementing this modification in your code from Sec. 1. Take your initial condition to be a Gaussian wave:
To simplify things, let us set and ; this is admissible by choosing these constants as the units of their respective quantities. (Correspondingly, we are certainly not talking about metres, seconds and kilograms in our simulations.) Instead of plotting , which is complex, use the square of its absolute value, . Initially, take the parameters which defines the initial wave to be and . 10) What kind of numerics do you now need to see a smooth, reasonable evolution? What happens when you wave hits the boundary of your numerical grid at ? 11) Rerun your simulation with different values for the three parameters , and . What do each these three physical parameters represent? 12) Try and construct a new initial wave function similar to the one in Eq. (12) - except this time, let it be a sum of two waves. These two waves should not overlap initially. They should also travel towards each other. What happenes when the waves meet?
3 The Schrödinger Equation for a Particle Exposed to a Potential
Thus far we have not included any potential in our simulation. We will now introduce it. Specifically, we will take it to be time-independent and have this form:
In our simulation, it may be introduced into Eq. (11) simply by adding a diagonal matrix with the -values on the grid to the matrix . 13) What will the -matrix in Eq. (11), which we call the Hamiltonian, now be? 14) With the same parameters as you used initially - both for the numerical grid and for the initial state, rerun your simulation with a nonzero potential. You can start with the parameter values and . If you manage to make a plot of the potential along with your wave function in your simulation, that would make it more informative. 15) Try and change the height and the width of your barrier. Also, run your simulation with various values for the initial velocity of your quantum particle. How does this affect the dynamcis? Any oddities? What happens if you use a negative value for ?
4 The Schrödinger Equation for a Particle Trapped in a Potential
You may have learned that systems which demonstrate wave-behaviour sometimes form standing waves. A guitar string would be a good example of such a system. Also a quantum system may - under certain conditions - form standing waves. 16) Suppose a wave function has this form:
For such a function to be a solution of the time-dependent Schrödinger equation, which we now will write compactly as
which equation must the time-independent part fulfill ?
2The answer to the above question is called the time-independent Schrödinger equation. It is, in effect, an eigenvalue equation. And the corresponding eigenvalues are the possible outcomes that an energy measurement on the system may produce. Correspondingly, we call the eigenvalues of the Hamiltonian eigenenergies.
If the quantum system is confined somehow, for instance by placing a quantum particle in a potential from which it cannot escape, these eigenenergies are quantized. By this we mean that they will form a discrete set of possible values. The electrons which are stuck to a nucleus in an atom, for instance, are subject to this phenomenon. Again, much could be said about this. However, we will not do so here.
Suppose that our Hamiltonian is time-dependent and that it only has a discrete set of eigenenergies and eigenstates, and , respectively, viz.,
It can be shown that any (admissible) wave function may be written as a linear combination of these eigenstates. Or, in other words, for any give wave function there will always exist a set of coefficients such that
- Suppose our wave function at time has the above form. Now, show that at time , our wave function reads
In many situations the energy states we are most interested in are the one with the lowest energy - the ground state. Chemists and physicists have developed a rich plethora of methods for determining ground states - and the ground state energy - of quantum systems such as atoms, ions and molecules. One of these methods consists in, artificially, replacing the time in the time-dependent Schrödinger equation, Eq. (8), by imaginary time:
- Show that if you make this replacement in Eq. (9), you actually recover Eq. (2). In order to solve this equation, you may simply replace with -it in your Crank-Nicolson implementation.
- For a system with a nozero potential in the Hamiltonian, explain why virtually any initial state becomes proportional to the ground state wave function when it evolves according to the Schrödinger equation with imaginary time.
- With some negative value for in Eq. (13), run your simulation with imaginary time in order to find the ground state. In doing so, the wave function will become extremely large in this case. To avoid this, we may renormalize it at each time step,
- Can you augment your implementation so that it also estimates the ground state energy ? Hint: If and , what is ?
5 Solutions
- From line 45 , uFunk0 - np.sqrt(np.abs(x)), we may see that
It has a “sharp point” at , one that quickly becomes smoother. 3) This equation seems to be rather flexible when it comes to the numerics - at least with this initial condition. For instance, with we get a simulation which seems reasonable with . These are very large values for and . The parameter essentially dictates how fast the diffusion happens; if we increase it, the initial distribution flattens out more quickly. Correspondingly, we also need a lower value for to get reasonable results. 4) Actually, for this case, the forward Euler-method preforms rather well. However, if we plot this solution with the Crank-Nicoloson approximation together - for comparatively large values of and we may see a certain difference. 5) For each time step, the forward Euler method introduces an error which is proportional to . With steps in total, where is the total duration, this adds a total error proportional to in total. In other words, if i reduce by a factor 10, also the total error is reduced by a factor 10. For the Cranck-Nicolson scheme, however, the local error is proportional to and the total, global error is proportional to . So if we reduce by a factor , the total error is reduced to of the original error. 6) To solve the Schrödinger equation, like any differential equation, we need an initial condition. We cannot determine the evolution of the wave function unless we know where to start from. The other aspect is more subtle: The absolute value squared of the wave function, , is a probability distribution telling us how likely it is that a position measurement will provide the result . Thus, it reflects an uncertainty; we do not know where the particle will be measured to be, we only have the probabilityes of the various possible outcomes. Suppose now that a measurement has been made and we found that the particle was positioned at . Then there is no longer any uncertainty about where the particle is. Correspondingly, the probability distribution should become a sharp peak centered at . We say that the wave function collapses when the particle is measured. This collapse is fundamentally random, and it is not governed by the Schrödinger equation. 7) We rewrite Eq. (9) so that we get the time derivative alone on the left hand side:
except for the names versus , which is not really any difference, this is identical to Eq. (2) with . So this is quite similar. But there is one big difference: si real, is imaginary. Correspondingly, the wave function is complex, while the temperature distribution is real. 8) We simply replace in Eq. (7) with :
or, in terms of the matrix :
- In this implementation, there are three things that need to be modified: 1) The above adjustment of the Cranck-Nicoloson matrix, 2) updating the initial condition and 3) making sure that the plot displays the absolute value squared of the solution, not the complex solution directly. Pyplot seems to dislike using the column vector np.abs(Psi)**2 as input in this regard. It could be fixed by a slightly more cumbersome construct: np.power(np.abs(Psi), 2).
- You may find that the Scrhödinger equation requires more precise numerics for a reasonable solution. And if you try to solve it with the forward Euler method, chances are you are bound to struggle. In case your simluation lasts long enough for the wave to hit the boundary at , you will see it bouncing back. This is because we have, in effect, implemented Dirichlet boundary conditions.
or, strictly speaking , when we implemented the numerical derivation as a matrix multiplication. In terms of physics, we can think of this as hard, massive walls placed at . 11) Hopefully, you will find that is the starting position of the wave, is its velocity and is a parameter that determines the width of the initial wave; the smaller is, the wider is . Actually, is the uncertainty of the momentum or velocity of the quantum particle. The joint uncertainty of the position and momentum is restricted by the Heisenberg uncertainty relation
which is fulfilled with equality for our particular choice of initial wave. 12) You could choose to have your initial state to be the sum of two wave of the same form as in Eq. (12) - except that and have opposite signs. The one with negative should have positive and vice versa. This way we will have waves traveling towards each other, and we will see an interference pattern emerge when they overlap - just like ripples in a pond when you throw two stones into it. 13) We simply add a diagonal matrix with the entries to :
Eq. (21) still apply. 15) If you compare the kinetic energy of the initial wave,
where , with the energy that it takes - classically - to get over the barrier, you may be surprised to see that large portions of the wave makes it to the other side. And with a comparatively narrow barrier, i.e., a small parameter, quite a significant part of the way may penetrate the barrier. This phenomenon is called tunneling. It is actually used technologically - for microscopy among other things. Another thing which you would not see with a classical particle, i.e., a particle that follows Newton’s laws rather then the Schrödinger equation, is the fact that parts of it may be reflected when it hits a well - a potential with negative . 16) We insert the separated wave function, into the time-depenent Schrödinger equation
In the last equation we have divided by and reversed the order of the equation. We define the Hamiltonian
and introduce the energy so that it may be written
This is the time-independent Schrödinger equation. It determines which energies are admissible for a quantum particle. Sometimes the possible energies, i.e., the -values that emerge from the Eq. (22), constitute a set of discrete energy levels. This phenomenon is called quantization. 17) We consider first the left hand side of the time-dependent Schrödinger equation, with on the form of Eq. (18):
The right hand side reads
where we have used Eq. (16). If we compare, we may see that the left and right hand sides are identical, and the time-dependent Schröginger equation is fulfilled. 18) If we make the replacement , the “time-derivative” becomes
Inserted into Eq. (9) we get
Since we arrive at
which is identical to Eq. (2) with the real and positive -value . 19) To find the solution of the imaginary-time Schrödinger equation
we can introduce the replacement in Eq. (19) directly into Eq. (18):
For any choice of initial state , there will always be a set of -coefficients such that Eq. (17) holds. In order to see how our state converges towards a state proportional to the ground state , we divide by :
Since for all , all exponential factors falls off to zero exponentially fast - except for the first one. Correspondingly, the only surviving term is the first one:
So, unless we are extremely unlucky with our initial guess and happens to be zero, we are left with a state proportional to our ground state in the end. 21) Suppose we have evolved our odd imaginary time-wave function long enough so that
Here we have assumed that is normalized,
We can disregard the coefficient and the factor after having done the renormalization in Eq. (20). Then the next step will take us to
which, in turn, gives
We find the ground state energy simply by solving this equation for :