Хорошо, насколько я могу судить, вы столкнулись с двумя проблемами.Одним из них была сумма, которую вы пытались взять на себя массив функций (и все остальное, связанное с npwavefxns).
Другая проблема заключалась в том, что она была настоящей вонючей и, как оказалось, привела к тому, что lambda было назначено для F3.Короткая версия заключается в том, что вы использовали переменные цикла j и k в вашем lambda.lambda s "перехватывает" переменные, чтобы вы могли использовать их при последующем вызове lambda.Уловка в том, что значение этих переменных может меняться, как и значения j и k на каждой итерации цикла.К тому времени, когда вы на самом деле вызывали эти lambdas, все они в итоге использовали одно и то же значение j и k (которое в данном случае было окончательным значением, которое они имели в каждом последнем цикле).
Я исправил проблему lambda, используя технику, называемую замыканием (подробнее см. в этой теме ).Краткое объяснение состоит в том, что он позволяет вам явно захватить текущее значение переменной для последующего использования.
Во всяком случае, вот полный рабочий пример вашего кода.Все, что находится над строкой wavefxns=[], осталось нетронутым:
from scipy import mat, matrix, integrate
from scipy.integrate import quad
import numpy as np
from numpy import linalg as la
import os
from matplotlib import pyplot as plt
# Defining variables and functions
MP=float(9.10938356e-31) #mass of electron in kg
WL=float(1e-10) #length of well in meters
CON=float(1.60218e-19) #constant height in joules
Hb = float(1.054571726e-34) #reduced planck's constant in J s
NB=int(input("Number of basis functions ")) #define number of basis sets to be used
#####Potential energy of initial state#####
PE=[]
for j in range(1,NB+1):
alist=[]
for k in range(1,NB+1):
F1=integrate.quad(lambda x:((2/WL)*np.sin((k*np.pi*x)/WL)*
((-CON)*np.sin(np.pi*x/WL))*np.sin((j*np.pi*x)/WL)),0,WL)[0]
if F1 < -1e-25:
F1=F1
elif F1 > 1e-25:
F1=F1
else:
F1=0
alist.append(F1)
PE.append(alist)
PEarray=np.asarray(PE)
#####Kinetic Energy of initial state#####
KE=[]
for j in range(1,NB+1):
blist=[]
for k in range(1,NB+1):
F2=integrate.quad(lambda x:(((((Hb**2)*(j**2)*(np.pi**2))/(MP*(WL**3)))*
((np.sin(j*np.pi*x/WL))*(np.sin(k*np.pi*x/WL))))),0,WL)[0]
if F2 < -1e-25:
F2=F2
elif F2 > 1e-25:
F2=F2
else:
F2=0
blist.append(F2)
KE.append(blist)
KEarray=np.asarray(KE)
#####Adding PE and KE to give the full hamiltonian of the initial state#####
#####Then convert the list to a numpy array#####
sum=[0]*NB
for i in range(NB):
sum[i]=[0]*NB
for y in range(len(PEarray)):
for z in range(len(PEarray[0])):
sum[y][z]=PEarray[y][z]+KEarray[y][z]
npsum=np.asarray(sum)
EVal, EVec=la.eigh(npsum)
wavefxns=[]
for j in range(0,NB):
clist=[]
for k in range(0,NB):
F3 = (lambda a,b: (lambda x: ((EVec.item(b-1, a-1)) * (np.sin((((b+1)*np.pi)/WL)*x)))))(j,k)
clist.append(F3)
wavefxns.append(clist)
gridspec_kw = {'wspace': 0, 'hspace': 0}
fig,axarr = plt.subplots(NB, sharex=True, squeeze=False, gridspec_kw=gridspec_kw, figsize=(3,7.5))
fig.subplots_adjust(left=0, bottom=0, right=1, top=1)
for j,ax in zip(range(0,NB), axarr.ravel()):
wfxn = lambda x: np.sum([wavefxns[i][j](x) for i in range(len(wavefxns))], axis=0)
if j==(NB - 1):
ax.set_xlabel("Box length")
ax.set_ylabel("energy")
x = np.linspace(0,WL,500)
ax.plot(x, wfxn(x), '--m')
, которое (при запуске и вводе 4) выдает:
Я настроил подпрограммы построения графиков так, чтобы все волновые функции отображались на одной фигуре (и поэтому у меня была только одна фигура для копирования / вставки в этот ответ).