Why is my adaptive quadrature worse than the composite Simpson rule?

Question

I wrote a basic adaptive quadrature that uses the difference between trapezoidal and Simpson rule to estimate the error on each intervall and makes the mesh finer if needed:

def adaptQuad(f, Meshes, tol): 
 M = Meshes[-1]
 h = np.diff(M) #Calculate intervall lengths h_i
 mp = (M[1:]+M[:-1])/2 #calculate midpoints
 fx = f(M); fm = f(mp) #evaluate f(x)
 trp = h*(fx[:-1]+2*fm+fx[1:])/4 #evaluate trapezoidal (on each interval)
 simp = h*(fx[:-1]+4*fm+fx[1:])/6 #evaluate Simpson (on each interval)
 err = abs(trp-simp) #Difference between the rules
 if sum(err) > tol: #are we done? 
    ref = np.nonzero(err > sum(err)/len(err))[0] #Refine mesh where needed
    Meshes.append(np.sort(np.append(M,mp[ref])))
    return adaptQuad(f, Meshes, tol) #recursion
 return sum(simp), Meshes #return composite simpson

I wanted to compare it to the composite simpson rule:

def compSimpson(f,a,b,N):
 x,h = np.linspace(a,b,N+1,retstep=True)
 midP = (x[1:]+x[:-1])/2
 x = x[1:-1]
 return h/6 *(f(a) + 2*sum(f(x))+4*sum(f(midP))+f(b))

But when I tried it on:

f = lambda x: 1/(1+25*x**2)

with the following:

Q,M = adaptQuad(f,[np.linspace(-1,1,5)],0.001)
Q_simp = compSimpson(f,-1,1,len(M[-1])-1) #Try Comp. Simpson with as many intervals
I = integrate.quad(f,-1,1)[0]
print("Error composite Simpson: ", abs(Q_simp-I))         
print("Error adaptive quad: ", abs(Q-I))

I get: Error composite Simpson: 5.274740111360643e-08 Error adaptive quad: 1.2452016172281866e-05

I tried some different total tolerances and most of the time the composite simpson rule outperforms the adaptive quadrature? This shouldn't be possible but I can't find my error in the implementation either, I'd be glad for some help:)

Edit: Here an illustration of how the Mesh evolves over time (the nodes are the dots, each line of same coloured dots shows one Meshes[i])

Answer 1

As commented, your refinement strategy is skewed. It prefers small intervals over large, or in other words, it is too eager to split large intervals.

What you really should want is that each interval contributes to the error according to their length, that is, a constant error density. For that purpose I would change your code to

def adaptQuad(f, Meshes, tol): 
  M = Meshes[-1]
  h = np.diff(M) #Calculate interval lengths h_i
  mp = (M[1:]+M[:-1])/2 #calculate midpoints
  fx = f(M); fm = f(mp) #evaluate f(x)
  trp = h*(fx[:-1]+2*fm+fx[1:])/4 #evaluate trapezoidal (on each interval)
  simp = h*(fx[:-1]+4*fm+fx[1:])/6 #evaluate Simpson (on each interval)
  err = abs(trp-simp) #Difference between the rules

  ref = err > tol*h
  if not ref.any():  #are we done? 
    return sum(simp), Meshes #return composite Simpson

  Meshes.append(np.sort(np.append(M,mp[ref])))
  return adaptQuad(f, Meshes, tol) #recursion

This equalizes the results somewhat,

Error composite Simpson:  1.4457746544405126e-07
Error adaptive quad:  7.845544380824521e-07

Both are far below the desired error bound 1e-3. But already for tol=1e-4 there is again the marked difference towards the composite method,

Error composite Simpson:  1.935168691957756e-10
Error adaptive quad:  3.361578637672835e-08

The second point to consider is a little more subtle, the adaptive strategy produces a subdivision for the error of the trapezoidal method, but returns the result for the Simpson method. While the result is better, the error bound tol*interval_length is not valid or very pessimistic.

Another effect is that a greater density of small segments is required for the trapezoidal method in "stiff" regions than for the Simpson method. This means that the adaptive strategy has unnecessary computations in these stiff regions, while the equidistant composite Simpson still has sufficient segment density there, and greater density in more flat regions, thus giving an overall better result.

Or in other words, if you want a fair comparison, take the composite trapezoidal method with the same number of function evaluations.

See Saving nodes in iterative Adaptive Simpson quadrature method (Matlab) for another approach where the recursion actually splits the interval and the requested error bound.