Integración en Python
Contents
5.6. Integración en Python#
En esta sección aprenderemos a integrar con SymPy y a aproximar integrales definidas usando NumPy.
5.6.1. Integración con SymPy#
Para calcular exactamente la integral de una función mediante SymPy, se emplea la función integrate. Por ejemplo, para calcular una primitiva de \(\sin(x)\), escribiremos
import sympy as sp
x = sp.symbols('x')
f = sp.sin(x)
I = sp.integrate(f,x)
print('Una primitiva de ',f, ' es = ',I)
Una primitiva de sin(x) es = -cos(x)
Para calcular la integral definida \(\displaystyle\int_0^\pi\sin(x)\,dx\), escribiremos
import sympy as sp
x = sp.symbols('x')
f = sp.sin(x)
Idef = sp.integrate(f,(x,0,sp.pi))
print('La integral de ',f, ' entre 0 y pi es = ',Idef)
La integral de sin(x) entre 0 y pi es = 2
SymPy no siempre es capaz de calcular una primitiva. En caso de no poder hacerlo, devuelve como salida la integral de partida:
I = sp.integrate(sp.sin(x*sp.cos(x)),x)
print(I)
Integral(sin(x*cos(x)), x)
Es posible calcular algunas integrales impropias, cuando los límites de integración son \(-\infty\) y/o \(+\infty\), es decir, integrales de la forma:
x = sp.symbols('x')
f = sp.exp(x)
I1 = sp.integrate(f,(x,-sp.oo,0))
print('Integral de ',f,' entre -oo y 0 es = ',I1)
g = 1/(x**2)
I2 = sp.integrate(g,(x,1,sp.oo))
print('Integral de ',g,' entre 1 y +oo es = ',I2)
Integral de exp(x) entre -oo y 0 es = 1
Integral de x**(-2) entre 1 y +oo es = 1
Si necesitamos calcular una integral impropia de segunda especie, como por ejemplo,
y hacemos:
h = 1/x
I2e = sp.integrate(h,(x,-1,2))
print('La integral vale = ', I2e)
La integral vale = nan
En este caso, para entender qué está ocurriendo, es mejor usar la definición:
donde
y
h = 1/x
a,b = sp.symbols('a:b', real=True)
L1 = sp.limit(sp.integrate(h,(x,-1,b)),b,0,'-')
L2 = sp.limit(sp.integrate(h,(x,a,2)),a,0,'+')
print('Valor del primer límite = ', L1)
print('Valor del segundo límite = ', L2)
print('La integral impropia es = ',L1+L2)
Valor del primer límite = -oo
Valor del segundo límite = oo
La integral impropia es = nan
5.6.2. Integración numérica. Fórmulas simples#
Como hemos visto, no siempre es posible calcular exactamente la integral de una función mediante SymPy. También puede ocurrir que la expresión de la primitiva sea demasiado costosa de evaluar o que solo conozcamos los valores de la función en un conjunto finito de puntos. En estos casos, se emplean técnicas de integración numérica.
Importamos a continuación los módulos necesarios.
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
En la Sección Integración numérica hemos visto varias fórmulas de integración numérica para aproximar una integral de la forma
donde \(a\) y \(b\) son números reales. Concretamente, hemos visto las fórmulas del punto medio, del trapecio y de Simpson.
Vamos a aproximar la integral de
mediante diferentes fórmulas de cuadratura.
Recuerda que, en este caso sencillo, podemos conocer el valor exacto de la integral:
x = sp.Symbol('x', real = True)
#Aproximación por punto medio
a = 0
b = 3
pm = (a+b)/2
f_exp = x**4+1
f = sp.lambdify(x,f_exp)
fpm = f(pm)
I_aprox = (b-a) * fpm
print('Valor aproximado de I mediante la fórmula del punto medio = ', I_aprox)
Valor aproximado de I mediante la fórmula del punto medio = 18.1875
Implementamos a continuación una función para la fórmula del trapecio simple:
def trapecio(a,b,fa,fb):
aprox_tr = (b-a) * (fa + fb)/2
return aprox_tr
La empleamos para aproximar \(I\):
x = sp.Symbol('x', real = True)
a = 0
b = 3
f_exp = x**4+1
f = sp.lambdify(x,f_exp)
fa = f(a)
fb = f(b)
aproximacion_trapecio = trapecio(a,b,fa,fb)
print('Valor aproximado de la integral por trapecio simple = ', aproximacion_trapecio)
Valor aproximado de la integral por trapecio simple = 124.5
5.6.3. Ejercicio 1#
Implementa una función que permita aproximar la integral definida de una función dada \(f\) en un intervalo \([a,b]\) mediante la fórmula de Simpson:
Empléala para aproximar el valor de \(I\).
# Escribe aquí tu código
5.6.4. Integración numérica. Fórmulas compuestas#
Como habrás podido ver en los apartados anteriores, las fórmulas simples pueden dar resultdos bastante… espantosos.
Para obtener una aproximación mejor se suelen usar fórmulas compuestas, como ya contamos en la Sección Integración numérica. Vamos a implementar ahora de manera eficiente estas fórmulas utilizando la función de ‘np.sum`.
x = sp.Symbol('x', real = True)
f_exp = x**4+1
f = sp.lambdify(x,f_exp)
a = 0; b = 3
n = 100
x1 = np.linspace(a,b,n+1) # aquí guardamos los x_{i}.
# Recuerda que, en Python, se guarda x1[0], x1[1], ..., x1[(n+1)-1] = x1[n]
y1 = f(x1)
h = (b-a)/n # el tamaño de cada subintervalo
aprox_trap = h/2 * (y1[0]+2*np.sum(y1[1:n])+y1[n])
aprox_medio = 2*h * np.sum(y1[1:n:2])
aprox_simpson = 2*h/6 * (y1[0] + 4*np.sum(y1[1:n:2])+2*np.sum(y1[2:n-1:2])+y1[n])
print('aprox_trap: ',aprox_trap)
print('aprox_medio: ',aprox_medio)
print('aprox_simpson: ',aprox_simpson)
print('Exacta: ',b**5/5+b)
aprox_trap: 51.608099919
aprox_medio: 51.583801133999984
aprox_simpson: 51.60000032399999
Exacta: 51.6
5.6.5. Ejercicio 2#
Escribe dos funciones para la programación de trapecio compuesto:
Una primera función haga el cálculo con
np.sum, como tienes ahí arriba.Una segunda función que haga lo mismo, pero con un sumatorio acumulado (será menos eficiente).
# Aquí tu código
5.6.6. Ejercicio 3#
Para obtener el volumen de un cono circular recto de altura \(h\) y radio de la base \(R\), se puede hacer girar la recta que pasa por el origen y el punto \((h,R)\) alrededor del eje \(OX\) entre \(x=0\) y \(x=h\).
Representamos a continuación un cono circular recto de altura \(h=10\) y radio de la base \(R=2\):
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
RR = 2.0
hh = 10.0
u = np.linspace(0, 6.5, 60)
v = np.linspace(0, 6.5, 60)
U, V = np.meshgrid(u, v)
X = U
Y1 = RR/hh*U*np.cos(V)
Z1 = RR/hh*U*np.sin(V)
ax.plot_surface(X, Y1, Z1, alpha=0.3, color='red', rstride=6, cstride=12)
plt.show()
Prepara un código que, utilizando la función
sp.integrate, calcule el volumen de un cono de altura \(h\) y radio de la base \(R\).Para el cono del dibujo anterior (\(h=10\), \(R=2\)), aproxima este volumen utilizando una punto medio, trapecio y Simpson compuestas con \(n=1000\).
# Tu código aquí