Obtener, mediante la regla de Simpson, el área de la función sqrt(1+x^2) entre 0 y 1 con un error menor a 0.001

En primer lugar, se definen las funciones y las constantes:

f(x):=sqrt(1+x^2);
a:0;
b:1;
E:0.001;

En virtud de la fórmula del error para la regla de Simpson:
E <= (b-a)^5*|max(f' ' ' '(x))|/(180*n^4)
Y a la vista de la gráfica de la derivada cuarta de f(x), que podemos ver gracias a Maxima y cuyo máximo, en valor absoluto, vale 3

define(g(x),diff(f(x),x,4));
plot2d(g(x),[x,-1,1],[y,-3,3]);
maxg:3;

Nos queda que para determinar el número de intervalos, sólo falta despejar la n:

solve(E=((b-a)^5*maxg)/(180*n^4),n);

Maxima ofrece las cuatro soluciones de esa ecuación de cuarto grado, de las cuales debemos tomar la entera y positiva, que nos ha salido la 4ª:

float(%[2]);

Lo que nos da un valor de n = 2.02. Por muy poco, con n=2 no tendremos la precisión deseada. Tomaremos pues n=4, ya que la Regla de Simpson sólo tiene sentido con un número par de intervalos.

n:4;

La anchura t de esos intervalos será, por tanto (b-a)/4, y los puntos quedarán recogidos en el vector puntos.

t:(b-a)/n;
puntos:[a,a+t,a+2*t,a+3*t,b];

Según la Regla de Simpson, la integral vale:
(b-a)*(f(a)+4*f(a+t)+2*f(a+2t)+2*f(a+3t)+f(b))/(3*n),
que para la función inicial y los valores de a y b:

float((b-a)*(f(a)+4*f(a+t)+2*f(a+2*t)+4*f(a+3*t)+f(b))/(3*n));

Devuelve un valor de 1.1477822. Haciendo la integral con la orden integrate, el valor ofrecido es:

float(integrate(f(x),x,0,1));

1.14779357
Que difiere del valor obtenido en mucho menos que 0.001.