70 lines
1.9 KiB
Python
70 lines
1.9 KiB
Python
import numpy as np
|
|
import pylab as pl
|
|
|
|
|
|
# recurrence relation is H_n(x) = 2*x*H_{n-1}(x) - 2*(n-1)*H_{n-2}(x)
|
|
# source https://en.wikipedia.org/wiki/Hermite_polynomials
|
|
def hermite(x: np.ndarray, n: int) -> np.ndarray:
|
|
if n == 0:
|
|
return 1.0 + 0.0 * x
|
|
if n == 1:
|
|
return 2.0 * x
|
|
return 2.0 * x * hermite(x, n - 1) - 2.0 * (n - 1) * hermite(x, n - 2)
|
|
|
|
|
|
def normalize(psi: np.ndarray) -> None:
|
|
psi /= np.sqrt(np.sum(psi * np.conjugate(psi)))
|
|
|
|
|
|
x_vals = np.linspace(-5, 5, 1000)
|
|
gaussian = np.exp(-x_vals ** 2)
|
|
|
|
psi_funcs = []
|
|
for n in range(4):
|
|
psi = hermite(x_vals, n) * gaussian
|
|
normalize(psi)
|
|
psi_funcs.append(psi)
|
|
|
|
print("===================")
|
|
print("Check Orthogonality")
|
|
print("===================")
|
|
for i in range(4):
|
|
for j in range(i, 4):
|
|
print(f'psi_{i} * psi_{j} = {np.dot(psi_funcs[i], psi_funcs[j])}')
|
|
|
|
print("==================")
|
|
print("Expectation Values")
|
|
print("==================")
|
|
expect_xs = []
|
|
for i, psi in enumerate(psi_funcs):
|
|
expect_x = np.sum(np.conjugate(psi) * x_vals * psi)
|
|
expect_xs.append(expect_x)
|
|
print(f'<x>_{i} = {expect_x}')
|
|
|
|
expect_x_sqrds = []
|
|
for i, psi in enumerate(psi_funcs):
|
|
expect_x_sqrd = np.sum(np.conjugate(psi) * x_vals ** 2 * psi)
|
|
expect_x_sqrds.append(expect_x_sqrd)
|
|
print(f'<x^2>_{i} = {expect_x_sqrd}')
|
|
|
|
for i, psi in enumerate(psi_funcs):
|
|
sigma_x = np.sqrt(expect_x_sqrds[i] - expect_xs[i] ** 2)
|
|
print(f'sigma_x_{i} = {sigma_x}')
|
|
|
|
pl.rcParams['figure.dpi'] = 300
|
|
fig, axs = pl.subplots(2, 2, tight_layout=True)
|
|
fig.tight_layout(pad=2.0)
|
|
axs[0, 0].plot(x_vals, psi_funcs[0])
|
|
axs[0, 0].set_title('ψ0')
|
|
axs[0, 0].grid()
|
|
axs[0, 1].plot(x_vals, psi_funcs[1], 'tab:orange')
|
|
axs[0, 1].set_title('ψ1')
|
|
axs[0, 1].grid()
|
|
axs[1, 0].plot(x_vals, psi_funcs[2], 'tab:green')
|
|
axs[1, 0].set_title('ψ2')
|
|
axs[1, 0].grid()
|
|
axs[1, 1].plot(x_vals, psi_funcs[3], 'tab:red')
|
|
axs[1, 1].set_title('ψ3')
|
|
axs[1, 1].grid()
|
|
pl.show()
|