在numpy中使用集成方案的问题

Question

我正在尝试执行函数的积分，以便找到平均位置。然而，当我使用quad执行积分时，我得到了尺寸不匹配的问题。当我单独运行该函数时，它可以正常工作。然而，当该函数被四重积分方案使用时，它给出了维数不匹配的误差。我将在下面发布完整的功能代码和错误消息，我希望有人能告诉我哪里出了问题，以及如何修复它。

如果有什么不清楚的地方，请让我知道，这样我就可以添加更多的信息。

import numpy as np
import time
from scipy.sparse.linalg import eigsh
from scipy.sparse.linalg import spsolve
from scipy.sparse import diags
from scipy.sparse import identity        
import scipy.sparse as sp        
from scipy.integrate import quad
x0     = 15
v      = 16
sigma2 = 5
tmax   = 4       
        
def my_gaussian(x, mu=x0, var=5):
    
    return np.exp(-((x - mu)**2 / (2*var))+(1j*v*x/2))

L = 100
N = 3200
dx = 1/32
x = np.linspace(0, L, N)
func = lambda x: my_gaussian(x)*my_gaussian(x).conjugate()
C,e = quad(func, 0, L)

def task3(x):
    
    psi_0 = (C**-(1/2))*my_gaussian(x)

    H = (dx**-2)*diags([-1, 2, -1], [-1, 0, 1], shape=(N, N))
    H = sp.lil_matrix(H)
    H[0,0]=0.5*H[0,0]
    H[N-1,N-1]=0.5*H[N-1,N-1]
    lam, phi = eigsh(H, 400, which="SM")

    a = phi.T.dot(psi_0)
    psi = phi.dot(a*np.exp(-1j*lam*tmax))
    return psi*psi.conjugate()

func = lambda x: task3(x)*x
N1,e = quad(func, 0, L)

​

​

Answer 1

最初，这似乎是一个非常简单的形状不匹配。如果将二维数组乘以一维数组，则只有当二维数组的最后一维与一维数组的大小相同时，numpy才会自动广播它们。否则，您必须显式地重塑1-d数组，使其可与2-d数组一起广播。

为了解决这个问题，我们需要知道哪些数组的形状不匹配。我将这段代码添加到您的代码中，以查看：

try:                                                                        
    psi = phi.dot(a * np.exp(-1j*lam*tmax).reshape(-1, 1))
except ValueError:
    print('phi.shape:', phi.shape)
    print('a.shape:', a.shape)
    print('lam.shape:', lam.shape)
    raise

结果是

phi.shape: (3200, 400)
a.shape: (400, 3200)
lam.shape: (400,)

因此，我们需要将lam项重塑为列向量：

np.exp(-1j*lam*tmax).reshape(-1, 1))

这解决了形状问题。但这并不能解决整个问题因为..。那么task3的输出就是一个(3200,3200)数组！当然，对于期望返回单个标量的函数的集成例程来说，这是没有用的。

最后一个问题是你必须自己解决的问题，因为我无法知道你的目标是什么。