Integración Numérica - Métodos del Trapezoide y Simpson

Integración Numérica - Métodos del Trapezoide y Simpson

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 p1(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 p1(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 [xj-1,xj] 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, x0=1, x1=1.5, x2=2. Ahora

Con n=4 tenemos h=(2-1)/4=0.25, x0=1, x1=1.25, x2=1.5, x3=1.75, x2=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:

n Tn(f) en=I(f)- Tn(f) en/ e2n

2 0.708333 -0.0151862 -----

4 0.697024 -0.00387663 3.91736

8 0.694122 -0.00097467 3.97738

16 0.693391 -0.000244022 3.99419

32 0.693208 -0.0000610277 3.99854

64 0.693162 -0.0000152583 3.99963

128 0.693151 -3.81467e-006 3.99991

256 0.693148 -9.53672e-007 3.99998

512 0.693147 -2.38418e-007 3.99999

1024 0.693147 -5.96046e-008 4.00000

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(h2) lo cual confirmaremos teoricamente más adelante.

Regla de Simpson: Utilizamos ahora un polinomio de interpolación cuadrático. Sea p2(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, x0=1, x1=1.5, x2=2. Ahora

Con n=4 tenemos h=(2-1)/4=0.25, x0=1, x1=1.25, x2=1.5, x3=1.75, x2=2, de modo que

...