CODEWITHSUNDEEP
pythonnumpyscipynumerical-integrationquad
samie
samie

Reputation: 191

how to fix "only length-1 arrays can be converted to Python scalars" when getting integral with an array argument

I am using quad from scipy.integrate to get an integral in a limited range from an object. suppose the target object is in the blow:

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

which both A and B are numpy matrix.

To solve this I have used blow code:

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))

But for the result I get this error:

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

I know that The error "only length-1 arrays can be converted to Python scalars" is raised when the function expects a single value but you pass an array instead. but my problem is based on array. so how can I fix it.

Upvotes: 0

Views: 462

Answers (2)

Nico Schlömer
Nico Schlömer

Reputation: 58721

quadpy does vectorized computation. The fact that expm only works on square matrices (and not on lists of square matrices) requires a bit of juggling with the matrix shapes, though.

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]]

Upvotes: 0

Tarifazo
Tarifazo

Reputation: 4343

As pointed in the comments, quad expects a scalar function. You can always pass the function to a scalar by adding the index as an output:

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]])

Note that this is very inefficient since you are calculating the integral 4 times, as long as this doesn't bother you.

Alternatively, use 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]])

Upvotes: 1

Related Questions