Integración numérica¶
Contenidos¶
Fórmulas de cuadratura¶
Las fórmulas de integración numérica o de cuadratura son de la forma:
$$ \int_{a}^{b}f(x)dx\approx \omega_{0}\;f(x_{0})+\omega_{1}\;f(x_{1})+\cdots +\omega_{n}\;f(x_{N}) $$donde $x_{0},x_{1},...,x_{N}$ (nodos) son $N+1$ puntos distintos pertenecientes al intervalo $[a,b]$ y $\omega_{0},\omega_{1},\ldots ,\omega_{N}$ (pesos) son números reales.
Fórmulas interpolatorias¶
Si $P_{N}$ es el polinomio que interpola a $f$ en los puntos distintos $x_{0},x_{1},...,x_{N} \in \left[ a,b\right]$ y
$$ \int_{a}^{b}f(x)dx\approx \int_{a}^{b}P_N(x)dx=\omega_{0}f(x_{0})+\omega_{1}f(x_{1})+\cdots +\omega_{n}f(x_{N}) $$decimos que la fórmula de cuadratura es de tipo interpolatorio.
Fórmulas de cuadratura simples y compuestas¶
Las fórmulas de cuadratura se llaman simples si la aproximación se hace en el intervalo completo $(a,b)$, y compuestas si, antes de aplicar la fórmula, dividimos el intervalo $(a,b)$, en $n$ subintervalos.
Grado de precisión¶
Una fórmula de cuadratura tiene grado de precisión $r$ si es exacta para
$$ f\left( x\right) =1,\quad f\left( x\right) =x,\quad f\left( x\right) =x^{2},...,f\left( x\right) =x^{r}$$pero no es exacta para $f\left( x\right) =x^{r+1}$
Fórmulas de cuadratura de Newton-Cotes simples¶
Son fórmulas de cuadratura de tipo interpolatorio, eligiendo los puntos de interpolación (nodos de la fórmula) igualmente separados de una de las dos formas siguientes:
- Fórmulas cerradas Los límites de integración $a$ y $b$ son nodos de la fórmula.
- Fórmulas abiertas Ninguno de los límites de integración es nodo de la fórmula.
Calculemos primero una integral de forma exacta (simbólica) para comparar los sucesivos resultados numéricos
%matplotlib inline
Cargamos los paquetes numpy y sympy
import numpy as np
import sympy as sym
La integral
$$ \int_0^3 e^x \,dx \qquad (1) $$se calcula de forma exacta
x = sym.Symbol('x', real=True)
f = sym.exp(x)
I_exacta = sym.integrate(f,(x,0,3))
print(I_exacta)
Que también podemos escribir
I_exacta = float(I_exacta)
print(I_exacta)
Fórmula del punto medio¶
La fórmula del punto medio es
$$\int_{a}^{b}f(x)dx\approx (b-a)\;f\left( \frac{a+b}{2}\right)$$Usar la fórmula del punto medio para integrar una función en un intervalo, equivale a sustituir, dentro de la integral, la función a integrar por el polinomio de interpolación de grado cero, es decir, una recta horizontal que pasa por el punto medio de la función. Sustituimos la función por una recta e integramos. Estamos entonces calculando el área de un rectángulo.
Escribir la función punto_medio(f,a,b) que tiene como argumentos de entrada la función f a integrar y los extremos del intervalo de integración a and b, y que devuelve el valor aproximado de la integral utilizando la Regla del Punto Medio.
Calcular la integral
$$ \int_0^3 e^x \,dx$$Escribir el valor exacto (calculado arriba) y el valor aproximado.
%run Ejercicio1.py
Fórmula de los trapecios¶
La fórmula de los trapecios simple es
$$\int_{a}^{b}f(x)dx\approx \frac{b-a}{2}\;\left( f\left( a\right) +f\left( b\right) \right)$$Usar la fórmula de los trapecios para integrar una función en un intervalo, equivale a sustituir, dentro de la integral, la función a integrar por el polinomio de interpolación de grado uno, que pasa por los puntos de la función de los extremos del intervalo. Es decir, sustituimos la función por una recta e integramos. Estamos entonces calculando el área de un trapecio.
Escribir la función trapecio(f,a,b) del ejemplo anterior que tiene como argumentos de entrada la función $f$ a integrar y los extremos del intervalo de integración $a$ y $b$, y que devuelve el valor aproximado de la integral utilizando la Regla del Trapecio.
Calcular la integral
$$ \int_0^3 e^x \,dx$$Escribir el valor exacto y el valor aproximado.
%run Ejercicio2.py
Fórmula de Simpson¶
La fórmula de Simpson simple es
$$\int_{a}^{b}f(x)dx\approx \frac{b-a}{6}\;\left( f\left( a\right) +4\;f\left( \frac{a+b}{2}\right) +f\left( b\right) \right)$$Usar la fórmula de Simpson para integrar una función en un intervalo, equivale a sustituir, dentro de la integral, la función a integrar por el polinomio de interpolación de grado dos, que pasa por los puntos de la función de los extremos y el punto medio del intervalo. Es decir, sustituimos la función por una parábola e integramos.
Escribir la función simpson(f,a,b) del ejemplo anterior que tiene como argumentos de entrada la función f a integrar y los extremos del intervalo de integración a y b, y que devuelve el valor aproximado de la integral utilizando la Regla del Simpson.
Calcular la integral
$$ \int_0^3 e^x \,dx$$Escribir el valor exacto y el valor aproximado.
%run Ejercicio3.py
Fórmulas de cuadratura de Newton-Cotes compuestas¶
Una forma de disminuir el error de las fórmulas anteriores es aumentar el número de nodos utilizando las fórmulas compuestas. Estas se obtienen dividiendo el intervalo $\left[ a,b\right] $ en $n$ subintervalos y aplicando a cada uno de estos subintervalos una fórmula de cuadratura sencilla.
Si dividimos el intervalo $\left[ a,b\right] $ en $n$ subintervalos de igual longitud usando los nodos $x_{0},x_{1},...,x_{n}$,
la longitud de un subintervalo es
$$h=\frac{b-a}{n}$$los nodos son
$$x_i=a+i\,h \quad i=0,1,\ldots,n$$y el punto medio de un intervalo es
$$\bar x_i =\frac{x_{i-1}+x_i}{2}$$entonces, las fórmulas compuestas se pueden escribir:
Fórmula del punto medio compuesta¶
$$ \int_a^b f dx \approx h\sum_{i=1}^{n}f(\bar x_i)$$Escribir la función punto_medio_comp(f,a,b,n) que tiene como argumentos de entrada la función f a integrar, los extremos del intervalo de integración a y b, y el número de subintervalos que vamos a usar en la fórmula compuesta n y devuelve el valor aproximado utilizando la Regla del Trapecio compuesta.
Calcular, con n = 5 subintervalos, la integral
Escribir el valor exacto y el valor aproximado.
%run Ejercicio4.py
Regla de los trapecios compuesta¶
$$ \int_a^b f dx \approx \frac{h}{2}(f(a)+f(b))+h\sum_{i=1}^{n-1}f(x_i)$$Escribir la función trapecio_comp(f,a,b,n) que tiene como argumentos de entrada la función f a integrar, los extremos del intervalo de integración a y b y el número de subintervalos que vamos a usar en la fórmula compuesta n y devuelve el valor aproximado utilizando la Regla del Trapecio compuesta.
Calcular, con n = 4 subintervalos, la integral
Escribir el valor exacto y el valor aproximado.
%run Ejercicio5.py
Fórmula de Simpson Compuesta¶
$$ \int_a^b f dx \approx \frac{h}{6}\sum_{i=1}^{n}\left(f(x_{i-1})+4f(\bar x_i)+f(x_i)\right)) $$Escribir la función simpson_comp(f,a,b,n) que tiene como argumentos de entrada la función f a integrar, los extremos del intervalo de integración a y b y el número de subintervalos que vamos a usar en la fórmula compuesta n y devuelve el valor aproximado utilizando la Regla del Punto Medio compuesta.
Calcular, con n = 4 subintervalos, la integral
Escribir el valor exacto y el valor aproximado.
%run Ejercicio6.py
Fórmulas de cuadratura gaussianas¶
En la fórmula de cuadratura:
$$ \int_{a}^{b}f(x)dx\approx \omega_{0}f(x_{0})+\omega_{1}f(x_{1})+\cdots +\omega_{N}f(x_{N}) $$¿Es posible calcular los pesos $\omega_i$ y los nodos $x_i$ de forma que la precisión de la fórmula sea lo mayor posible? Sí, pero entonces los nodos no estarán equiespaciados.
Si $[a,b]=[-1,1]$ los pesos y nodos son
Así, por ejemplo, la fórmula gaussiana para un punto es
$$ \int_{-1}^{1}f(x)dx\approx 2\;f(0) $$Para dos puntos
$$ \int_{-1}^{1}f(x)\;dx\approx f(-0.5773502692)+f(0.5773502692) $$Para tres puntos $$ \int_{-1}^{1}f(x)\;dx\approx 0.555555556\;f(-0.774596692)+0.8888888889\;f(0)+ 0.555555556\;f(0.774596692) $$
Y así sucesivamente.
Estos resultados se pueden generalizar a cualquier intervalo $[a,b]$ cambiando los $x_i$ por $y_i$ de acuerdo con la fórmula
$$y_i=\frac{b-a}{2}x_i+\frac{a+b}{2}$$Y entonces la fórmula de cuadratura es
$$ \int_{a}^{b}f(x)dx\approx\frac{b-a}{2}\left( \omega_{0}\;f(y_{0})+\omega_{1}\;f(y_{1})+\cdots +\omega_{n}\;f(y_{n})\right) $$Escribir la función gauss(f,a,b,n) que tiene como argumentos de entrada la función f a integrar, los extremos del intervalo de integración a y b y el número de nodos n y devuelve el valor aproximado utilizando la fórmulas gaussianas con n = 1, n = 2 y n = 3 nodos la integral
Escribir el valor exacto y el valor aproximado.
%run Ejercicio7.py
Integración numérica con órdenes python¶
El módulo scipy.integrate provee varias funciones con diferentes técnicas de integración. Entre sus funciones está la función de propósito general quad.
import scipy.integrate as integrate
f = lambda x : np.exp(x)
a = 0.; b = 3;
I = integrate.quad(f,a,b)
print 'El valor aproximado es %f' % (I[0])
print 'El valor exacto es %f' % (I_exacta)