How to make my python integration faster?

Question

Hi i want to integrate a function from 0 to several different upper limits (around 1000). I have written a piece of code to do this using a for loop and appending each value to an empty array. However i realise i could make the code faster by doing smaller integrals and then adding the previous integral result to the one just calculated. So i would be doing the same number of integrals, but over a smaller interval, then just adding the previous integral to get the integral from 0 to that upper limit. Heres my code at the moment:

import numpy as np                              #importing all relevant modules and functions
from scipy.integrate import quad
import pylab as plt
import datetime
t0=datetime.datetime.now()                      #initial time
num=np.linspace(0,10,num=1000)                  #setting up array of values for t
Lt=np.array([])                                 #empty array that values for L(t) are appended to
def L(t):                                       #defining function for L
    return np.cos(2*np.pi*t)
for g in num:                                   #setting up for loop to do integrals for L at the different values for t
    Lval,x=quad(L,0,g)                          #using the quad function to get the values for L. quad takes the function, where to start the integral from, where to end the integration
    Lv=np.append(Lv,[Lval])                     #appending the different values for L at different values for t

What changes do I need to make to do the optimisation technique I've suggested?

Should Lt be called Lv? Otherwise Lv is not initialised before append method is called. — Moustache, Aug 16 '17 at 11:40
You might also think about using ctypes, which I think should speed up integration as well: https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#faster-integration-using-low-level-callback-functions — wordsforthewise, Jul 17 '18 at 02:18

Sam Mussmann · Accepted Answer · 2013-12-10T15:25:05.663

5

Basically, we need to keep track of the previous values of Lval and g. 0 is a good initial value for both, since we want to start by adding 0 to the first integral, and 0 is the start of the interval. You can replace your for loop with this:

last, lastG = 0, 0
for g in num:
  Lval,x = quad(L, lastG, g)
  last, lastG = last + Lval, g
  Lv=np.append(Lv,[last])

In my testing, this was noticeably faster.

As @askewchan points out in the comments, this is even faster:

Lv = []
last, lastG = 0, 0
for g in num:
  Lval,x = quad(L, lastG, g)
  last, lastG = last + Lval, g
  Lv.append(last)
Lv = np.array(Lv)

edited Dec 10 '13 at 15:25

answered Dec 10 '13 at 02:13

Sam Mussmann

5,883
2
29
43

1

Nice trick, more than *noticeably* faster :P By the way, you can increase the speed by another 25% or so by starting with `Lv = []`, and at the end of the loop `Lv.append(last)`, and after the loop `Lv = np.array(Lv)`. – askewchan Dec 10 '13 at 02:47
`Lv1` is not initialized. – Warren Weckesser Dec 10 '13 at 02:51
@WarrenWeckesser: Oops, `Lv1` should be `Lv`. I've edited the answer. – Sam Mussmann Dec 10 '13 at 15:20
This works well but Scipy has an in-built function (see below) which generally outperforms any native-Python implementation. – Moustache May 29 '18 at 19:15

Moustache · Answer 2 · 2017-08-16T11:37:31.343

Using this function:

scipy.integrate.cumtrapz

I was able to reduce time to below machine precision (very small).

The function does exactly what you are asking for in a highly efficient manner. See docs for more info: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.cumtrapz.html

The following code, which reproduces your version first and then mine:

# Module Declarations
import numpy as np                             
from scipy.integrate import quad
from scipy.integrate import cumtrapz
import time

# Initialise Time Array                      
num=np.linspace(0,10,num=1000)

# Your Method    
t0 = time.time()                  
Lv=np.array([])                                
def L(t):                                      
    return np.cos(2*np.pi*t)
for g in num:                                  
    Lval,x=quad(L,0,g)                        
    Lv=np.append(Lv,[Lval])   
t1 = time.time()
print(t1-t0)

# My Method
t2 = time.time()
functionValues = L(num)
Lv_Version2 = cumtrapz(functionValues, num, initial=0)
t3 = time.time()
print(t3-t2)

Which consistently yields:

t1-t0 = O(0.1) seconds

t3-t2 = 0 seconds

How to make my python integration faster?

2 Answers2