I have a code to make about the 2D schrodinger equation time dependant. I was able to figure how to use the finites differences and I guess I was able to solve the 1D schrodinger equation. But know I have to plot my wave function in spherical coordinate (so in 2D?) and create an animation to observe a wave packet at certain energies (3.48 Hartree). It is suggested that E=k^2/2 So I used this to make my schrodinger equation energy dependant. The potential and deltak are given. I also need to make the wave packet propagatinf from r=15 to r=0, I think that I just need to insert a minus in my exponential in my psi0 which I did in the code below. Any idea or even explanation would be welcome
import numpy as np
from scipy import sparse
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from matplotlib.animation import FuncAnimation
from matplotlib import animation
from IPython import display
dx = 0.02 # spatial separation
a= 15
x = np.arange(0.1, a, dx) # spatial grid points
deltak = 0.2 # center of initial gaussian wave-packet
E=3.48
x0= 2
A=1.0 / (deltak * np.sqrt(np.pi)) # normalization constant
hbar = 1
dt = 0.1 # time interval for snapshots
t0 = 0.0 # initial time
tf = 1.0 # final time
t_eval = np.arange(t0, tf, dt)
k = np.sqrt(2 * E)
Vx= 7.5 * x**2 * np.exp(-x)
def psi_t2(t,psi):
return (-1j * (((-0.5 * D2.dot(psi)) + V * psi)))
# Initial Wavefunction
psi0_ = np.sqrt(A) * np.exp(-(x-x0)**2 / (2.0 * deltak**2)) * np.exp(1j * -k * x)
# Solve the Initial Value Problem
sol_ = integrate.solve_ivp(psi_t2, t_span = [t0, tf], y0 = psi0_, t_eval = t_eval,method="RK23")
for i,t in enumerate (sol_.t):
plt.plot(x, np.abs(sol_.y[:,i])**2)
fig=plt.figure()
ax=plt.subplot(1,1,1)
#fig,ax = plt.subplots(figsize=(8,4))
ax.set_xlim(-5, 16)
ax.set_ylim(0, 5)
title = ax.set_title('')
line1, = ax.plot([], [], "k--")
line2, = ax.plot([], [])
def init():
line1.set_data(x, V)
return line1,
def animate(i):
line2.set_data(x, np.abs(sol_.y[:,i])**2)
#line2.set_data(x, np.real(sol.y[:, i]))
#line2.set_data(x, np.abs(psi))
title.set_text('Time = {0:1.3f}'.format(sol_.t[i])) #permet d'afficher le temps
return line1,
anim = animation.FuncAnimation(fig, animate, init_func=init,frames=len(sol_.t), interval=50, blit=True)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close()
