4

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.

kennytm
  • 510,854
  • 105
  • 1,084
  • 1,005
Rohin Kumar
  • 730
  • 2
  • 10
  • 24
  • I remember long ago solving this problem using np.vectorize method. I can't quite recall how I did it... But it seemed like a universal solution at the time and worked for me. Anybody can throw light on it solving in similar direction? – Rohin Kumar May 01 '17 at 12:23
  • `np.vectorize` just wraps the iteration in a function call. It doesn't speed up your code. – hpaulj May 01 '17 at 15:00
  • I used it to work on numpy arrays somehow... not for speeding up – Rohin Kumar May 01 '17 at 17:52

3 Answers3

5

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
kennytm
  • 510,854
  • 105
  • 1,084
  • 1,005
  • Interesting work-around. But having to have a '0' value thing for initial condition is kind of limiting and doesn't seem universal to me. Say, if I want to integrate with z varying from 0.1 to 0.5 for e.g. – Rohin Kumar May 01 '17 at 12:17
  • @RohinKumar You could just always [prepend](http://stackoverflow.com/questions/36998260/prepend-element-to-numpy-array) a 0. – kennytm May 01 '17 at 13:27
3

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.

Rohin Kumar
  • 730
  • 2
  • 10
  • 24
0

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.

John Zwinck
  • 239,568
  • 38
  • 324
  • 436
  • May be I didn't quite follow your solution. 0 to 1, 0 to 2 etc. is an example case. It need not always be equally spaced or be integer. So, say z=[0.1, 0.15, 0.367, 0.265...] What happens then? – Rohin Kumar May 01 '17 at 12:21
  • Well you said `z = np.arange(0,100)` but if you have some other `z` then just sort it first and my solution should still work--simply integrate from 0 to the smallest `z` value, then from that `z` to the next-smallest `z`, etc. – John Zwinck May 01 '17 at 12:28