En esta lección comenzamos el estudio de métodos numéricos para el cálculo numérico de integrales de la forma

Un método común para aproximar I(f) es reemplazando f(x) con un polinomio de interpolación. Este procedimiento se conoce como las reglas de Cuadratura de Newton. Examinamos los primeros dos casos de este método donde se usan polinomios de interpolación lineales y cuadráticos.

Método del trapezoide: Sea p₁(x) el polinomio lineal que interpola a f(x) en x=a y x=b, i.e.,

Usando la fórmula para el area de un trapezoide o integrando p₁(x) directamente se obtiene que

Asi que podemos escribir la aproximación:

(*)

Más adelante análizamos en detalles el error en esta aproximación. Por el momento basta observar que la aproximación es buena siempre que f sea aproximadamente lineal. En el caso general, dividimos el intervalo [a,b] en subintervalos más pequeños y aplicamos la fórmula anterior en cada subintervalo. Si los subintervalos son suficientemente pequeños, entonces f es aproximadamente lineal en cada subintervalo y la aproximación es buena. Definimos el largo de los subintervalos por:

El j-esimo subintervalo esta dado por [x_j-1,x_j] donde

Podemos escribir ahora que:

Usando la aproximación (*) podemos escribir

Usando esto en la fórmula anterior, obtenemos que

Esto se conoce como la regla (compuesta) del trapezoide para aproximar I(f).

Ejemplo 1: Usando la regla del trapezoide con n=2 y n=4 aproximamos:

cuyo valor exacto es correcto al número de cifras mostradas. Para n=2 tenemos que h=(2-1)/2=0.5, x₀=1, x₁=1.5, x₂=2. Ahora

Con n=4 tenemos h=(2-1)/4=0.25, x₀=1, x₁=1.25, x₂=1.5, x₃=1.75, x₂=2, de modo que

Estos calculos los podemos realizar también utilizando la función trapz de MATLAB. En el siguiente programa no solo calculamos los dos resultados de arriba sino que generamos una tabla de errores (exactos) para varios valores de n aprovechando que en este ejemplo tenemos el valor exacto del integral:

iexacto=log(2);
n=2;
error1=0;
for i=1:10
x=linspace(1,2,n+1);
y=1./x;
iaprox=trapz(x,y);
error=iexacto-iaprox;
ratio=error1/error;
disp(['n=' num2str(n) ', iaprox=' num2str(iaprox,6) ',error=' num2str(error,6) ',ratio=' num2str(ratio,6)])
n=2*n;
error1=error;
end

Los resultados fueron como sigue:

Estos resultados confirman claramente la convergencia del método del trapezoide en este ejemplo particular. Podemos ver que cada ves que se duplica la n, lo cual equivale a dividir la h entre dos, el error disminuye por un factor de cuatro aproximadamente (última columna de la tabla) esto es característico de convergencia O(h²) lo cual confirmaremos teoricamente más adelante.

Regla de Simpson: Utilizamos ahora un polinomio de interpolación cuadrático. Sea p₂(x) el polinomio de grado (a lo más) dos que interpola a f(x) en x=a, x=(a+b)/2, x=b. Este polinomio se puede escribir como:

Tenemos ahora que

Pero con h=(b-a)/2 y u=x-a tenemos que

En forma similar se obtiene que

Tenemos pues que

(**)

Argumentando en forma similar a en método del trapezoide, tenemos que si n es un entero par (¿por qué?) entonces

Usando la fórmula (**) podemos aproximar

Ahora

Esta fórmula se conoce como la regla (compuesta) de Simpson para aproximar a I(f).

Ejemplo 2: Usando la regla de Simpson con n=2 y n=4 aproximamos:

cuyo valor exacto es correcto al número de cifras mostradas. Para n=2 tenemos que h=(2-1)/2=0.5, x₀=1, x₁=1.5, x₂=2. Ahora

Con n=4 tenemos h=(2-1)/4=0.25, x₀=1, x₁=1.25, x₂=1.5, x₃=1.75, x₂=2, de modo que

MATLAB no tiene una rutina simp equivalente a trapz. ¡Tiene una mejor llamada quad! La subrutina quad utiliza una regla de Simpson adaptativa donde el valor de h se ajusta para que el error en la aproximación satisfaga una tolerancia especificada por el usuario. También MATLAB tiene la subrutina quad8 que al igual que quad usa un método adaptativo pero con una fórmula de aproximación de grado mayor. En lugar de usar estas rutinas que hacen las comparaciones un tanto complicadas, implementamos nuestra versión de simp equivalente a trapz:

function q=simp(x,y);
n=length(x)-1;
if (n/2)~=floor(n/2)
disp('n tiene que ser par');
break;
end
h=x(2)-x(1);
v=2*ones(n+1,1);
v2=2*ones(n/2,1);
v(2:2:n)=v(2:2:n)+v2;
v(1)=1;v(n+1)=1;
q=(h/3)*y*v;

Esta subrutina implementa una forma vectorizada del método de Simpson que ejecuta eficientemente en MATLAB. Note que se requiere que n sea par. Recuerde también que en MATLAB los indices de los arreglos corren empezando en uno. El mismo programa del Ejemplo 1 lo podemos usar aqui ahora reemplazando la llamada a trapz por simp. Obtuvimos los siguientes resultados:

Estos resultados confirman claramente la convergencia de la regla de Simpson en este ejemplo particular. Podemos ver que cada ves que se duplica la n, lo cual equivale a dividir la h entre dos, el error disminuye por un factor de 16 aproximadamente (última columna de la tabla) esto es caracteristico de convergencia O(h⁴) lo cual confirmaremos teoricamente más adelante.

Ejercicios:

1. Usando las reglas del trapezoide y de Simpson y los programas descritos en esta lección, aproxime el siguiente integral:

   El valor exacto de este integral es p/4. Use esto para generar una tabla con las aproximaciones y los errores (exactos) y estime el orden de convergencia.

2. La regla del punto medio se puede usar para obtener la siguiente aproximación de I(f):

   Usando esta fórmula diseñe una fórmula compuesta llamada la regla (compuesta) del punto medio. Use esta fórmula en el ejercicio anterior y estime el orden de convergencia de la misma.

3. Utilizando un polinomio cúbico de Hermite para interpolar a f(x) en [a,b], desarrolle una fórmula para aproximar a I(f). Obtenga la fórmula compuesta correspondiente.