A numerical solution of the time-dependent Schroedinger equation in one dimension for a particle confined to a region 0 £ x £ L.

Reference: S. E. Koonin, Computational Physics, Addison-Wesley Publishing Company Inc, (1986).

The program n1.txt is a simple program written in QuickBasic.  It demonstrates a method for solving the time-dependent Schroedinger equation in one dimension.  The evolution of the wave packet is represented graphically.  The wave packet is observed to spread with time.  You are encouraged to modify the program, add a user interface or translate it to another programming language.

Explanation:

We want to solve the Schroedinger equation

,

or

.

The particle is confined to the region 0 £ x £ L so we need

.

Let x_i denote points on a grid in space and tⁿ denote points on a grid in time.

Define y_iⁿ=y(x_i,tⁿ).  To solve the Schroedinger equation numerically, we convert it into a difference equation.

,

,

.

Here we have used

(1+A)^-1(1-A)=(1+A)^-1-(1+A)^-1A+(1+A)^-1-(1+A)^-1=2(1+A)^-1-(1+A)^-1(1+A)=2(1+A)^-1-1.

Let us now look at a Taylor series expansion.

.

This yields

.

We can therefore approximate the second derivative of the wave function with respect to x in the following way:

.

This yields the following expression for the Hamiltonian:

.

If we choose to represent y by a column matrix of its values y_i at the grid points, then H will be represented by a tridiagonal matrix.  To illustrate, let f =Hy .  Then

f _i=H_i,i-1y _i-1+H_i,iy _i+H_i,i+1y _i+1.

,

.

To solve for the time evolution of a wave packet, we need to know the shape of the wave packet at t=0.  Assume y_i⁰ is given with y₀⁰=y_N⁰=0.  How do we calculate the shape of the wave packet at a later time?

Let us define the function c_iⁿ at the grid points through

.

Then

.

.

We can write out this equation explicitly using 

.

.

.

This equation is of the form  

with A_i^- = A_i⁺ = 1.

1. Assume that the solution is of the form  .
2. Then  ,
3. or   with .

Comparing expressions (1.) and (3.) we find

.

We can now determine g_iⁿ and b_iⁿ at each grid point.

The boundary conditions are y_Nⁿ=y_Nⁿ⁺¹=0.  Therefore   .

can be satisfied by setting aⁿ_N-1=gⁿ_N=0.  Then bⁿ_N-1=0.

We can now calculate gⁿ_N-1, gⁿ_N-3, gⁿ_N-3, … recursively using .

.

Once we have calculated g_iⁿ each grid point we calculate the b_iⁿ recursively using .

Re(b_i-1ⁿ)=Re(g_iⁿ)Re(b_iⁿ-b_iⁿ)-Im(g_iⁿ)Im(b_iⁿ- b_iⁿ)

Im(b_i-1ⁿ)=Re(g_iⁿ)Im(b_iⁿ-b_iⁿ)+Im(g_iⁿ)Re(b_iⁿ- b_iⁿ)

Now that we have calculated all the b_iⁿ and g_iⁿ we can calculate the wave function at a later time.

Starting with y₀ⁿ=c₀ⁿ=0 we calculate c₁ⁿ, y₁ⁿ⁺¹, c₂ⁿ, y₂ⁿ⁺¹, …, recursively.

cⁿ_i+1=gⁿ_i+1c_iⁿ+b_iⁿ , y_iⁿ⁺¹=c_iⁿ-y_iⁿ.

Re(cⁿ_i+1)=Re(gⁿ_i+1)Re(c_iⁿ)-Im(gⁿ_i+1)Im(c_iⁿ)+Re(b_iⁿ)

Im(cⁿ_i+1)=Re(gⁿ_i+1)Im(c_iⁿ)+Im(gⁿ_i+1ⁿ)Re(c_iⁿ) )+Im(b_iⁿ)

Let the particle be an electron.  Let us measure all distances in units of 10^-10m and all times in units of 10^-16s.  Let V_i=0 for all i.

Then A_i⁰=-2+i3.457(Dx)²/Dt, b_iⁿ=i2Im(A_i⁰)y_iⁿ.