2011/11/23 Eraldo Pomponi <eraldo.pomponi@gmail.com>:
Hi folks, I'm working on the integration of a function like: K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] over Voronoi cells (r is the distance from the point at which is associated the cell). I googled a lot and I found this two useful hints: http://stackoverflow.com/questions/5941113/looking-for-python-package-for-nu... http://mathforum.org/kb/message.jspa?messageID=4963570&tstart=0 but I'm still not able to understand how I should do this integration. I have the function that returns the segments ( list of [(x_star,y_start),(x_sop,y,stop)] ) necessary to construct the Voronoi cells associated to a set of points. Could someone suggest how to proceed ? There's also a numerical problem connected with the integration due to the singularity in r==0. Could you suggest which is a reasonable stable integration method available in scipy that could handle the function [1]?
Well, AISI, your function is radially symmetric, so the integration over 1/r can be done analytically within a small disk where the Gaussian is approximately one. Since: 2 \pi r * [1 / (2 \pi D r)] is just 1/D constantly. So you're left with [the integration up to] R_0/D when R_0 is the radius of that disk (where the integration is done analytically). [No more divergence, since you cut out the centre region.] You could also just integrate with the full Gaussian, since the kernel just integrates up the lengthes of the disk perimeters, if you understand. All you need to do, AISI atm, is to calculate a function that gives you the perimeter length of the circle which is inside your Voronoi cell. This is zero at r = 0 and will stay bounded for all finite radii. So integration should be fairly straightforward. You could even just do a sum over a 1D grid for the radius, since the function varies slowly, this should be fast and easy (both in runtime as in coding time). You might not even need scipy for this. Additionally, the Gaussian just introduces a suppression of radii which are farther outside. It makes the complete function bounded even for an infinitely large Voronoi cell :-) The kernel 1/r makes things easy and convenient, instead of making things troublesome: 1) because it exhibits rotational invariance; 2) because it converges nicely in 2D for an integral over a circle line I might overlook something obvious. Is the function you gave really the full function or only the "kernel", the weighting function? Friedrich P.S.: You divide your tringles resulting from the centre point and the boundary lines into two parts, separated by the closest point on the boundary, which is always in the bounds of that boundary line. Then you just have a perimeter length which is linear in r until you touch the line, and then it'll be a littel more complicated. I think you'll figure it out. P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)