创建一维数组，其中元素是存储函数的二维 numpy 数组的总和

Posted 2023-03-27

【中文标题】创建一维数组，其中元素是存储函数的二维 numpy 数组的总和

【英文标题】：Creating 1D array where elements are the sum of a 2D numpy array of stored functions

【发布时间】：2018-11-12 22:57:00

【问题描述】：

如果之前有人问过这个问题，我深表歉意。我是 python 和一般编程的新手，所以如果有人问我，请指出正确的方向。我正在使用 python 3.7。

我有一个 2D numpy 数组，其中每个元素都是一个存储函数。我想在每一列中添加函数以获得一维数组，其中一维数组的元素是单个函数。我不确定为什么 np.sum() 函数不起作用。我得到一个一维数组，但函数仅来自“npwavefxns”数组的第一列。

即

[[X₀₀,X₀₁,...,X_0n]

[X₁₀, X₁₁,...,X_1n]

...

[X_n0,X_n1,...[X_nn]]

-> [[X₀₀+X₁₀+...+X_n0, X₀₁ +X₁₁+...+X_n1, X_0n+X_1n+...+X _nn]]

np.sum() 函数似乎适用于整数，所以我不确定为什么当元素是函数时它不起作用。下面给出了我的代码示例。如果此代码正常工作，我怀疑在使用“4”个基函数时会得到这 4 个图。

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(1,NB+1):
    clist=[]
    for k in range(1,NB+1):
        F3 = (lambda x: ((EVec.item(k-1, j-1))*
        (np.sin((((k+1)*np.pi)/WL)*x))))
        clist.append(F3)
    wavefxns.append(clist)
npwavefxns=np.asarray(wavefxns)

EQS=[]
for j in range(0,NB):
    F4=np.sum(npwavefxns.item(j))
    EQS.append(F4)
npEQS=np.asarray(EQS)

for j in range(0,NB):
    wfxn1=(lambda x: ((npEQS.item(j))(x)))
    plt.xlabel("Box length")
    plt.ylabel("energy")
    x = np.linspace(0,WL,500)
    plt.plot(x, wfxn1(x), '--m')
    plt.show()

【问题讨论】：

所以事情是这样的：numpy 数组最适合用于存储数字，也许还有字符串。数组足够灵活，您可以存储函数（和其他任意对象），但一般来说，您可能不应该这样做。例如，你不能对函数数组求和：通常你会得到像TypeError: unsupported operand type(s) for +: 'function' and 'function' 这样的错误。您的代码没有出现该错误的原因是 F4=np.sum(npwavefxns.item(j)) 行仅获取单个函数的总和，而该函数只会返回该函数不变。

啊，我明白了。如果我将它们存储在嵌套列表中（这是正确的术语吗？）我仍然会遇到不支持的操作数错误吗？

你会的，因为在他们评估的某个时刻sum([f1, f2, f2]) 和np.array([f1, f2, f3]).sum() 最终都会调用基本的加法运算符。这相当于f1 + f2，函数不支持它。我不得不说，你对 lambda 函数总和的惰性求值的想法很有趣。它不是您可以在 Python 中轻松实现的。您最好的选择是在执行过程中评估函数并仅存储值。或者，您可以尝试调整您拥有的内容，以便最终将功能和输入结合在一起。

【参考方案1】：

好的，据我所知，您遇到了两个问题。其中之一是您试图接管一系列函数（以及与npwavefxns 相关的所有其他内容）的总和。

另一个问题是一个真正的问题，结果导致回到您分配给F3 的lambda。简短的版本是您在lambda 中使用了循环变量j 和k。 lambdas “捕获”变量，以便稍后在调用 lambda 时使用它们。问题是这些变量的值可以改变，就像循环的每次迭代中j 和k 的值一样。当您实际调用这些 lambdas 时，它们最终都使用了完全相同的值 j 和 k（在这种情况下，这是它们在最后一个循环中的最终值）。

我使用称为闭包的技术修复了lambda 问题（请参阅this thread for more details）。简短的解释是，它允许您显式捕获变量的当前值以供以后使用。

无论如何，这是您的代码的完整工作示例。 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）产生输出：

我调整了绘图例程，以便所有波函数都绘制到一个图形上（所以我只有一个图形可以复制/粘贴到这个答案中）。

【讨论】：

谢谢！我想一个更简单的方法可能是绘制数组的主对角线，而不是将所有函数加在一起，只要非对角线元素没有大的扰动，这应该给出相同的结果。该程序的最终目标是计算初始波函数乘以最终波函数的积分，这将为每一行提供一个值。正如您在上一条评论中所建议的那样，我可能会尝试将其构建为一个循环。我只是希望能够在途中绘制和检查功能。