[Users] Integration in Multipole

[Ian Hinder]

On 9 Oct 2013, at 20:14, Erik Schnetter <schnetter at cct.lsu.edu> wrote:

> On 2013-10-09, at 13:16 , Ian Hinder <ian.hinder at aei.mpg.de> wrote:
> 
>> 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.
> 
> The domain in the phi direction is periodic. Introducing an arbitrary boundary in the phi direction does not make sense. The only sensible integration method is to take the sum over all points in the phi direction (and to divide by N).

In our case, we interpolate from 0 to 2 Pi, so we need to exclude the last (repeated) point.  Come to think of it, this is fairly dumb for the phi direction, but it makes the code simpler to deal with theta and phi uniformly.  Maybe this should be changed (at the cost of making all waveforms computed previously unreproducible).

> If you do not want to take this as granted, then consider the following:
> - choose are arbitrary (!) integration method in the phi direction
> - perform this integration N times, moving the integration boundary by 1 point for each step. this is possible since the domain is symmetric. each point then becomes the boundary for one of these steps.
> - take the average of these N integrations
> - by linearity, these N integrations can be collapsed to a single integration
> - by symmetry, the coefficient for each grid point for this collapsed integration has to be the same
> - qed

This just shows that there exists an integration method which has all the weights the same, which is already known.  I don't think it says anything about the individual methods.

> Put differently: In a periodic grid, no point is special. Hence, the integration weight for each point must be the same.
> 
> This corresponds to a spectral method (using a Fourier basis), with exponential convergence.
> 
> We should remove all other phi integration methods.


It looks like you know much more about this than I do.  I didn't know about this result for quadrature of periodic functions, but now that I think about it some more, it makes sense.  The integration methods in Multipole are the classical quadrature methods for \int_a^b f(x) dx with no additional information about f.

My worry about using a highly accurate (i.e. high order) method with a small number of points is that noise at some of these small number of points could significantly affect the result.  Do you know anything about this?  Am I worrying for nothing?  

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

-------------- 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/0fe1ebfb/attachment-0001.bin