Time Evaluation of WavePacket: Program 1

Consider a one-dimensional quantum system described by the time-dependent Schrödinger equation.

\[ i\hbar \frac{\partial \psi(x,t)}{\partial t} = \left( -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right) \psi(x,t) \]

Time evolution of a wave packet moving in free space. A Gaussian wave packet is defined as:

$$\psi(x, t=0) = A \cdot e^{-\frac{(x - x_0)^2}{2 \sigma^2}} \cdot e^{i k_0 x}$$

The energy of the wave packet is defined as:

$$E = \frac{\hbar^2 k_0^2}{2m} + \frac{\hbar^2}{2m \sigma^2}$$

Analyze Time evolution of the wave packet by Crank-Nicolson Method

'''Time evolution of Wavepacket'''
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import simps

# Parameters
L, x0, k, a, N = 40.0, 4.0, 5.0, 2.0, 120
xs, h = np.linspace(0, L, N, retstep=True)
us = np.exp(-(xs-x0)**2/a) * np.exp(1.j*k*xs)
Vs = xs[1:-1]*0.0
dt = 0.1
mu = (0. + 1.j)*dt/(2*h**2)
t, tf = 0, 2.0

# Matrix for evolution
A1 = np.eye(N-2)*(1.0 + 2.0*mu + 1.j*dt/2*Vs) - mu*np.eye(N-2, k=1) - mu*np.eye(N-2, k=-1)
A2 = np.eye(N-2)*(1.0 - 2.0*mu + 1.j*dt/2*Vs) + mu*np.eye(N-2, k=1) + mu*np.eye(N-2, k=-1)
B = np.dot(np.linalg.inv(A1), A2)

# Set up the figure and axis
fig, ax = plt.subplots()
line1, = ax.plot([], [], label='Real part')
line2, = ax.plot([], [], label='Probability density')

def init():
    ax.set_xlim(xs[0], xs[-1])
    ax.set_ylim(-0.6, 0.6)
    ax.legend()
    return line1, line2

def update(frame):
    global us, t
    us[1:-1] = np.dot(B, us[1:-1])
    t += dt
    normus = us / np.sqrt(simps(us * np.conj(us), xs))
    line1.set_data(xs, np.real(normus))
    line2.set_data(xs, abs(normus)**2)
    return line1, line2

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=np.arange(0, int(tf/dt)), init_func=init, blit=True, interval=10)
ani.save('animation3.gif')

TDSE Solution

Time Evaluation of WavePacket: Program 1

Output 1