Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Fast Numerical Integration in Python

I have a program that involves computing a definite integral many times, and have been struggling to find a way to do so quickly. The integrals I need to solve are of the following form:

\int^{b}_{a(r)}f(x)*g(x-r)dx

I have to solve this integral for many different values of r, which affects both the limits of integration and also the integrand (through the function g). Because of this, I have not found a way to vectorize the problem and must instead rely on loops. This significantly slows down the problem, because I need to make function calls in each loop. Below is one way to do it using loops (using made up data and functions):

import numpy as np 

f = lambda x: x**2
g = lambda x: np.log(x)

b=1000
r = np.arange(10,500,10)
a = 1.1*r+r**-1

def loop1(r,a):
    integration_range=[np.linspace(a[i],b,1000) for i in range(len(a))]
    out=np.zeros(len(r))
    i=0
    while i<len(r):
        out[i]=np.trapz(f(integration_range[i])*a_pdf(integration_range[i]-r[i]),integration_range[i])
        i=i+1
    return out  

This takes approximately 17.7 ms, which is too slow for my current needs. I don't care too much about getting the integrals to be super precise; I would be happy with a solution that gave approximations within 1% of the true value. Any help would be greatly appreciated!

like image 749
shadowprice Avatar asked Nov 30 '13 22:11

shadowprice


People also ask

How do you numerically integrate in Python?

Numerical integrationis provided by the quad() function of the scipy. integrate module. It takes as input arguments the function f(x) to be integrated (the “integrand”), and the lower and upper limits a and b.

Which method is best for numerical integration?

If the functions are known analytically instead of being tabulated at equally spaced intervals, the best numerical method of integration is called Gaussian quadrature. By picking the abscissas at which to evaluate the function, Gaussian quadrature produces the most accurate approximations possible.

Can Python solve integrals?

The scipy. integrate sub-package has several functions for computing integrals. The trapz takes as input arguments an array of function values f computed on a numerical grid x.

How do you double integrate in Python?

dblquad() method, we can get the double integration of a given function from limit a to b by using scipy. integrate. dblquad() method. Return : Return the double integrated value of a polynomial.


1 Answers

If you have lot's of these to do and f is more complicated than your example, you could get some benefits from memoizing f and possibly g.

What is memoization and how can I use it in Python?

Basically, anywhere you can, cache a computation and trade memory for cpu.

like image 122
Fred the Magic Wonder Dog Avatar answered Sep 30 '22 08:09

Fred the Magic Wonder Dog