Gregzor - 1 year ago 128
Python Question

# Double integral with function and sampled data Python

I am looking for a way to perform double integration on sampled data using numpy trapz or a similar function from the scipy stack.

In particular, I would like to calculate function:

where

`f(x',y')`
is the sampled array and
`F(x, y)`
is an array of the same dimension.

This is my attempt, which gives incorrect results:

``````def integrate_2D(f, x, y):
def integral(f, x, y, x0, y0):
F_i = np.trapz(np.trapz(np.arcsinh(1/np.sqrt((x-x0+0.01)**2+(y-y0+0.01)**2)) * f, x), y)
return F_i
sigma = 1.0
F = [[integral(f, x, y, x0, y0) for x0 in x] for y0 in y]
return F

xlin = np.linspace(0, 10, 100)
ylin = np.linspace(0, 10, 100)
X,Y = np.meshgrid(xlin, ylin)

f = 1.0 * np.exp(-((X - 8.)**2 + (Y - 8)**2)/(8.0))
f += 0.5 * np.exp(-((X - 1)**2 + (Y - 9)**2)/(10.0))

F = integrate_2D(f, xlin, ylin)
``````

The output array seems to be oriented towards the diagonal of the resulting grid, while it should rather return an array that looks like blurred input array.

I can see what you're trying to do, but the nesting is hiding the logic. Try something like this,

``````def int_2D( x, y, xlo=0.0, xhi=10.0, ylo=0.0, yhi=10.0, Nx=100, Ny=100 ):

# construct f(x,y) for given limits
#-----------------------------------
xlin = np.linspace(xlo, xhi, Nx)
ylin = np.linspace(ylo, yhi, Ny)
X,Y = np.meshgrid(xlin, ylin)

f = 1.0 * np.exp(-((X - 8.)**2 + (Y - 8)**2)/(8.0))
f += 0.5 * np.exp(-((X - 1)**2 + (Y - 9)**2)/(10.0))

# construct 2-D integrand
#-----------------------------------
m = np.sqrt( (x - X)**2 + (y - Y)**2 )
y = 1.0 / np.arcsinh( m ) * f

# do a 1-D integral over every row
#-----------------------------------
I = np.zeros( Ny )
for i in range(Ny):
I[i] = np.trapz( y[i,:], xlin )

# then an integral over the result
#-----------------------------------
F = np.trapz( I, ylin )

return F
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download