quadpack.error: Supplied function does not return a valid float - when trying to double integrate over a matrix

Question

我有一堆矩阵，我必须使用它们来计算下面显示的代码中的最终 4x4 矩阵 [Ke]：

import numpy as num
from numpy import *
import scipy as sci
from scipy.integrate import *
from scipy import *
import sympy as sym
from sympy import *

x1 = 3; x2 = 0; x3 = 0; x4 = 3
y1 = 3; y2 = 3; y3 = 0; y4 = 0

K = []

s = symbols('s')
t = symbols('t')
k = symbols('k')

N1 = (s+1)*(t+1)/4
N2 = (1-s)*(t+1)/4
N3 = (1-s)*(1-t)/4
N4 = (s+1)*(1-t)/4

dN1s = diff(N1, s); dN1t = diff(N1, t)
dN2s = diff(N2, s); dN2t = diff(N2, t)
dN3s = diff(N3, s); dN3t = diff(N3, t)
dN4s = diff(N4, s); dN4t = diff(N4, t)

x = (x1*N1)+(x2*N2)+(x3*N3)+(x4*N4)
y = (y1*N1)+(y2*N2)+(y3*N3)+(y4*N4)

dxs = sym.diff(x, s)
dxt = sym.diff(x, t)
dys = sym.diff(y, s)
dyt = sym.diff(y, t)

print dxs, dxt, dys, dyt

J = sym.Matrix([[dxs, dys], [dxt, dyt]])
print J

J_inv = J.inv()
print J_inv

J_det = J.det()
print J_det

ST = sym.Matrix([[dN1s, dN2s, dN3s, dN4s], [dN1t, dN2t, dN3t, dN4t]])
print ST

XT = J_inv * ST
print XT

XT_trans = XT.transpose()
print XT_trans

k_small = sym.Matrix([[5, 0], [0, 5]])
print k_small

Ke = num.matrix(XT_trans*k_small*XT*J_det)
print Ke

for i in range(16):
    K.append([])
    K.append(dblquad(lambda t, s: Ke.item(i), -1, 1, lambda s: -1, lambda s: 1))

print K

下一步是对矩阵 [Ke] 关于 't' 和 's' 在 -1 到 1 的范围内进行双重积分。我尝试循环 [Ke] 矩阵，因为 scipy.integrate.dblquad 不直接支持矩阵。但是，当我运行代码时，出现以下错误：

quadpack.error: Supplied function does not return a valid float.

问题出在我尝试集成的那一行，但我不太明白问题出在哪里。任何帮助将不胜感激。谢谢。

PS: [Ke] 总是一个 4x4 矩阵，[k_small] 是一个对角矩阵，它的对角线可以有任何值。

Answer 1

终于找到我哪里做错了。它更多的是我用于集成的模块而不是语法。在我问题的代码中，我使用 scipy 进行集成，但没有用。原因是，scipy 中的积分只有在我指定我在 dblquad 函数中积分的方程时才有效。由于我一直在尝试遍历矩阵 Ke 并按索引调用元素，因此出现错误。现在，我修改了我的集成以使用来自 sympy 的 integrate 函数，因为它可以处理一个需要方程的变量。修改后的代码如下。

import numpy as num
from scipy import integrate
import sympy as sym

x1 = 3; x2 = 7; x3 = 6; x4 = 3
y1 = 3; y2 = 3; y3 = 4; y4 = 5

K = []

s = sym.symbols('s')
t = sym.symbols('t')
k = sym.symbols('k')

N1 = (s+1)*(t+1)/4
N2 = (1-s)*(t+1)/4
N3 = (1-s)*(1-t)/4
N4 = (s+1)*(1-t)/4

dN1s = sym.diff(N1, s); dN1t = sym.diff(N1, t)
dN2s = sym.diff(N2, s); dN2t = sym.diff(N2, t)
dN3s = sym.diff(N3, s); dN3t = sym.diff(N3, t)
dN4s = sym.diff(N4, s); dN4t = sym.diff(N4, t)

x = (x1*N1)+(x2*N2)+(x3*N3)+(x4*N4)
y = (y1*N1)+(y2*N2)+(y3*N3)+(y4*N4)

dxs = sym.diff(x, s)
dxt = sym.diff(x, t)
dys = sym.diff(y, s)
dyt = sym.diff(y, t)

print dxs, dxt, dys, dyt

J = sym.Matrix([[dxs, dys], [dxt, dyt]])
print J

J_inv = J.inv()
print J_inv

J_det = J.det()
print J_det

ST = sym.Matrix([[dN1s, dN2s, dN3s, dN4s], [dN1t, dN2t, dN3t, dN4t]])
print ST

XT = J_inv * ST
print XT

XT_trans = XT.transpose()
print XT_trans

k_small = sym.Matrix([[k, 0], [0, k]])
print k_small

Ke = num.matrix(XT_trans*k_small*XT*J_det)
print Ke

for i in range(16):
    K.append(integrate(Ke.item(i), (t, -1, 1), (s, -1, 1)))

print K

请注意，我实际上在矩阵 k_small 中使用了变量 k 并且它仍然有效（这是我的实际问题，因为集成应该能够处理变量代替 k 和 return 一个本身就是变量的答案）。

quadpack.error: Supplied function does not return a valid float - when trying to double integrate over a matrix

quadpack.error: Supplied function does not return a valid float - when trying to double integrate over a matrix

python

numpy

matrix

sympy

scipy