restart;
# Figure out weights for Simpson's rule which is
# exact for quadratic polynomials
p0:=x->1;
p1:=x->x;
p2:=x->x^2;
I0:=int(p0(x),x=-1..1);
I1:=int(p1(x),x=-1..1);
I2:=int(p2(x),x=-1..1);
q:=f->w0*f(-1)+w1*f(0)+w2*f(1);
Q0:=q(p0);
Q1:=q(p1);
Q2:=q(p2);
W:=solve({I0=Q0,I1=Q1,I2=Q2},{w0,w1,w2});
# Thus the method is
subs(W,q(f));
# Find a similar rule exact for cubics
p3:=x->x^3;
q:=f->w0*f(-1)+w1*f(-1/3)+w2*f(1/3)+w3*f(1);
Q0:=q(p0);
Q1:=q(p1);
Q2:=q(p2);
Q3:=q(p3);
I3:=int(p3(x),x=-1..1);
W:=solve({I0=Q0,I1=Q1,I2=Q2,I3=Q3},{w0,w1,w2,w3});
# Thus the method is
subs(W,q(f));
# Find yet another rule that is exact for quartics
p4:=x->x^4;
I4:=int(p4(x),x=-1..1);
q:=f->w0*f(-1)+w1*f(-1/2)+w2*f(0)+w3*f(1/2)+w4*f(1);
Q0:=q(p0);
Q1:=q(p1);
Q2:=q(p2);
Q3:=q(p3);
Q4:=q(p4);
W:=solve({I0=Q0,I1=Q1,I2=Q2,I3=Q3,I4=Q4},{w0,w1,w2,w3,w4});
# Thus the method is
subs(W,q(f));
restart;
# Note that the form of the quadrature can place the location
# of the x_k's in arbitrary places as long as they are unique.
q:=f->w0*f(-sqrt(2))+w1*f(sqrt(2));
p0:=x->1;
p1:=x->x;
I0:=int(p0(x),x=-1..1);
I1:=int(p1(x),x=-1..1);
Q0:=q(p0);
Q1:=q(p1);
W:=solve({I0=Q0,I1=Q1},{w0,w1});
# Thus, an alternative to the trapezoid method is
subs(W,q(f));