SymPy多维数组计算报错：NotImplementedError问题咨询

我帮你排查了这个报错的根源，你遇到的NotImplementedError本质是SymPy的MutableDenseNDimArray元素是SymPy符号类型（哪怕值是1，也是sympy.Integer(1)），而scipy.integrate.dblquad要求被积函数必须返回纯数值。

当你用sp.lambdify生成被积函数时，H[i,j,k,l]是SymPy的整数对象，导致生成的函数返回的是SymPy表达式而非数值，SciPy无法处理这种非数值输出，因此抛出错误。下面给你几种可行的解决方案：

方案1：把SymPy数组转成数值数组（推荐） #

如果你的H数组里都是固定数值（不是符号表达式），直接转换成NumPy数组就可以解决问题——这样每个元素都是Python/numpy的数值类型，lambdify生成的函数会返回纯数值，SciPy就能正常处理了：

import scipy as sc
import numpy as np
import sympy as sp
from sympy import symbols, sin
from scipy import integrate

theta, phi = symbols('theta, phi')

# 用NumPy创建全1的四维数组，元素为浮点型
H = np.full((3,3,3,3), 1.0)
P = np.zeros_like(H)

for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                # 指定用numpy作为lambdify的后端，确保生成数值函数
                integrand = sp.lambdify([theta, phi], H[i,j,k,l]*sin(phi), modules=['numpy'])
                # dblquad返回(积分结果, 误差估计)，我们只取结果部分
                P[i,j,k,l], _ = sc.integrate.dblquad(integrand, 0, np.pi, 0, 2*np.pi)

方案2：显式转换SymPy数组元素为数值 #

如果因为某些需求必须保留SymPy数组，那在使用元素时显式转成数值类型就行：

import scipy as sc
import numpy as np
import sympy as sp
from sympy import symbols, sin, MutableDenseNDimArray
from scipy import integrate

theta, phi = symbols('theta, phi')

# 保留原SymPy数组初始化
H = MutableDenseNDimArray([1 for _ in range(81)], (3,3,3,3))
P = MutableDenseNDimArray([0 for _ in range(81)], (3,3,3,3))

for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                # 把SymPy元素转成Python浮点数值
                h_val = float(H[i,j,k,l])
                integrand = sp.lambdify([theta, phi], h_val*sin(phi), modules=['numpy'])
                result, _ = sc.integrate.dblquad(integrand, 0, np.pi, 0, 2*np.pi)
                P[i,j,k,l] = result

额外优化：针对全1数组的高效写法 #

如果你的H确实是全1的四维数组，完全没必要写四层循环浪费时间——只需要计算一次积分，再把结果赋值给整个P数组就好：

import scipy as sc
import numpy as np
import sympy as sp
from sympy import symbols, sin
from scipy import integrate

theta, phi = symbols('theta, phi')

# 计算一次积分即可
integrand = sp.lambdify([theta, phi], sin(phi), modules=['numpy'])
total_integral, _ = sc.integrate.dblquad(integrand, 0, np.pi, 0, 2*np.pi)
# 直接生成全是积分结果的数组
P = np.full((3,3,3,3), total_integral)

内容的提问来源于stack exchange，提问作者Donald Zeka