Comment (by Roland Haas):

I dug a bit deeper and it turns out the lower convergence order reported in repos/einsteinanalysis/Multipole/test/test_22/test_midpoint_convergence_order..asc is a bit of a red herring. The order is indeed lower for the coded up midpoint rule if given an arbitrary polynomial (which the code to compute order uses), however the physics code does not use an arbitrary polynomial but instead a function of the form f(theta) = sin(theta) g(theta) in which case (adjusting for double counting in phi) the convergence order of the coded up “midpoint” rule actually happens to be 2. This is because the code happens to implement a rule:

int(f, x) = f(x_0)*dx/2 + sum(f(x_i)*dx,i=1..N-1)) + f(x_N)*dx/2

ie a midpoint rule for the interval [x_{1/2},x_{N-1/2}] with just Riemann sums using lower / upper values for the endpoints x_0 and x_N (with incorrect weights but since the contribution is always 0 this does not matter). This also turns out to be exactly the formula for the trapezoidal rule in the theta direction.

So, looking at this now, I would say:

1. in the theta direction, given the currently fixed evaluation pointsx_i = i * pi/ntheta (used in other places of Multipole than just integration so possibly not so straightforward to change) the code does as good as it can. Depending on the point of view one takes it either implements a strange midpoint rule with Riemann sum endcaps (and incorrect weights) or the trapezoidal rule (with the same incorrect weights).
2. integration in phi is off because it double counts the phi=0=2pi point. A fix there is to reduce the integration range to “< ny” (phi is “y”). This seems to be source of convergence failure.

So, since in the phi direction, for periodic data, the trapezoidal rule and midpoint rule are identical and, as explained above, the “midpoint” rule implemented in the theta direction, for the case of a sin(theta) factor pre

--
Ticket URL: https://bitbucket.org/einsteintoolkit/tickets/issues/2398/multipole-default-midpoint-integration