From 4f5ab5e8ee04ab4fd9135b3eaa249d49d627f0d9 Mon Sep 17 00:00:00 2001 From: orosmatthew Date: Wed, 24 Apr 2024 12:05:43 -0400 Subject: [PATCH] Numerov animation --- numerov-2-ani.py | 66 ++++++++++++++++++++++++++++++++++++++++++++ numerov-ani.py | 64 ++++++++++++++++++++++++++++++++++++++++++ numerov-animation.py | 65 +++++++++++++++++++++++++++++++++++++++++++ numerov-thing.py | 60 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 255 insertions(+) create mode 100644 numerov-2-ani.py create mode 100644 numerov-ani.py create mode 100644 numerov-animation.py create mode 100644 numerov-thing.py diff --git a/numerov-2-ani.py b/numerov-2-ani.py new file mode 100644 index 0000000..566a091 --- /dev/null +++ b/numerov-2-ani.py @@ -0,0 +1,66 @@ +import pylab as pl +import numpy as np +import matplotlib.animation as animation + +iterations = 60000 # iterations for approximation + +step = 0.0001 # step size for x +step_sqrd = pow(step, 2) # square of step size + +n = 2 +epsilon = (n * np.pi) ** 2 / (12 ** 2) # 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 = 0 # epsilon + pos - 2 * step # k_0, potential energy at two steps before current +potential_past_1 = 0 # 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 # 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 + +fig, ax = pl.subplots() +ax.set_xlim(-6, 6) +ax.set_ylim(-1, 1) +line, = ax.plot(x_out, y_out, label=f'epsilon = {epsilon}') + + +def update(frame): + global pos, potential_past_1, psi_past_1, potential_past_2, psi_past_2, ax + # Numerov integration loop + # count += 1 + for i in range(0, 1000, 1): + pos += step + potential = epsilon + # potential = epsilon # 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 + + ax.set_xlim(min(x_out), max(x_out)) + ax.set_ylim(min(y_out), max(y_out)) + line.set_data(x_out, y_out) + return line, + + +# Plot +ani = animation.FuncAnimation(fig, update, frames=np.arange((-1 * iterations + 2) // 1000, (iterations + 2) // 1000), + blit=True, repeat=False) +pl.show() diff --git a/numerov-ani.py b/numerov-ani.py new file mode 100644 index 0000000..6049298 --- /dev/null +++ b/numerov-ani.py @@ -0,0 +1,64 @@ +import pylab as pl +import numpy as np +import matplotlib.animation as animation + +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 + +fig, ax = pl.subplots() +ax.set_xlim(-6, 6) +ax.set_ylim(-1, 1) +line, = ax.plot(x_out, y_out, label=f'epsilon = {epsilon}') + + +def update(frame): + global pos, potential_past_1, psi_past_1, potential_past_2, psi_past_2, ax + # Numerov integration loop + # count += 1 + for i in range(0, 1000, 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 + + ax.set_xlim(min(x_out), max(x_out)) + ax.set_ylim(min(y_out), max(y_out)) + line.set_data(x_out, y_out) + return line, + + +# Plot +ani = animation.FuncAnimation(fig, update, frames=np.arange((-1 * iterations + 2) // 1000, (iterations - 2) // 1000), + blit=True, repeat=False) +pl.show() diff --git a/numerov-animation.py b/numerov-animation.py new file mode 100644 index 0000000..c0046b8 --- /dev/null +++ b/numerov-animation.py @@ -0,0 +1,65 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation + +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 + +fig, ax = plt.subplots() +line, = ax.plot([], [], label=f'epsilon = {epsilon}') + +ax.set_xlim(-1, 1) +ax.set_ylim(-1, 1) +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.set_title("Schrodinger Eqn in Harmonic Potential") +ax.legend(loc=1) + + +def init(): + line.set_data([], []) + return line, + + +def update(frame): + global pos, potential, psi, potential_past_2, potential_past_1, psi_past_2, psi_past_1 + for i in range(frame * 1000, (frame * 1000) + 1000): + print(i) + 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 + + line.set_data(x_out, y_out) + return line, + + +ani = animation.FuncAnimation(fig, update, frames=np.arange(0, (iterations - 2)), init_func=init, blit=True) +plt.show() diff --git a/numerov-thing.py b/numerov-thing.py new file mode 100644 index 0000000..2556329 --- /dev/null +++ b/numerov-thing.py @@ -0,0 +1,60 @@ +import pylab as lab +import numpy as np + +iterations = 60000 # iterations for approximation + +step = 0.0001 # step size for x +step_sqrd = pow(step, 2) # square of step size + +well_width = 2 + +# try values 1, 2, and 3 for various solutions +n = 1 +epsilon = (n * np.pi) ** 2 / (well_width ** 2) # 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 + # close enough to infinity :) + potential = 100000000 + if abs(pos) < well_width / 2: + potential = epsilon # 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()