ENH: improve integrate.simpson for even number of points

[andyfaff]

Closes #5618

This PR improves the accuracy of integrate.simpson when there is an even number of points (odd number of intervals).

Historical behaviour has been to use trapezoidal integration for the last/first intervals and adding it on to the _basic_simpson sum of the other intervals. This is what the even kwd was used for. When even='avg' the mean of last and first approaches was returned.

It was pointed out in #18151 that there is a more accurate way of dealing with this situation, as outlined on wikipedia. This uses a parabolic segment over the last three points to calculate the sum of the last interval, thereby using simpson behaviour for the last interval as well.

This PR:

• adds the machinery for performing the improved calculation for the last interval when there is an even number of samples.
• this machinery is accessed by adding a 'simpson' option to the even kwd.
• Since even='simpson' should offer improved behaviour in most cases (i.e. it continues the quadratic nature of the rest of the summation) this now becomes the default for the even kwd.
• Because even='simpson' offers this improved behaviour it shouldn't be necessary to access the {'first', 'last', 'avg'} options anymore. It's therefore suggested that the even kwd be deprecated and removed in a future version of scipy. After that point the 'simpson' behaviour will be used all the time.
• To achieve this deprecation the even kwdis removed from the function signature and replaced by**kwds. Existing code that uses even` will continue to work, but will receive a warning.
• If the integration axis only has two points (i.e. impossible to form a parabolic segment with the last three points) then the integration is done trapezoidally.
• Tests are added to check simpson integration over multiple axes, for odd and even number of points in the integration axis. These tests were missing and strengthen the test suite.

If there is downstream code using simpson they will:

• start getting the DeprecationWarning. This may come as a surprise as the even keyword for simpson has been present for a long time.
• get different results if they integrate with an even number of points. The new behaviour is more accurate, but people may complain that their code isn't giving the same numbers any more.

[andyfaff]

The macOS-meson workflow has started to fail, not due to changes introduced in this PR.

[j-bowhay]

Thanks, some comments around on the deprecation. I think there is a backwards compatibility issue for anyone using positional arguments (see my comment below)

[j-bowhay]

Thanks @andyfaff all the maths looks in order to me and this gives a nice improvement

[j-bowhay]

I cannot remember, was there a post on the mailing list about the deprecation aspect?

[andyfaff]

Have just written an email. I'd give it a couple of days before merging. Do you have any concerns about implementation?

[andyfaff]

@j-bowhay

[j-bowhay]

Thanks had another skim and tried the example from the original issue on this branch. This gives a clear improvement so let get it in