0

我正在使用 scipy.integrate 中的 quad 从对象获取有限范围内的积分。假设目标对象处于打击中：

∫expm(A*X).expm(B*X)dx

其中 A 和 B 都是 numpy 矩阵。

为了解决这个问题，我使用了打击代码：

from scipy.integrate import quad
from scipy.linalg import expm
import numpy as np

def integrand(X, A, B):
    return np.dot(expm(A*X),expm(B*X))


A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 2], [3, 4]])

I= quad(integrand, 0, 1, args=(A,B))

但是对于结果，我得到了这个错误：

TypeError: only length-1 arrays can be converted to Python scalars

我知道当函数需要单个值但您传递一个数组时，会引发错误“只有长度为 1 的数组可以转换为 Python 标量”。但我的问题是基于数组。那么我该如何解决它。

4

2 回答 2

1

正如评论中所指出的，quad需要一个标量函数。您始终可以通过添加索引作为输出将函数传递给标量：

def integrand(X, A, B, ix=None):
    """ pass ix=None to return the matrix, ix = 0,1,2,3 to return an element"""
    output = np.dot(expm(A*X),expm(B*X))
    if ix is None:
        return output
    i, j = ix//2, ix%2
    return output[i,j]
I= np.array([quad(integrand, 0, 1, args=(A,B, i))[0]
for i in range(4)]).reshape(2,2)
I
>>array([[1031.61668602, 1502.47836021],
       [2253.71754031, 3285.33422634]])

请注意，这是非常低效的，因为您要计算 4 次积分，只要这不会打扰您。

或者，使用trapz：

x_i = np.linspace(0,1,60)
np.trapz([integrand(x, A, B) for x in x_i], x=x_i, axis=0)
>>array([[1034.46472361, 1506.62915374],
   [2259.94373062, 3294.40845422]])
于 2019-02-18T13:31:59.270 回答
0

quadpy 进行矢量化计算。expm不过，仅适用于方阵（而不适用于方阵列表）的事实需要对矩阵形状进行一些处理。

from quadpy import quad
import numpy as np
from scipy.linalg import expm


A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 2], [3, 4]])


def integrand(X):
    expAX = np.array([expm(A * x) for x in X])
    expAX = np.moveaxis(expAX, 0, -1)
    #
    expBX = np.array([expm(B * x) for x in X])
    expBX = np.moveaxis(expBX, 0, -1)
    return np.einsum("ij...,jk...->ik...", expAX, expBX)


val, err = quad(integrand, 0, 1)
print(val)

[[1031.61668602 1502.47836021]
 [2253.71754031 3285.33422633]]
于 2020-02-02T23:37:20.113 回答