INTEGRACION

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 área de un trapezoide o integrando p1(x) directamente se obtiene que

Asi que podemos escribir la aproximación:

(*)

Más adelante analizamos 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 cálculos 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 teóricamente 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

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:

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

2 0.694444 -0.00129726 -----

4 0.693254 -0.000106788 12.1481

8 0.693155 -7.35009e-006 14.5288

16 0.693148 -7.35009e-006 14.5288

32 0.693147 -2.97299e-008 15.885

64 0.693147 -1.86151e-009 15.9708

128 0.693147 -1.16398e-010 15.9927

256 0.693147 -7.27562e-012 15.9983

512 0.693147 -4.54747e-013 15.9993

1024 0.693147 -2.84217e-014 16.0000

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 característico de convergencia O(h4) lo cual confirmaremos teóricamente

...