Autores: Iván Sánchez, Samuel Ochoa.
La integración numérica constituye una amplia gama de algoritmos para calcular el valor numérico de una integral definida. Para integrales unidimensionales, existen las llamadas fórmulas de Newton-Cotes, que son métodos de interpolación con polinomios evaluados en puntos igualmente separados.
Algunos métodos de la familia de fórmulas de Newton-Cotes son la regla del rectángulo, la regla del trapecio, la regla del trapecio compuesta, la regla de Simpson y la regla de Simpson 1/3 compuesta. A continuación se presenta la implementación de cada uno de estos algoritmos en Python (excepto la regla del rectángulo, la más básica):
def trapecio_simple(f, A, B):
return (B-A) * ((f(A) + f(B))/2)
def trapecio_compuesto(f, A, B, n):
k = 1
valor = 0
compuesta = 0
while k <= n - 1:
c = A + k*(B-A)/n
valor += f(c)
k += 1
compuesta = ((B-A)/n)*((f(A) + f(B))/2.0 + valor)
return compuesta
def simpson(f, A, B):
h = (B-A)/6
return h * (f(A) + 4*f((A+B)/2) + f(B))
def simpson13(f, A, B, n):
h = (B-A) / n
suma = 0.0
for i in range(1, n):
#x = a - h + (2 * h * i)
x = A + i * h
if i % 2 == 0:
suma += 2 * f(x)
else:
suma += 4 * f(x)
suma += f(A) + f(B)
return suma * (h / 3)
La asignación consiste de una serie de ejercicios de integración utilizando métodos numéricos, los cuales se presentan a continuación:
Resuelva la siguiente integral de forma analítica y luego compare el resultado obtenido por la regla Simpson para n = 6
Para simplificar la resolución analítica, nos ahorramos los cambios en los límites de integración resolviendo primero sin ellos:
$$\int_{0}^{1}\!\left(1 + e^{-x}\sin(4x)\right)\,dx = F(1) - F(0)$$$$ \begin{align*} F(x) &= \int\!\left(1 + e^{-x}\sin(4x)\right)\,dx \\ &= \int\!dx + \int\!e^{-x}\sin(4x)\,dx \\ &= x + \int\!e^{-x}\sin(4x)\,dx \end{align*} $$Para resolver la integral de la derecha se utiliza integración por partes, derivando la función seno e integrando la función exponencial: $$ \begin{align*} u = \sin(4x) &\Rightarrow du = 4 \cos(4x)\,dx \\ dv = e^{-x}\,dx &\Rightarrow v = -e^{-x} \end{align*} $$
$$ \begin{align*} \int\!e^{-x}\sin(4x)\,dx &= uv - \int\!v\,du \\ &= -e^{-x}\sin(4x) - \int\!\left(-4 e^{-x} \cos(4x)\right)\,dx \\ &= -e^{-x}\sin(4x) + 4 \int\!e^{-x} \cos(4x)\,dx \end{align*} $$Esta integral también requiere integración por partes, derivando la función coseno e integrando la función exponencial: $$ \begin{align*} u = \cos(4x) &\Rightarrow du = -4 \sin(4x)\,dx \\ dv = e^{-x}\,dx &\Rightarrow v = -e^{-x} \end{align*} $$
$$ \begin{align*} \int\!e^{-x}\sin(4x)\,dx &= -e^{-x}\sin(4x) + 4 \int\!e^{-x} \cos(4x)\,dx \\ &= -e^{-x}\sin(4x) + 4 \left[uv - \int\!v\,du\right] \\ &= -e^{-x}\sin(4x) + 4 \left[-e^{-x}\cos(4x) - \int\!4e^{-x}\sin(4x)\,du\right] \\ &= -e^{-x}\sin(4x) - 4 e^{-x}\cos(4x) - 16\int\!e^{-x}\sin(4x)\,du \end{align*} $$La integral es cíclica, así que se despeja para obtener su resultado: $$ \begin{align*} \int\!e^{-x}\sin(4x)\,dx &= -e^{-x}\sin(4x) - 4 e^{-x}\cos(4x) - 16\int\!e^{-x}\sin(4x)\,du \\ 17\int\!e^{-x}\sin(4x)\,dx &= -e^{-x}\sin(4x) - 4 e^{-x}\cos(4x) \\ \int\!e^{-x}\sin(4x)\,dx &= \frac{-e^{-x}\sin(4x) - 4 e^{-x}\cos(4x)}{17} \\ \int\!e^{-x}\sin(4x)\,dx &= -\frac{\sin(4x) + 4\cos(4x)}{17 e^x} \end{align*} $$
Reemplazando en $F(x)$ tenemos que: $$ \begin{align*} F(x) &= x + \int\!e^{-x}\sin(4x)\,dx \\ &= x - \frac{\sin(4x) + 4\cos(4x)}{17 e^x} \end{align*} $$
Ahora podemos evaluar esta expresión para encontrar el valor de la integral definida que buscamos: $$ \begin{align*} \int_{0}^{1}\!\left(1 + e^{-x}\sin(4x)\right)\,dx &= F(1) - F(0) \\ &= 1 - \frac{\sin(4) + 4\cos(4)}{17 e} - 0 + \frac{\sin(0) + 4\cos(0)}{17 e^0} \\ &= 1 - \frac{\sin(4) + 4\cos(4)}{17 e} + \frac{4}{17} \\ &= \frac{21}{17} - \frac{\sin(4) + 4\cos(4)}{17 e} \\ &\approx 1.308250605 \end{align*} $$
Ahora, para obtener el resultado de esta integral usando la regla de Simpson compuesta con n = 6
, simplemente definimos el integrando como una función en Python y aplicamos el método simpson13
que definimos anteriormente:
from numpy import exp, sin
def f1(x):
return 1 + exp(-x)*sin(4*x)
simpson13(f1, 0.0, 1.0, 6)
Se puede observar que el resultado por integración numérica es igual al obtenido por integración analítica hasta el tercer decimal. Si el número de subintervalos (n) fuese mayor, el resultado sería más preciso. El error porcentual de esta estimación es de:
$$ E_\% = \left|\frac{1.3084732793517873 - 1.308250605}{1.308250605}\right|\cdot100 \approx 0.017\% $$Para fines ilustrativos, a continuación se presenta la gráfica de la función y su área de integración. Primero creamos una función general para graficar el área bajo una curva y luego la aplicamos:
from matplotlib.patches import Polygon
from numpy import linspace
import matplotlib.pyplot as plt
def graficar(f, a, b, xmin, xmax):
x = linspace(xmin, xmax)
y = f(x)
fig, ax = plt.subplots()
ax.plot(x, y)
ix = linspace(a, b)
iy = f(ix)
verts = [(a, 0.0), *zip(ix, iy), (b, 0.0)]
poly = Polygon(verts, facecolor='0.9', edgecolor='0.5')
ax.add_patch(poly)
ax.axhline(0.0, color='k')
ax.axvline(0.0, color='k')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()
graficar(f1, 0.0, 1.0, -0.5, 2.0)
Aproxime las siguientes integrales usando la regla de Simpson y trapecio para n = 10
en cada caso.
Primero hallaremos la integral definida de manera analítica: $$ \begin{align*} \int_{-1}^{1}\!\left(1 + x^2\right)^{-1}\,dx &= \int_{-1}^{1}\!\frac{1}{1 + x^2}\,dx \\ &= \left[\arctan(x)\right|_{-1}^{1} \\ &= \arctan(1) - \arctan(-1) \\ &= \frac{\pi}{4} + \frac{\pi}{4} \\ &= \frac{\pi}{2} \\ &\approx 1.570796327 \end{align*} $$
Ahora, para obtener el resultado de esta integral usando la regla de Simpson compuesta con n = 10
, simplemente definimos el integrando como una función en Python y aplicamos el método simpson13
que definimos anteriormente. De igual manera, para aplicar la regla del trapecio compuesta usamos la función trapecio_compuesto
:
def f2a(x):
return (1 + x**2)**-1
print(f'Simpson:\t{simpson13(f2a, -1.0, 1.0, 10)}')
print(f'Trapecio:\t{trapecio_compuesto(f2a, -1.0, 1.0, 10)}')
Se puede observar que el resultado por la regla de Simpson es igual al obtenido por integración analítica hasta el quinto decimal. El error porcentual de esta estimación es de:
$$ E_\% = \left|\frac{1.5707953880911876 - 1.570796327}{1.570796327}\right|\cdot100 \approx 0,00000059\% $$En cambio, el método del trapecio es mucho menos preciso, pues solo es igual al valor real hasta el primer decimal. El error porcentual de esta estimación es de:
$$ E_\% = \left|\frac{1.567463056905495 - 1.570796327}{1.570796327}\right|\cdot100 \approx 0,21\% $$Para fines ilustrativos, a continuación se presenta la gráfica de la función y su área de integración:
graficar(f2a, -1.0, 1.0, -2.0, 2.0)
Primero hallaremos la integral definida de manera analítica: $$\int_{0}^{1}\!\left(2 + \sin\left(2\sqrt{x}\right)\right)\,dx = F(1) - F(0)$$ $$ \begin{align*} F(x) &= \int\!\left(2 + \sin\left(2\sqrt{x}\right)\right)\,dx \\ &= \int\!2\,dx + \int\!\sin\left(2\sqrt{x}\right)\,dx \\ &= 2x + \int\!\sin\left(2x^{\frac{1}{2}}\right)\,dx \end{align*} $$
A la integral de la derecha se le aplica cambio de variable: $$ \begin{align*} w &= 2x^{\frac{1}{2}} \\ dw &= 2\cdot\frac{1}{2}\cdot x^{\frac{1}{2} - 1}\,dx = x^{-\frac{1}{2}}\,dx \\ \Rightarrow dx &= x^{\frac{1}{2}}\,dw = \frac{w}{2}\,dw \end{align*} $$
$$ \int\!\sin\left(2x^{\frac{1}{2}}\right)\,dx = \int\!\frac{w}{2}\sin(w)\,dw = \frac{1}{2} \int\!w\sin(w)\,dw $$Luego se aplica integración por partes, derivando $w$ e integrando la función seno: $$ \begin{align*} u = w &\Rightarrow du = dw \\ dv = \sin(w) dw &\Rightarrow v = -\cos(w) \end{align*} $$
$$ \begin{align*} \frac{1}{2} \int\!w\sin(w)\,dw &= \frac{1}{2} \left[uv - \int\!v\,du\right] \\ &= \frac{1}{2} \left[-w\cos(w) - \int\!\left(-\cos(w)\right)\,dw\right] \\ &= \frac{1}{2} \left[-w\cos(w) + \int\!\cos(w)\,dw\right] \\ &= \frac{1}{2} \left[-w\cos(w) + \sin(w)\right] \end{align*} $$Devolviendo el cambio de variable tenemos que: $$ \begin{align*} \int\!\sin\left(2x^{\frac{1}{2}}\right)\,dx &= \frac{1}{2} \int\!w\sin(w)\,dw \\ &= \frac{1}{2} \left[-w\cos(w) + \sin(w)\right] \\ &= \frac{1}{2} \left[-2x^{\frac{1}{2}}\cos\left(2x^{\frac{1}{2}}\right) + \sin\left(2x^{\frac{1}{2}}\right)\right] \\ &= -\sqrt{x}\cos\left(2\sqrt{x}\right) + \frac{1}{2}\sin\left(2\sqrt{x}\right) \end{align*} $$
Reemplazando en F(x) para obtener la integral indefinida: $$ \begin{align*} F(x) &= 2x + \int\!\sin\left(2x^{\frac{1}{2}}\right)\,dx \\ &= 2x - \sqrt{x}\cos\left(2\sqrt{x}\right) + \frac{1}{2}\sin\left(2\sqrt{x}\right) \end{align*} $$
Para obtener el valor de la integral definida que buscamos, se evalúa esta expresión en los límites de integración: $$ \begin{align*} \int_{0}^{1}\!\left(2 + \sin\left(2\sqrt{x}\right)\right)\,dx &= \begin{aligned}[t] &F(1) - F(0) \end{aligned}\\ &= \begin{aligned}[t] &2\cdot1 - \sqrt{1}\cos\left(2\sqrt{1}\right) + \frac{1}{2}\sin\left(2\sqrt{1}\right) \\ &- 2\cdot0 + \sqrt{0}\cos\left(2\sqrt{0}\right) - \frac{1}{2}\sin\left(2\sqrt{0}\right) \end{aligned}\\ &= \begin{aligned}[t] &2 - \cos(2) + \frac{1}{2}\sin(2) - 0 + 0 - \frac{1}{2}\sin(0) \end{aligned}\\ &= \begin{aligned}[t] &2 - \cos(2) + \frac{1}{2}\sin(2) \end{aligned}\\ &\approx \begin{aligned}[t] &2.87078555 \end{aligned}\\ \end{align*} $$
Ahora, para obtener el resultado de esta integral usando la regla de Simpson compuesta con n = 10
, simplemente definimos el integrando como una función en Python y aplicamos el método simpson13
que definimos anteriormente. De igual manera, para aplicar la regla del trapecio compuesta usamos la función trapecio_compuesto
:
from numpy import sqrt
def f2b(x):
return 2.0 + sin(2*sqrt(x))
print(f'Simpson:\t{simpson13(f2b, 0.0, 1.0, 10)}')
print(f'Trapecio:\t{trapecio_compuesto(f2b, 0.0, 1.0, 10)}')
Se puede observar que ambos resultados son iguales al valor real hasta la primera cifra decimal, aunque sus errores no son iguales.
El error porcentual de la estimación por la regla de Simpson es de:
$$ E_\% = \left|\frac{2.8656007191058452 - 2.87078555}{2.87078555}\right|\cdot100 \approx 0,18\% $$En cambio, el método del trapecio produce un error porcentual de:
$$ E_\% = \left|\frac{2.8574088477106523 - 2.87078555}{2.87078555}\right|\cdot100 \approx 0,47\% $$Para fines ilustrativos, a continuación se presenta la gráfica de la función y su área de integración:
graficar(f2b, 0.0, 1.0, 0.0, 2.0)
Primero hallaremos la integral definida de manera analítica: $$\int_{0}^{2}\!x^2\sin\left(x^3\right)\,dx = F(2) - F(0) \\ F(x) = \int\!x^2\sin\left(x^3\right)\,dx $$
Esta integral indefinida se resuelve haciendo un cambio de variable: $$ \begin{align*} u &= x^3 \\ du &= 3x^2 \Rightarrow x^2 = \frac{1}{3} du \end{align*} $$ $$ \begin{align*} F(x) &= \int\!x^2\sin\left(x^3\right)\,dx \\ &= \int\!\frac{1}{3}\sin(u)\,du \\ &= \frac{1}{3}\int\!\sin(u)\,du \\ &= -\frac{1}{3}\cos(u) \\ &= -\frac{1}{3}\cos\left(x^3\right) \end{align*} $$
Para obtener el valor de la integral definida que buscamos, se evalúa $F(x)$ en los límites de integración: $$ \begin{align*} \int_{0}^{2}\!x^2\sin\left(x^3\right)\,dx &= F(2) - F(0) \\ &= -\frac{1}{3}\cos\left(2^3\right) + \frac{1}{3}\cos\left(0^3\right) \\ &= -\frac{1}{3}\cos(8) + \frac{1}{3}\cos(0) \\ &= -\frac{1}{3}\cos(8) + \frac{1}{3} \\ &\approx 0.3818333446 \end{align*} $$
Ahora, para obtener el resultado de esta integral usando la regla de Simpson compuesta con n = 10
, simplemente definimos el integrando como una función en Python y aplicamos el método simpson13
que definimos anteriormente. De igual manera, para aplicar la regla del trapecio compuesta usamos la función trapecio_compuesto
:
def f2c(x):
return x**2 * sin(x**3)
print(f'Simpson:\t{simpson13(f2c, 0.0, 2.0, 10)}')
print(f'Trapecio:\t{trapecio_compuesto(f2c, 0.0, 2.0, 10)}')
Se puede observar que el resultado por la regla de Simpson es mucho menos preciso que el del trapecio en este caso, pues desde el primer decimal es diferente. Esto probablemente sea porque la función posee oscilaciones relativamente bruscas.
El error porcentual de la estimación por la regla de Simpson es de:
$$ E_\% = \left|\frac{0.2888112633521003 - 0.3818333446}{0.3818333446}\right|\cdot100 \approx 24.36\% $$El error porcentual de la estimación por la regla del trapecio es de:
$$ E_\% = \left|\frac{0.38019406976081865 - 0.3818333446}{0.3818333446}\right|\cdot100 \approx 0,43\% $$Para fines ilustrativos, a continuación se presenta la gráfica de la función y su área de integración:
graficar(f2c, 0.0, 2.0, 0.0, 2.2)
Desarrolle un código usando la implementación de la regla de Simpson que acá se proporciona, que le permita generar la siguiente tabla de distribución normal y estimar los valores faltantes de la misma. Utilice la librería pandas para generar la tabla como un dataframe.
El código requerido consistirá en generar iterativamente el área bajo la curva de distribución normal estándar entre la media y valores de x positivos con respecto a ella.
La distribución normal es una distribución de probabilidad contínua para variables aleatorias de valor real. Está descrita por una función que toma como parámetros el valor x, el promedio $\mu$ y la desviación estándar $\sigma$, la cual tiene la siguiente fórmula: $$ f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}} $$
Para este ejercicio se nos pide usar la distribución normal estándar, que es el caso más simple de esta distribución, en el cual $\mu = 0$ y $\sigma = 1$. Por lo tanto, está descrita por la siguiente función: $$ \varphi (x)={\frac {1}{\sqrt {2\pi }}}e^{-{\frac {1}{2}}x^{2}} $$
Esta función es la que debemos integrar usando la regla de Simpson compuesta (usamos un valor de n = 10
). Como se mencionó anteriormente, se deben generar iterativamente los valores de la tabla, tomando como límite inferior el cero, y como límite inferior la suma del valor correspondiente de la primera columna con el valor correspondiente de la primera fila. Por ejemplo, para hallar el último valor, se integraría entre 0 y 3.09 (3.0 + 0.9).
Una vez obtenidos estos valores, se redondean a cuatro cifras decimales y se almacenan en un dataframe de la librería Pandas. Esto es de utilidad porque Jupyter puede mostrar los dataframes como tablas automáticamente.
A continuación se presenta la implementación de este algoritmo en Python:
from math import pi
import pandas as pd
def phi(x):
return exp(-0.5 * x**2) / sqrt(2 * pi)
def generar_tabla(f):
columnas = [round(i * 0.01, 2) for i in range(10)]
filas = [round(j * 0.1, 1) for j in range(31)]
datos = []
for fila in range(31):
xs = []
for columna in range(10):
val_fila = filas[fila]
val_columna = columnas[columna]
limite_sup = round(val_fila + val_columna, 2)
x = simpson13(f, 0.0, limite_sup, 10)
xs.append(round(x, 4))
datos.append(xs)
return pd.DataFrame(datos, index=filas, columns=columnas)
generar_tabla(phi)
Con esta tabla podemos determinar los valores que faltan en la imagen: $$ \begin{align*} x = 0.23 &\Rightarrow \varphi(0.23) = 0.0910 \\ x = 1.24 &\Rightarrow \varphi(1.24) = 0.3925 \\ x = 1.25 &\Rightarrow \varphi(1.25) = 0.3944 \\ x = 1.40 &\Rightarrow \varphi(1.40) = 0.4192 \\ x = 1.50 &\Rightarrow \varphi(1.50) = 0.4332 \\ x = 1.60 &\Rightarrow \varphi(1.60) = 0.4452 \\ x = 2.49 &\Rightarrow \varphi(2.49) = 0.4936 \\ x = 2.80 &\Rightarrow \varphi(2.80) = 0.4974 \\ x = 2.83 &\Rightarrow \varphi(2.83) = 0.4977 \end{align*} $$
Para concluir, utilizaremos la función graficar
que definimos al principio para mostrar el área entre la media y $x = 2.24$, como en la imagen original del ejercicio:
graficar(phi, 0.0, 2.24, -3.09, 3.09)