<html>#2398: Multipole: Default (midpoint) integration method incorrectly implemented, yielding wrong results... especially with m=odd modes
<table style='border-spacing: 1ex 0pt; '>
<tr><td style='text-align:right'> Reporter:</td><td>Zach Etienne</td></tr>
<tr><td style='text-align:right'> Status:</td><td>open</td></tr>
<tr><td style='text-align:right'>Milestone:</td><td></td></tr>
<tr><td style='text-align:right'> Version:</td><td></td></tr>
<tr><td style='text-align:right'> Type:</td><td>bug</td></tr>
<tr><td style='text-align:right'> Priority:</td><td>major</td></tr>
<tr><td style='text-align:right'>Component:</td><td>EinsteinToolkit thorn</td></tr>
</table>
<p>Comment (by Roland Haas):</p>
<p>I had a look at the source code as well, and in broad strokes agree with the assessment above. There may be two different issues going on. One for the theta direction, one for the phi direction.</p>
<p>In the phi direction the functions integrated are periodic so any equidistant sampling of the [0,2pi] interval with unit weight given to the values at the sample point is a “midpoint rule” and, because of the periodicity converges exponentially well, since it is essentially the Fourier series. However the current code in the lines quoted in the description double counts the point phi=0=2pi namely one of the loops (not sure which one on top of my head) should only go to “< nx” (or “< ny”).</p>
<p>In the theta direction there is no periodicity and a weight of <code>sin(theta)</code> shows up (due to the surface element on a sphere in the physicists usual spherical coordinates). A typical midpoint rule would have ntheta segments with sample points at <code>pi/ntheta/2 + itheta * pi/ntheta</code> where the offset is to make it a midpoint. The current code is lacking the offset as well as actually using <code>ntheta+1</code> sample points a <code>itheta * pi/ntheta</code>. Sticking to an interpretation as a midpoint rule this means it includes cells centered on the poles that reach to negative theta. This would indicate a double counting this (negative theta) area in the integral. Since the weights of these segments are always zero they actually do not show up and instead the midpoint rule used ends up missing a bit of the integral (or being “open” instead of closed).</p>
<p>Note that neither one of these affects the convergence order of the method (there’s a test for this in integration_convergence_order output to <code>repos/einsteinanalysis/Multipole/test/test_22/test_midpoint_convergence_order..asc</code> and I had checked the order as well using a Maple spreadsheet before proposing #2210). </p>
<p>However it, as found here, affects the accuracy of the midpoint rule, which is (much) poorer than it ought to
<p>--<br/>
Ticket URL: <a href='https://bitbucket.org/einsteintoolkit/tickets/issues/2398/multipole-default-midpoint-integration'>https://bitbucket.org/einsteintoolkit/tickets/issues/2398/multipole-default-midpoint-integration</a></p>
</html>