[SciPy-User] Integration over Voronoi cells
Eraldo Pomponi
eraldo.pomponi@gmail....
Thu Nov 24 08:46:30 CST 2011
Dear Friedrich,
First of all let me thank you for your suggestions. It took a little bit to
get into them
but I still have some BIG doubts.
> 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:
>
Not so clear why the Gaussian is approximately one (it is true just when
r is small enough)
> 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.
Let me rewrite the function (it is just the kernel not the full function):
K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1]
This function is rotational invariant [a] so a clever way to integrate it
is to consider
a disk of radius R0. The perimeter of this disk is equal to:
l=2*pi*R0 [2]
The term 1/(2*pi*D*r) gives a constant contribution along the path eq.2,
equal to 1/D so
what is left to do is to integrate the term:
exp[-r^2/Dt]
in the interval [R0/D,R0/D+delta] and multiply by the length of the chose
path, i.e. eq.2. and sum
on a 1D grid in the interval [0,R0].
------>>>>> Is it correct? <<<<<-----------
The term exp[-r^2/Dt] can be analytically integrated (removing t for
convenience) :
g(r) = sqrt(pi)*erf(r*sqrt(1/D))/(2*sqrt(1/D)) [3]
Its defined integral in [R0/D,R0/D+delta], multiplied by the length of the
choose path,
give us:
g'(r) = 2*pi*r * [ g(r/D + delta) - g(r/D)] [4]
summing in the grid k=1..R0, we obtain the integral of eq.1 in [0,R0]:
integral([1]) = sum([g'(k) for k in np.arange(0,R0,0.001)] )
What do you think ?
>
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.
I didn't understood what you mean here. Sorry
> 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.
>
>
Yes , I hope so. When the path in eq.2 is no more circular, the integration
along it is more complicated. A bit difficult to figure out how to do that
programmatically but this is the next problem.
> P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)
>
I guess so but I think that I will fully understand this sentence just when
my
doubts will be cleared. Thanks a lot for your help.
Cheers,
Eraldo
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20111124/1fd374c3/attachment.html
More information about the SciPy-User
mailing list