CODEWITHSUNDEEP
pythonarraysnumpyscipynumerical-integration
Rohin Kumar
Rohin Kumar

Reputation: 810

Integrate a function with each element of numpy arrays as limit of integration

I have a function in python (using scipy and numpy also) defined as 

import numpy as np
from scipy import integrate
LCDMf = lambda x: 1.0/np.sqrt(0.3*(1+x)**3+0.7)

I would like to integrate it from 0 to each element in a numpy array say z = np.arange(0,100)

I know I can write a loop for each element iterating through like

an=integrate.quad(LCDMf,0,z[i])

But, I was wondering if there is a faster, efficient (and simpler) way to do this with each numpy element.

Upvotes: 5

Views: 4518

Answers (3)

Rohin Kumar
Rohin Kumar

Reputation: 810

After tinkering with np.vectorize I found the following solution. Simple - elegant and it works!

import numpy as np
from scipy import integrate
LCDMf = lambda x: 1.0/math.sqrt(0.3*(1+x)**3+0.7)
np.vectorize(LCDMf)

def LCDMfint(z):
    return integrate.quad(LCDMf, 0, z)

LCDMfint=np.vectorize(LCDMfint)
z=np.arange(0,100)

an=LCDMfint(z)
print an[0]

This method works with unsorted float arrays or anything we throw at it and doesn't any initial conditions as in the odeint method.

I hope this helps someone somewhere too... Thanks all for your inputs.

Upvotes: 3

kennytm
kennytm

Reputation: 523634

You could rephrase the problem as an ODE.

enter image description here

The odeint function can then be used to compute F(z) for a series of z.

>>> scipy.integrate.odeint(lambda y, t: LCDMf(t), 0, [0, 1, 2, 5, 8])
array([[ 0.        ],    # integrate until z = 0 (must exist, to provide initial value)
       [ 0.77142712],    # integrate until z = 1
       [ 1.20947123],    # integrate until z = 2
       [ 1.81550912],    # integrate until z = 5
       [ 2.0881925 ]])   # integrate until z = 8

Upvotes: 5

John Zwinck
John Zwinck

Reputation: 249552

It can definitely be done more efficiently. In the end what you've got is a series of calculations:

  1. Integrate LCDMf from 0 to 1.
  2. Integrate LCDMf from 0 to 2.
  3. Integrate LCDMf from 0 to 3.

So the very first thing you can do is change it to integrate distinct subintervals and then sum them (after all, what else is integration!):

  1. Integrate LCDMf from 0 to 1.
  2. Integrate LCDMf from 1 to 2, add result from step 1.
  3. Integrate LCDMf from 2 to 3, add result from step 2.

Next you can dive into LCDMf, which is defined this way:

1.0/np.sqrt(0.3*(1+x)**3+0.7)

You can use NumPy broadcasting to evaluate this function at many points instantly:

dx = 0.0001
x = np.arange(0, 100, dx)
y = LCDMf(x)

That should be quite fast, and gives you one million points on the curve. Now you can integrate it using scipy.integrate.trapz() or one of the related functions. Call this with the y and dx already computed, using the workflow above where you integrate each interval and then use cumsum() to get your final result. The only function you need to call in a loop then is the integrator itself.

Upvotes: 0

Related Questions