下面是我用来尝试查找能量水平的代码，但是在尝试运行我的最后一段代码时，我一直收到此错误。当我将势函数更改为仅返回 0 的 infintesquare 井函数时，我运行完美，但是当我将势更改为相对于 x 变化的谐波势时，代码不起作用。
我不知道如何修复这段代码。有人可以告诉我它有什么问题以及如何解决它：
m = 9.109383702 * 10**-31 # kg, electron mass
hbar = 1.054571817 * 10**-34 # J.s constant
e = 1.602176634 * 10**-19 ## A⋅s, electron charge
d = 5 * 10**-9 # side length of cubic quantum dot, meters
c = 2.998e8 # speed of light in m/s
xstart = -d/2 # start position,meters
xend = d/2 # end position
N = 2000 # number of points for Runge-Kutta
h = (xend - xstart)/N # size of Runge-Kutta steps
tpoints = np.arange(xstart, xend, h)
def V(x, potential_function):
'''
This function returns the potential energy of a particle in a 1D potential well.
Parameters:
x : float
The position of the particle in meters.
potential_function : callable
A function that takes the position x as an argument and returns the potential energy at that position.
Returns:
float
The potential energy of the particle at the given position x.
'''
return potential_function(x)
def f(r, x, E, potential_function):
psi = r[0]
phi = r[1]
d_psi = phi
d_phi = ((2*m)/hbar**2)*(V(x, potential_function)-E)*psi
return np.array([d_psi, d_phi])
def RungeKutta2d(r, x, f, E, potential_function):
xpoints = []
ypoints = []
for t in tpoints:
xpoints.append(r[0])
ypoints.append(r[1])
k1 = h * f(r, x, E, potential_function)
k2 = h * f(r + 0.5 * k1, x + 0.5 * h, E,potential_function)
k3 = h * f(r + 0.5 * k2, x + 0.5 * h, E, potential_function)
k4 = h * f(r + k3, x + h, E, potential_function)
r = r + (k1 + 2 * k2 + 2 * k3 + k4) / 6
xpoints.append(r[0])
ypoints.append(r[1])
return np.array([xpoints, ypoints])
V0 = 700*e
def harmonic_potential(x):
return V0*(x**2)/((d/2)**2)
fig, ax = plt.subplots(1, 2, figsize=(12, 6), gridspec_kw={'wspace':0.3})
# Find energy states for n = 2, 3, 4, 5
ax[0].axvline(x=-d/2,c='#5f5f5f',ls='-',lw=2.5)
ax[0].axvline(x=d/2,c='#5f5f5f',ls='-',lw=2.5)
ax[1].axvline(x=-d/2,c='#5f5f5f',ls='-',lw=2.5)
ax[1].axvline(x=d/2,c='#5f5f5f',ls='-',lw=2.5)
for n in range(1,4):
# Define initial energy guesses
E_guess = (np.pi**2 * hbar**2 * n**2) / (2*m*d**2)
E1 = E_guess - 3 * 10**-21
E2 = E_guess
# Solve for the first guess
# Here we set the initial conditions in a separate array r
r = np.array([0, ϕ])
psi1 = RungeKutta2d(r, tpoints, f, E1,harmonic_potential)[0, N]
psi2 = RungeKutta2d(r, tpoints, f, E2,harmonic_potential)[0, N]
### now for the secant method to converge on the right answer:
tolerance = e/100000 # set the tolerance for convergence
while abs(E2-E1) > tolerance:
E3 = E2 - psi2*(E2-E1)/(psi2-psi1)
# update initial energies for the next iteration
E1 = E2
E2 = E3
# and recalculate wavefunctions
psi1 = RungeKutta2d(r, tpoints, f, E1,harmonic_potential)[0, N]
psi2 = RungeKutta2d(r, tpoints, f, E2,harmonic_potential)[0, N]
soln = RungeKutta2d(r, tpoints, f, E3,harmonic_potential)
psi = soln[0][:len(tpoints)]
# Normalize wavefunction and plot result
if n % 2 == 0:
I = h*(0.5*psi[0]**2 + np.sum(psi[1:-1]**2) + 0.5*psi[-1]**2)
psi_norm = psi/np.sqrt(I)
ax[0].plot(tpoints, psi_norm, label=f'n = numerical {n}')
ax[0].set_xlabel('Position (m)')
ax[0].set_ylabel('Wavefunction')
ax[0].set_title("Even wavefunctions")
ax[0].legend()
else:
I = h*(0.5*psi[0]**2 + np.sum(psi[1:-1]**2) + 0.5*psi[-1]**2)
psi_norm = psi/np.sqrt(I)
ax[1].plot(tpoints, psi_norm, label=f'n = numerical {n}')
ax[1].set_xlabel('Position (m)')
ax[1].set_ylabel('Wavefunction')
ax[1].set_title("Odd wavefunctions")
ax[1].legend()
print(f"Energy level n = {n}: E = {E3}"
/tmp/ipykernel_692/532686950.py:37: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
return np.array([d_psi, d_phi])
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/tmp/ipykernel_692/532686950.py in <cell line: 68>()
75 r = np.array([0, ϕ])
76
---> 77 psi1 = RungeKutta2d(r, tpoints, f, E1,harmonic_potential)[0, N]
78 psi2 = RungeKutta2d(r, tpoints, f, E2,harmonic_potential)[0, N]
79
/tmp/ipykernel_692/532686950.py in RungeKutta2d(r, x, f, E, potential_function)
45 k1 = h * f(r, x, E, potential_function)
46 k2 = h * f(r + 0.5 * k1, x + 0.5 * h, E,potential_function)
---> 47 k3 = h * f(r + 0.5 * k2, x + 0.5 * h, E, potential_function)
48 k4 = h * f(r + k3, x + h, E, potential_function)
49 r = r + (k1 + 2 * k2 + 2 * k3 + k4) / 6
ValueError: operands could not be broadcast together with shapes (2,) (2,2000)