Ver el código de «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» - Software Libre para la Enseñanza y el Aprendizaje de las Matemáticas (2010-11)

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 <br/>que difiere del valor obtenido en mucho menos que 0.001