[Users] Integration in Multipole

[Ian Hinder]

On 7 Oct 2013, at 21:55, Ian Hinder <ian.hinder at aei.mpg.de> wrote:

> 
> On 7 Oct 2013, at 19:28, Roland Haas <roland.haas at physics.gatech.edu> wrote:
> 
>> Ps4 mode decomposition:
>> * Jim Healy sees some odd behaviour in 2,1 mode of psi4 when using
>> Multipole thorn but not when using old GT mode decomposition thorn
>> * mode that is expected to be zero is not and leads to BH kick
>> * run uses rotatingsymmetry180 and reflectionsymmetry (in z direction)
>> * no physics change in multipole since Oersted release
>> * suggest running in full mode
>> * will work with Jim for a while then possibly take issue to mailing
>> list again
> 
> The other aspect was that this worked OK with the Oersted release, but fails with the Gauss release.  Since there were no significant changes in Multipole or AEILocalInterp between those two releases, the cause probably lies elsewhere (e.g. Carpet).
> 
> Jim, could you try running without symmetries, to see if the problem is related to symmetries?  Also, you could try setting the parameter CarpetLib::interpolate_from_buffer_zones = yes in the Gauss version?  The default changed from "yes" to "no" between these two releases. 


According to Jim, the problem was not in fact newly-introduced; the same results are obtained with earlier releases. There must have been some confusion in the telecon.  I have investigated, and the culprit appears to be in one of the integration methods used by Multipole.  Multipole interpolates data onto a sphere at regular intervals in theta and phi, and then performs a 2D numerical quadrature on these points.  It has several integration methods available, which give results of various orders of accuracy.  The default (first order accurate) integration method, labelled "midpoint" [** see below], is implemented in an unusual way.  In the following, I will discuss 1D integration, but the generalisation to 2D is straightforward.

The natural first order accurate approximation of an integral is given by summing the integrand at all the points from a to b-dx and multiplying by dx, i.e. by adding up the areas of the rectangles.  This is the analogue of the Riemann sum integral definition.  However, in Multipole, the sum runs from a to b; i.e. it "over counts" by one.  The effect of this is to introduce an additional first order error.  So the resulting error is still first order (which is what was originally tested when Multipole was developed), but it is not the most natural definition.  Additionally, this method does not have the property that a function antisymmetric about (a+b)/2 integrates numerically exactly to zero, though the integral will converge to zero as O(h)), just as the integral of any other function converges to the exact result as O(h).

In Jim's case, increasing the number of angular points in the phi direction should make the odd-m modes converge to zero if the underlying data is Pi-symmetric.  Alternatively, the other integration methods could be used, as they satisfy the symmetry property.

Typically, we set the number of points on the sphere to be so large that these issues should be negligible in comparison to other (accumulating) sources of error, but Jim's experience reminds us that this should be checked.  I have tended to worry that higher order integration methods might lead to more noise, but I haven't checked this.  It might be wise to use the higher order methods.

** Note on naming: the Wikipedia page on the "rectangle" or "midpoint" method, http://en.wikipedia.org/wiki/Rectangle_method, is very confused I think.  The expression given there is not a sum over the midpoints of the intervals, but is a sum over the left points of each interval.  The "Error" section then goes on to claim that the error is O(h^2).  This is true if you sum over the midpoints, but NOT if you sum over the left points (you can check this using a Taylor expansion).  If you sum over the left points, you get an O(h) accurate approximation of the integral.  I suspect that the first order accurate method in Multipole was named "midpoint" after a lazy wikipedia reference (probably by me).

Should we do anything about this?  The default "midpoint" method in Multipole is not only very badly named, it also does not have the desired behaviour of integrating an antisymmetric integrand exactly to zero.  

I think we should deprecate this integration method, with a warning that the default is going to change to "simpson" in future.  We could optionally implement a "rectangle" (is this the right name?) method which is the same as the current first order method but excluding the endpoints.

-- 
Ian Hinder
http://numrel.aei.mpg.de/people/hinder

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.einsteintoolkit.org/pipermail/users/attachments/20131009/366de675/attachment.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 203 bytes
Desc: Message signed with OpenPGP using GPGMail
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20131009/366de675/attachment.bin