[SciPy-User] Double integration with integrate.trapz

[bhanukiran perabathini]

Hi all,

- My task is to perform a double integral (with two free variables)
and make a filled-contour plot of the logarithm of its absolute
square.

- I'm kind of sure that the code has no major mistakes. Nevertheless,
I've badly got to accelerate my double-integration algorithm.

- I bet there must be ways to increase the speed substantially. My
guess is that the over use of scipy.vectorize is the culprit.
Here I'm attaching my code and the integral, any suggestions on how
to improve the speed would be deeply appreciated.

- I've checked out integrate.dblquad and it doesn't seem to handle
complex types easily. But, no matter what, I primary intension is to
improve on the trapz method.

Thanks a tonne for every second of yours,
- Bhanukiran,
FOSSEE,
IITB, India

[Pauli Virtanen]

Hi,

On Sat, 26 Feb 2011 23:44:55 +0530, bhanukiran perabathini wrote:
[clip]

> - My task is to perform a double integral (with two free variables) and
> make a filled-contour plot of the logarithm of its absolute square.

Avoid using numpy.vectorize: it will not increase the speed. Instead, try
to vectorize your operations.

So, suppose I would need to compute the integral

I = \int_{0}^{1} dp \int_{-p}^{p**2} dq exp(-q*p*p)

I would first change the variables in the inner integral
to get rid of the varying bounds

q(x) = ((p**2) + p) * x - p
x = 0...1
dq = ((p**2) + p) * dx

I = \int_{0}^{1} dp \int_0^1 dx (dq/dx) exp(-q(x)*p*p)

Then, I would compute the integral using trapz:

---------------
from numpy import linspace, newaxis, trapz, exp

p = linspace(0, 1, 200)
x = linspace(0, 1, 200)

# instead of meshgrid, use broadcasting
pp = p[:,newaxis]
xx = x[newaxis,:]
qq = (pp**2 + pp)*xx - pp
dq_per_dx = (pp**2 + pp)

# compute the value of the integral on the square
integrand = exp(-qq*pp*pp) * dq_per_dx

# compute the result
value = trapz(trapz(integrand, pp, axis=0), x, axis=0)

# compare to scipy.integrate.dblquad
from scipy.integrate import dblquad
value2 = dblquad(lambda q, p: exp(-q*p*p), 0., 1.,
lambda p: -p, lambda p: p**2)
print value, value2[0]
# -> 0.899982875122 0.899972594941
----------------

The same trick be used for *all* double integrals.

--
Pauli Virtanen

_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Giacomo Boffi]

bhanukiran perabathini writes:

> - My task is to perform a double integral (with two free variables)
> and make a filled-contour plot of the logarithm of its absolute
> square.

notwistanding the other excellent advice you've received, please note
that your integral is symmetric with respect to the free variables, so
that you could double the resolution of your plot using the same
amount of computational resources plotting the first quadrant only

--
anch'io la penso come me, ma -- SteO153, in IHC

--
This message has been scanned for viruses and
dangerous content by MailScanner, and is
believed to be clean.

[bhanukiran perabathini]

Hi Pauli Virtanen & Giacomo Boffi,

Both of your suggestions helped me. thanks so much.

But I guess I'm having some troubles understanding and using
broadcasting to its total power. Could you suggest some tutorial where
I can read up learn this stuff?

Manipulating with the axes seems a little tricky to me. I want to
master them so well that I can think in terms of them.

Thanks again for your time,
Bhanukiran

[bhanukiran perabathini]

oh looks like this is what i've been looking for..

http://www.scipy.org/EricsBroadcastingDoc

[bhanukiran perabathini]

Hi again!

let's say we have a siple function that goes like

def func(x, y):
return x + y

now I want to get this value for every point on the x-y plane, here's what I do

x = sp.array([-2, -1, 0, 1, 2])
y = sp.array([-2, -1, 0, 1, 2])

xx = x[:, sp.newaxis]
yy = y[sp.newaxis, :]

>>> func(xx, yy)
>>>
array([[-4, -3, -2, -1, 0],
[-3, -2, -1, 0, 1],
[-2, -1, 0, 1, 2],
[-1, 0, 1, 2, 3],
[ 0, 1, 2, 3, 4]])

works as expected right?

How do we get the following function working with broadcasted arrays?

def func2(x, y):
if x > y:
return x + y
else:
return x - y

func2(xx, yy) raises the error
ValueError: The truth value of an array with more than one element is
ambiguous. Use a.any() or a.all()

How do we handle such situations?

Thanks so much,
Bhanukiran

On Mon, Feb 28, 2011 at 2:48 PM, bhanukiran perabathini

[eat]

Hi

For example

In []: x= arange(5)

In []: y= arange(5)

In []: xx= x[:, None]

In []: yy= y[None, :]

In []: A= xx+ yy

In []: A[xx< yy]= (xx- yy)[xx< yy]

In []: A

Out[]:

array([[ 0, -1, -2, -3, -4],

       [ 1,  2, -1, -2, -3],

       [ 2,  3,  4, -1, -2],

       [ 3,  4,  5,  6, -1],

       [ 4,  5,  6,  7,  8]])

Regards,

eat 

[bhanukiran perabathini]

Thanks for the reply,

I am aware of this trick though. In light of the previous discussion,
it is about 'using broadcasted arrays as inputs of an already defined
function' I'm concerned about. It was with that intention that I made
up that easy example.

Look at the following function for instance,

def func(x, y):
res = 0
for j in xrange(x):
res += y + j
return res

now what if I want to use apply this function on arrays?

by doing something as simple as func(xx, yy).
In which shape are we supposed to manufacture xx and yy for this to work?

-bhanukiran