[Octave-bug-tracker] [bug #62468] quadgk.m should be able to do array-valued integration

[Michael Leitner]

Follow-up Comment #8, bug #62468 (project octave):

For the purely real case, waypoints are misleadingly named, they are just used
as dividing points for the first division of the total interval into
subintervals. The purpose is that if the integrand has integrable
singularities in the interior, it will work much better if these are at the
boundaries of the subintervals, so this is an opportunity for the user to give
their positions and have better efficiency. Thus, it is sensible to use them
in sorted order. It also would in principle work if the algorithm used them in
any random order (also if some of them are outside [a,b]), as the
corresponding contributions would just cancel (in exact arithmetics). In
finite-precision arithmetics, it is of course sensible to have them always
sorted and between a and b.

For the complex case, it is different. Here their purpose is to specify the
integration path through the complex plane. Due to the residue theorem, the
specific path you go is not relevant, all that matters is whether the poles
are on the one or the other side of the path. So if you want to do some
integral over the complex plane, you have to set the waypoints so that the
connecting line from one to the next passes on the correct side (the one you
want) of each pole. Specifically, if you do a closed contour integral around a
pole, if you flip the sense of rotation, you get a minus in the result. So no,
matlab's results are perfectly understandable, you would have seen a
difference in your non-closed contour if you had had a choice of waypoints
where you went under the pole at zero, but you didn't. Therefore you can
contract your integration path to the direction connection between a and b in
all cases by the residue theorem, and you get the same result.

For instance these two should give different results -- once you go
effectively around the pole anticlockwise but miss the way from -1+i to 1+i,
and once you effectively go once around the pole clockwise plus from -1+i to
1+i
++
quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [-1-.5i,1-.5i])
quadgk(@(x) 1./x, -1+i,1+i, 'waypoints', [1-.5i,-1-.5i])
--

Rik: waypoints are completely orthogonal to arrayvalued or not. We should add
the sorting of waypoints in the purely real case, this should make the
integration much more efficient (if they are not already sorted). And with
respect to h0: I am still reading.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62468>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/