<p dir="auto"><strong>### Problem</strong><br>
As I used simpson in scipy, I found there some issue about it. Simpson's 1/3 rule use parabola eq. to approximate every curve in the segment, to determine the parabola, it needs 3 points, so that is the constraint of 1/3 rule. It just apply to odd number well, but inaccurate as the number is even.<br>
<strong>### Solution</strong><br>
Combine 1/3 and 3/8 rule, so that simpson's rule apply every number of point.<br>
<strong>### Result</strong></p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy as np
from scipy.integrate import simps
def f(x):
    return x**3


x = np.linspace(1.0, 4.0, 20)
y = f(x)
a = simps(y,x)

x = np.linspace(1.0, 4.0, 19)
y = f(x)
b = simps(y,x)

print (a, b)

correct result is 63.75"><pre class="notranslate"><code class="notranslate">import numpy as np
from scipy.integrate import simps
def f(x):
    return x**3


x = np.linspace(1.0, 4.0, 20)
y = f(x)
a = simps(y,x)

x = np.linspace(1.0, 4.0, 19)
y = f(x)
b = simps(y,x)

print (a, b)

correct result is 63.75
</code></pre></div>
<p dir="auto">result1 should be more accurate than result2 because of the smaller h,<br>
but result1 = 63.7549, result2 = 63.7499<br>
I just write a C version and I uploaded it in my github, there is also the derivation of simpson rule.<br>
<a href="https://github.com/JrPhy/numerical/tree/master/Integration">https://github.com/JrPhy/numerical/tree/master/Integration</a></p>
<h1 dir="auto">problem</h1>
<p dir="auto">The Simpson's rule is suitable for arrays with an odd number of points. The current implementation allows for even number of points by integration the first or last interval using the trapezium rule, which is highly inaccurate.</p>
<h1 dir="auto">solution</h1>
<p dir="auto">In case of an even number n points the integral over the first n-1 points can be calculated using the normal Simpson scheme, i.e. the weights are 1/12_(2, 8, 2) over all intervals of three points. To calculated the contribution of the last interval, the 2nd order polynomial fit over the last 3 points is integrated over the interval between the last 2 points. This results in a weight scheme of 1/24_[-1, 8, 5] for the last inverval.</p>
<h1 dir="auto">result</h1>
<p dir="auto">Accuracy with n+1 points is better than with n points, even if n+1 is even. This is demonstrated for int_0^pi sin(x) dx in the attached code.</p>
<p dir="auto"><a href="https://github.com/scipy/scipy/files/70330/BetterSimps.txt">BetterSimps.txt</a></p>