Problems with computing Integral in Python

Question

So I've been trying to use the general integration method (quad) from Scipy on numpy on a specific expression; however, I get the following error

TypeError: unsupported operand type(s) for *: 'function' and 'float'

Here is the function I'm trying to integrate (seems like theres no mathjax here):

integral from 0 to a of t*y(t) * J0(u_{i} * t/a) dt where a is the length of y(t), J0 is a zeroth-order Bessel Function and u_{i} is the roots of J0(u) = 0

My code:

import numpy as np
import scipy.integrate as integrate
from scipy.special import *
y = np.load('dataset.npy')  # 1000 pts long
a = 1000
roots = jn_zeros(0, 1000)
func = lambda t: t * jv(0, t * (roots/a) )
result = integrate.quad(func * y , 0, a)

kabanus · Accepted Answer · 2018-08-13 18:10:07Z

0

The first problem is you are not calling your function, what you meant was

func(y)*y

otherwise you are just multiplying a function name by a number/array, which is the error you are seeing. Beyond that, the initial argument to integrate.quad is a callable, which is the function you want to integrate. It does not accept a vector, so your initial y vector is completely ignored. If anything, you want this:

result = integrate.quad(lambda y: func(y) * y , 0, a)

but this ignored the dataset completely. To integrate on samples use trapz, or simps, which accept an x array, and a y array. You probably need to generate your y array using np.linspace or something similar.

edited Aug 13, 2018 at 18:10

answered Aug 13, 2018 at 6:19

kabanus

26.3k7 gold badges48 silver badges79 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

kabanus Over a year ago

I noticed the arr typo, but let's make it a pun intended.

Sam Over a year ago

Thanks! but now it seems like I got the following error: ValueError: operands could not be broadcast together with shapes (5,) (1000,)

kabanus Over a year ago

@Sam Looks like func is returning a matrix. You did not post enough info to debug this (minimal reproducible example). I suggest you write a new question for your new problem, but only after you looked around at what func is returning and if it makes sense to you.

Collectives™ on Stack Overflow

Problems with computing Integral in Python

1 Answer 1

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related