Saving nodes in iterative Adaptive Simpson quadrature method (Matlab)

Question

Implementing on Matlab the Simpson adaptive rule to approximate an integral (the following code), I am struggling with saving nodes correctly. I have tried different solutions, but none of them seems to be the right one. Any suggestions?

[S, nodes] = SimpsonAdapt (f, a, b, tol, hmin)

h = b - a;
c = (a + b) / 2;
d = (a + c) / 2;
e = (c + b) / 2;
fa = feval (f, a);
fb = feval (f, b);
fc = feval (f, c);
fd = feval (f, d);
fe = feval (f, e);
nodes = [a b];


S1 = (h/6)*(fa + 4*fc + fb);
S2 = (h/(12))*(fa + 4*fd + 2*fc + 4*fe + fb);

E2 = (1/(15))*(S2 - S1);

if abs(E2) <= tol && (b - a) < hmin
    S = S2 + E2;
    nodes =[];
else
    [Q1, nodi1] = SimpsonAdapt (f, a, c, tol/2, hmin, nodes);
    [Q2, nodi2] = SimpsonAdapt (f, c, b, tol/2, hmin, nodes);
    nodes = [nodes, nodi1];
    nodes = [nodes, nodi2];
    
    S = Q1 + Q2;
end
nodes = unique(nodes);
end 

Answer 1

The nodes list that is returned is the list of points used to evaluate the quadrature formula. You should decide which nodes these are and use that decision uniformly. You first give only the interval end-points. However, in the branching you return the 6th order Richardson extrapolation of the Simpson values, which uses all 5 points [a,d,c,e,b].

In the composition, you simply have to remove the repeated middle point c once,

nodes = [ node1(1:end-1) node2 ]

See also https://stackoverflow.com/questions/58892251/building-a-recursive-function-in-matlab-to-draw-adaptive-trapezoidal-quadrature for a similar question.

Using the adapted python code (which does not care about saving function evaluations)

def adapt_simpson_quad(f,a,b,tol):
    def recurse_simp(a,b,tol):
        x = np.linspace(a,b,5)
        y = f(x)
        s1 = (b-a)*(y[0]+4*y[2]+y[4])/6
        s2 = (b-a)*(y[0]+4*y[1]+2*y[2]+4*y[3]+y[4])/12
        if abs(s2 - s1) < tol: return list(x[::2]), s1
        if abs(s2 - s1) < 15*tol: return list(x), s2
        node1, approx1 = recurse_simp(a,x[2],tol/2);
        node2, approx2 = recurse_simp(x[2],b,tol/2);
        return node1[:-1] + node2, approx1+approx2;
    return recurse_simp(a,b,tol);

gives reduced node lists without any clean-up and results in plots like

for the integration task

f = lambda x: (cos(x**0.25))/(1+x**1.5)
a=1; b=1.71;

or more demonstrative

for the task

f = lambda x: sin(pi*x**2);
a=0; b=1;