Referencia:
En análisis de ingeniería, normalmente debemos evaluar integrales definidas sobre un dominio complejo o en un espacio de dimensión alta.
Por ejemplo, podríamos querer calcular:
- la deflexión en una viga de geometría complicada,
- el volumen de una parte tridimensional de una aeronave,
- o evaluar alguna medida de rendimiento (rentabilidad) en algún proceso que sea expresada como una integral de alguna función sin antiderivada primitiva (que se pueda expresar en términos de funciones elementales).
A la mano tenemos herramientas de integración analítica cuando tanto el espacio de integración como la función a integrar son simples. Cuando la función a integrar es difícil (incluso, imposible) de integrar podemos aún recurrir a métodos numéricos de integración.
Desafortunadamente, los métodos determinísiticos de integración fallan cuando:
- la región es demasiado compleja para discretizarla,
- o la función a integrar es demasiado irregular,
- o la convergencia es demasiado lenta debido a la alta dimensionalidad del espacio de integración (ver Maldición de la dimensionalidad).
Por eso en esta clase veremos una técnica alternativa de integración numérica: Integración Montecarlo.
De su curso de cálculo integral seguro recordarán (o estarán viendo) que existen funciones cuya integral no tiene primitiva. Es decir, que no podemos encontrar una función que se pueda expresar en forma de funciones elementales cuya derivada sea tal función.
Esto no significa que dicha función no se pueda integrar, ya que sabemos que cualquier función continua es integrable (y la mayoría de funciones que vemos a ese nivel, lo son). Lo que ocurre es que no podemos expresar dicha integral de una forma sencilla (por ejemplo, en función de exponenciales, senos, cosenos, logaritmos...).
Algunas integrales que no son elementales son:
Referencia:
In [41]:
from IPython.display import YouTubeVideo
YouTubeVideo('Ti5zUD08w5s')
Out[41]:
In [2]:
YouTubeVideo('jmsFC0mNayM')
Out[2]:
Se basa en la definición de valor promedio de una función y en el valor esperado de una variable aleatoria uniforme.
Presentamos esto mediante un ejemplo.
Ejemplo. Aproxime el área bajo la curva $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
In [1]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
In [2]:
def parab(x):
return x**2
In [3]:
x = np.linspace(0,1)
y = parab(x)
plt.fill_between(x,y)
plt.text(0.8,0.2,'$\mathcal{D}$',fontsize=20)
plt.show()
Entonces, lo que queremos es aproximar el área de la región $\mathcal{D}$. Llamaremos esta área $A(\mathcal{D})$.
De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}y\text{d}x=\int_{0}^{1}x^2\text{d}x$$.
Por definición, el valor promedio de una función $f:\left[a,b\right]\to\mathbb{R}$ en un intervalo $\left[a,b\right]$ es
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x.$$Entonces, el área bajo la curva $y=x^2$ es exactamente el valor promedio de $f(x)=x^2$ en $\left[0,1\right]$. Este valor promedio puede aproximarse mediante el promedio de los valores de la función en puntos aleatorios uniformemente distribuidos en el intervalo $\left[0,1\right]$. Es decir,
$$A(\mathcal{D})=\int_{0}^{1}x^2\text{d}x=\int_{0}^{1}f(x)\text{d}x\approx \frac{1}{N}\sum_{i=1}^{N}f(u_i)=\frac{1}{N}\sum_{i=1}^{N}u_i^2$$,
donde $u_i$ son realizaciones de la variable aleatoria $U\sim\mathcal{U}\left[0,1\right]$ ($U$ distribuye uniformemente en el intervalo $\left[0,1\right]$).
¿Cómo construit vectores de números aleatorios?
En este caso necesitamos $N$ números aleatorios uniformemente distribuidos...
In [4]:
help(np.random.uniform)
In [17]:
N = 100000
x = np.random.uniform(0, 1, N)
A_Dapprox = np.sum(parab(x))/N
A_Dapprox
Out[17]:
En este caso, la integral se puede hacer fácilmente. Comparemos el resultado con el valor real:
$$A(\mathcal{D})=\int_{0}^{1}x^2\text{d}x=\left.\frac{x^3}{3}\right|_{x=0}^{x=1}=\frac{1}{3}$$Hagamos una tabla viendo:
In [18]:
import pandas as pd
In [23]:
A_D = 1/3
N = np.logspace(1,7,7)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion', 'Error_relativo'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
x = np.random.uniform(0, 1, n.astype(int))
df.loc[n,"Valor_aproximacion"] = np.sum(parab(x))/n
df.loc[n,"Error_relativo"] = np.abs(df.loc[n,"Valor_aproximacion"]-A_D)/A_D
df
Out[23]:
Ver que los resultados son distintos cada vez (¿porqué?). Sin embargo, se aproximan más o menos en la misma medida.
Aproximación de integrales en intervalos distintos a $\left[0,1\right]$.
Sin embargo, no todas las integrales que hacemos son en el intervalo $\left[0,1\right]$. En general, podemos integrar cualquier función continua en el intervalo $\left[a,b\right]$, donde $a,b\in\mathbb{R}$ con $a<b$.
Sea $f:\left[a,b\right]\to\mathbb{R}$ una función continua en el intervalo $\left(a,b\right)$ (por lo tanto es integrable endicho intervalo). Queremos resolver:
$$\int_{a}^{b}f(x)\text{d}x.$$¿Cómo podemos usar la idea del valor promedio para resolver esto?
El valor promedio de $f$ en $\left[a,b\right]$ es:
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x.$$Este valor promedio puede aproximarse mediante el promedio de $N$ valores de la función en puntos aleatorios uniformemente distribuidos en el intervalo $\left[a,b\right]$. Es decir,
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x\approx \frac{1}{N}\sum_{i=1}^{N}f(u_i)$$,
donde $u_i$ son realizaciones de la variable aleatoria $U\sim\mathcal{U}\left[a,b\right]$ ($U$ distribuye uniformemente en el intervalo $\left[a,b\right]$).
Finalmente, la aproximación montecarlo tipo 1 con $N$ términos es
$$\int_{a}^{b}f(x)\text{d}x\approx \frac{b-a}{N}\sum_{i=1}^{N}f(u_i)$$,
Escribamos una función que tenga como entradas:
y que devuelva la aproximación montecarlo tipo 1 de la integral $\int_{a}^{b}f(x)\text{d}x$.
In [26]:
# Escribir la función acá
def int_montecarlo1(f, a, b, N):
return (b-a)/N*np.sum(f(np.random.uniform(a,b,N)))
Actividad. Utilizar la anterior función para realizar las siguientes integrales. Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos) y cuyas columnas correspondan a las funciones.
In [27]:
# Resolver
def func1(x):
return np.exp(x**2)
def func2(x):
return 1/np.log(x)
def func3(x):
return np.sin(x)/x
In [29]:
a, b = 4, 5
N = np.logspace(1,5,5)
df = pd.DataFrame(index=N,columns=['Funcion1', 'Funcion2', 'Funcion3'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
df.loc[n,"Funcion1"] = int_montecarlo1(func1, a, b, n.astype(int))
df.loc[n,"Funcion2"] = int_montecarlo1(func2, a, b, n.astype(int))
df.loc[n,"Funcion3"] = int_montecarlo1(func3, a, b, n.astype(int))
df
Out[29]:
Con la integración montecarlo tipo 1 pudimos aproximar integrales de funciones continuas de una variable en un intervalo dado. En realidad este mismo análisis se puede ampliar para aproximar integrales definidas de funciones continuas de varias variables (integrales sobre áreas, volúmenes e hipervolúmenes) dado que la noción de valor promedio de una función se extiende a cualquier dimensión.
Este es en realidad el caso interesante, pues las integrales de funciones complicadas también se pueden aproximar por métodos numéricos clásicos, pero cuando la dimensión aumenta es cuando montecarlo se vuelve una herramienta relevante. Dado que no lo veremos en clase por la limitación de que la mayoría no han visto cálculo en varias variables, este tema puede ser elegido como proyecto de módulo, donde se exploraría también como mejorar la aproximación de integrales montecarlo.
Como vimos en el ejemplo (y como debe ser claro de su curso de cálculo integral) una de las aplicaciones más importantes de la integración es hallar áreas. Y no solo el área bajo una curva, sino áreas entre curvas y áreas de regiones más complicadas.
Antes de ver la integración montecarlo tipo 2, ¿cómo podemos usar la integración montecarlo tipo 1 para aproximar el área entre curvas?
Ejemplo. Aproxime el área entre las curvas $y=x$, y $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
In [44]:
x = np.linspace(-0.1,1.1)
y = parab(x)
plt.plot(x,x,'k--',label='$y=x$')
plt.plot(x,y,'k',label='$y=x^2$')
plt.fill_between(x,x,y)
plt.text(0.5,0.4,'$\mathcal{D}$',fontsize=20)
plt.legend(loc='best')
plt.show()
De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}x-x^2\text{d}x.$$Entonces...
In [31]:
# Usar la funcion int_montecarlo1
def f(x):
return x-x**2
A_Daprox = int_montecarlo1(f, 0, 1, 100000000)
A_Daprox
Out[31]:
De modo que si la región se puede describir fácilmente, diría el ferras 'no hay pedo, lo pago' (podemos usar montecarlo tipo 1).
In [42]:
YouTubeVideo('G8fOTMYDPEA')
Out[42]:
Pero, ¿qué pasa si la geometría de la región no se puede describir fácilmente?
Como en el caso anterior, motivaremos el método con un caso conocido. Vamos a aproximar el valor de $\pi$ usando el área de un círculo unitario.
Dibujemos el círculo unitario en la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$.
In [33]:
def circ_arriba(x, r):
return np.sqrt(r**2-x**2)
def circ_abajo(x, r):
return -np.sqrt(r**2-x**2)
In [36]:
x = np.linspace(-1,1,100)
y1 = circ_arriba(x, 1)
y2 = circ_abajo(x, 1)
plt.figure(figsize=(5,5))
plt.plot(x,y1,'k')
plt.plot(x,y2,'k')
plt.fill_between(x,y1,y2)
plt.text(0,0,'$\mathcal{D}$',fontsize=20)
plt.text(0.8,0.8,'$\mathcal{R}$',fontsize=20)
plt.show()
Si aproximamos $A(\mathcal{D})$ aproximamos el valor de $\pi$, pues el área del círculo unitario es:
$$A(\mathcal{D})=\pi(1)^2=\pi.$$Por otra parte es claro que el área de la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$ es
$$A(\mathcal{R})=4.$$Ahora, haremos uso de nuestro generador de números aleatorios. Supongamos que escogemos un punto aleatorio en la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$. Describimos este punto como $(X,Y)$ para $X$ e $Y$ variables aleatorias uniformes sobre el intervalo $\left[-1,1\right]$.
¿Cómo generamos puntos aleatorios en un rectángulo?
In [50]:
N = 1000000
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
In [43]:
X, Y = np.meshgrid(x,y)
plt.figure(figsize=(5,5))
plt.scatter(X,Y)
plt.show()
La probabilidad de que el punto $(X,Y)$ esté en el círculo unitario $\mathcal{D}$ es
$$P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}=\frac{\pi}{4}.$$Luego, definimos una variable aleatoria de Bernoulli $B$ de manera que
$$B=\left\lbrace\begin{array}{ccc}0 & \text{si} & (X,Y)\notin\mathcal{D}\\1 & \text{si} & (X,Y)\in\mathcal{D} \end{array}\right.=\left\lbrace\begin{array}{ccc}0 & \text{si} & X^2+Y^2>1\\1 & \text{si} & X^2+Y^2\leq 1 \end{array}\right..$$Entonces, el valor esperado de la variable aleatoria $B$ es
$$E\left[B\right]=\theta=P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}.$$De lo anterior, una estimación de theta se puede obtener como
$$\theta=\frac{A(\mathcal{D})}{A(\mathcal{R})}\approx \frac{1}{N}\sum_{i=1}^{N}b_i,$$donde
$$b_i=\left\lbrace\begin{array}{ccc}0 & \text{si} & x_i^2+y_i^2>1\\1 & \text{si} & x_i^2+y_i^2\leq 1 \end{array}\right.$$son realizaciones de la variable aleatoria $B$, que a su vez es producto de las realizaciones $x_i$ e $y_i$ de las variables aleatorias $X$ e $Y$, respectivamente.
Finalmente, la aproximación montecarlo tipo 2 con $N$ términos es
$$A(\mathcal{D})\approx \frac{A(\mathcal{R})}{N}\sum_{i=1}^{N}b_i.$$
In [46]:
def reg_circ(x,y):
return x**2+y**2<=1
In [52]:
A_R = 4
A_Dapprox = A_R*np.sum(reg_circ(x,y))/N
A_Dapprox
Out[52]:
De nuevo, comparemos con el valor real
In [53]:
A_D = np.pi
N = np.logspace(1,7,7)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion', 'Error_relativo'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
x = np.random.uniform(-1, 1, n.astype(int))
y = np.random.uniform(-1, 1, n.astype(int))
df.loc[n,"Valor_aproximacion"] = A_R*np.sum(reg_circ(x,y))/n
df.loc[n,"Error_relativo"] = np.abs(df.loc[n,"Valor_aproximacion"]-A_D)/A_D
df
Out[53]:
Escribamos una función que tenga como entradas:
y que devuelva la aproximación montecarlo tipo 2 del area de la region.
In [54]:
# Escribir la función acá
def int_montecarlo2(region, a1, b1, a2, b2, N):
A_R = (b1-a1)*(b2-a2)
x = np.random.uniform(a1, b1, N.astype(int))
y = np.random.uniform(a2, b2, N.astype(int))
return A_R*np.sum(region(x,y))/N
Actividad. Utilizar la anterior función para aproximar el área de la región descrita por
$$4(2x-1)^4+8(2y-1)^8<1+2(2y-1)^3(3x-2)^2$$Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos).
In [56]:
N = 100
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
In [57]:
def region(x,y):
return 4*(2*x-1)**4+8*(2*y-1)**8 < 1+2*(2*y-1)**3*(3*x-2)**2
In [58]:
X, Y = np.meshgrid(x,y)
plt.figure(figsize=(5,5))
plt.scatter(X,Y,c=~region(X,Y),cmap='bone')
plt.show()
In [59]:
# Resolver
a1, a2, b1, b2 = 0, 0, 1, 1
N = np.logspace(1,5,5)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
df.loc[n,"Valor_aproximacion"] = int_montecarlo2(region, a1, b1, a2, b2, n)
df
Out[59]: