Comment numerov

This commit is contained in:
orosmatthew 2024-04-22 11:52:38 -04:00
parent af9d25ffe8
commit ab94dcc39a

62
numerov.py Normal file
View File

@ -0,0 +1,62 @@
"""
--- Title: Schrodinger Equation in Harmonic Potential
--- Author: Matt Evans (mtdevans.com)
--- Copyright: Do what you want with it.
This program basically just draws the functions for a
given value of epsilon. Epsilon should be an integer n+1/2
for good solutions.
"""
import pylab as lab
iterations = 60000 # iterations for approximation
step = 0.0001 # step size for x
step_sqrd = pow(step, 2) # square of step size
epsilon = 2.5 # energy level, should be integer n+1/2 for good solutions
psi = 0.0 # initial value of wave function
potential = 0.0 # initial value of the potential energy function
pos = -1 * (iterations - 2) * step # initial value of the position
potential_past_2 = epsilon + pos - 2 * step # k_0, potential energy at two steps before current
potential_past_1 = epsilon + pos - step # k_1, potential energy at one step before current
amplitude = 0.1 # initial amplitude of wave function
psi_past_2 = 0 # y_0, wave function at two steps before current
psi_past_1 = amplitude # y_1, wave function at one step before current
x_out = [] # list to store x values for plotting
y_out = [] # list to store y values for plotting
count = -1 * iterations + 2 # counter for the loop
# Numerov integration loop
while count < iterations - 2:
count += 1
pos += step
potential = 2 * epsilon - pow(pos, 2) # potential energy function
b = step_sqrd / 12 # constant used for Numerov
# Numerov method to calculate psi at current step
psi = ((2 * (1 - 5 * b * potential_past_1) * psi_past_1 - (1 + b * potential_past_2) * psi_past_2)
/ (1 + b * potential))
# Save for plotting
x_out.append(pos)
y_out.append(psi)
# Shift for next iteration
psi_past_2 = psi_past_1
psi_past_1 = psi
potential_past_2 = potential_past_1
potential_past_1 = potential
# Plot
lab.figure(1)
lab.plot(x_out, y_out, label=f'epsilon = {epsilon}')
lab.xlabel("x")
lab.ylabel("y")
lab.title("Schrodinger Eqn in Harmonic Potential")
lab.legend(loc=1)
lab.show()