El framework SciPy se basa en el marco NumPy de bajo nivel para matrices multidimensionales y proporciona un gran número de algoritmos científicos de nivel superior. Algunos de los temas que SciPy cubre son:
Cada uno de estos submódulos proporciona una serie de funciones y clases que pueden utilizarse para resolver problemas en sus respectivos temas.
Para acceder al paquete SciPy en un programa Python, empezamos por importar todo desde el módulo scipy.
from scipy import *
# Esta línea configura matplotlib para mostrar las figuras incrustadas en el jupyter notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['figure.figsize'] = (15, 8)
Evaluación numérica de una función del tipo
$$\int_{a}^{b} f(x) dx$$Se llama cuadratura numérica, o simplemente cuadratura. SciPy proporciona una serie de funciones para diferentes tipos de cuadratura, por ejemplo, quad, dblquad y tplquad para integrales individuales, dobles y triples, respectivamente.
Por ejemplo calculemos la integral definida de
$$\int_{2}^{5} x^{2} dx$$from scipy.integrate import quad, dblquad, tplquad
# define a simple function for the integrand
def f(x):
return x**2
x_lower = 2 # the lower limit of x
x_upper = 5 # the upper limit of x
val, abserr = quad(f, x_lower, x_upper)
print("El valor de la integral = {} , error absoluto = {}".format(val,abserr))
El valor de la integral = 39.00000000000001 , error absoluto = 4.3298697960381115e-13
Como se muestra en el ejemplo anterior, también podemos usar 'Inf' o '-Inf' como límites integrales.
Para las integrales dobles:
$$ \int\limits_{\vphantom{b}0}^{10}\int\limits_2^{20} e^{-x^2-y^2} \mathrm{d}x \mathrm{d}y $$def integrand(x, y):
return exp(-x**2-y**2)
x_lower = 2
x_upper = 20
y_lower = 0
y_upper = 10
val, abserr = dblquad(integrand, x_lower, x_upper, lambda x: y_lower, lambda x: y_upper)
print(val,abserr)
0.003673884462971745 1.4031303125082134e-08
SciPy proporciona dos formas diferentes de resolver EDO's: Una API basada en la función odeint, y la API orientada a objetos basada en la clase ode. Por lo general, odeint es más fácil de empezar, pero la clase ode ofrece un nivel de control más fino.
Aquí vamos a utilizar las funciones odeint. Para obtener más información acerca de la clase ode, pruebe help (ode). Hace casi lo mismo que odeint, pero de una manera orientada a objetos.
Para integrar EDOs vamos a usar la función odeint del paquete integrate.
Los problemas EDO son importantes en física computacional, por lo que vamos a mirar un ejemplo más: la oscilación armónica amortiguada. Este problema está bien descrito en la página wiki: https://es.wikipedia.org/wiki/Oscilador_arm%C3%B3nico
La ecuación de movimiento para el oscilador amortiguado es:
$$ \frac{d^2 x}{dt^2} + 2 \zeta \omega_0\frac{dx}{dt} + \omega_0^2 x = 0$$Donde $x$ es la posición del oscilador, $\omega_0$ es la frecuencia y $\zeta$ es el coeficiente de amortiguación. Para escribir la ecuación diferencial de segundo orden de forma estandar introducimos $p = \frac{dx}{dp}$
$$ \frac{dp}{dt} = -2 \zeta \omega_0 p - \omega_0^2 x $$$$\frac{dx}{dt} = p $$En la implementación de este ejemplo agregaremos argumentos adicionales a la función RHS para la EDO. Como consecuencia de los argumentos adicionales a la RHS, necesitamos pasar un argumento de palabras clave args a la función odeint:
def dy(y, t, zeta, w0):
"""
El lado derecho del oscilador amortiguado EDO
"""
x, p = y[0], y[1]
dx = p
dp = -2 * zeta * w0 * p - w0**2 * x
return [dx, dp]
from scipy.integrate import odeint, ode
# Estado inicial:
y0 = [1.0, 0.0]
# tiempo para resolver el EDO
t = linspace(0, 10, 1000)
w0 = 2*pi*1.0
# Resolviendo el problema de EDO para los 3 diferentes coeficientes de amortiguación
y1 = odeint(dy, y0, t, args=(0.0, w0)) # sin amortiguar
y2 = odeint(dy, y0, t, args=(0.2, w0)) # Subamortiguado
y3 = odeint(dy, y0, t, args=(1.0, w0)) # Amortiguamiento critico
y4 = odeint(dy, y0, t, args=(5.0, w0)) # Sobreamortiguado
# Graficando
fig, ax = plt.subplots()
ax.plot(t, y1[:,0], 'k', label="Sin amortiguar", linewidth=0.25)
ax.plot(t, y2[:,0], 'r', label="Subamortiguado")
ax.plot(t, y3[:,0], 'b', label=r"Amortiguamiento critico")
ax.plot(t, y4[:,0], 'g', label="Sobreamortiguado")
ax.legend();
El módulo de álgebra lineal contiene una gran cantidad de funciones relacionadas con matrices, incluyendo la resolución de ecuaciones lineales, los solucionadores de valores propios, las funciones de matriz (por ejemplo matriz-exponenciación), varias descomposiciones (DVS, LU, cholesky), etc.
Documentación más detallada en: http://docs.scipy.org/doc/scipy/reference/linalg.html
Aquí veremos cómo usar algunas de estas funciones:
Donde A es la matrz y x,b son los vectores
import scipy.linalg as scla
import numpy as np
A = np.array([[1,2,3], [4,5,6], [7,8,9]])
b = np.array([1,2,3])
x = scla.solve(A, b)
x
array([-0.33333333, 0.66666667, 0. ])
# Checkeemos si el resultado es correcto
np.dot(A, x) - b
array([ 0.00000000e+00, -2.22044605e-16, 0.00000000e+00])
Para hallar los valores propios de una matriz $A$ donde:
$$ Av_{n} = \lambda_{n}v_{n} $$Donde $v_{n}$ es el enésimo autovector y $\lambda_{n}$ es el enésimo autovector.
Para calcular los valores propios de una matriz, utilice eigvals y para calcular tanto los autovalores como los vectores propios, utilice la función eig:
A = np.array([[1,2,0],
[0,2,0],
[-2,-2,-1]])
evals = scla.eigvals(A)
print(evals)
[-1.+0.j 1.+0.j 2.+0.j]
evals, evecs = scla.eig(A)
print(evals)
print()
print(evecs)
[-1.+0.j 1.+0.j 2.+0.j] [[ 0. 0.70710678 0.66666667] [ 0. 0. 0.33333333] [ 1. -0.70710678 -0.66666667]]
A = np.array([[1,2,3],
[0,1,4],
[5,6,0]])
A_inv = scla.inv(A)
print(A_inv)
[[-24. 18. 5.] [ 20. -15. -4.] [ -5. 4. 1.]]
print(scla.det(A))
0.9999999999999964
# Normas matriciales
norma_2 = scla.norm(A, ord=2)
norma_inf = scla.norm(A, ord=np.Inf)
print(norma_2,norma_inf)
8.2788392305 11.0
La optimización (encontrar mínimos o máximos de una función) es un campo grande en matemáticas, y la optimización de funciones complicadas o en múltiples variables puede ser algo complicada. Aquí sólo veremos algunos casos muy simples. Para una introducción más detallada a la optimización con SciPy, consulte:
http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html
Para utilizar el módulo de optimización en scipy primero incluya el módulo optimize:
from scipy import optimize
Veamos primero cómo encontrar los mínimos de una función simple de una sola variable:
# Definimos la función:
def f(x):
return 4*x**3 + (x-2)**2 + x**4
# Graficamos
fig, ax = plt.subplots()
x = np.linspace(-5, 3, 100)
ax.plot(x, f(x));
Podemos usar la función fmin_bfgs para encontrar los mínimos de una función:
x_min = optimize.fmin_bfgs(f, -2)
x_min
Optimization terminated successfully. Current function value: -3.506641 Iterations: 5 Function evaluations: 24 Gradient evaluations: 8
array([-2.67298151])
Para encontrar la raíz de una función de la forma $f(x)= 0$ podemos usar la función fsolve. Requiere de una aproximación inicial:
# Las raices de esta función son -2,-0.25,3
def f(x):
return 4*x**3 - 3*x**2 - 25*x - 6
fig, ax = plt.subplots(figsize=(10,4))
x = np.linspace(-4, 7, 1000)
y = f(x)
ax.plot(x, y)
ax.plot([-4, 7], [0, 0], 'k')
ax.set_ylim(-50,50);
ax.set_xlim(-4,7)
ax.grid(True)
optimize.fsolve(f, -1)
array([-0.25])
optimize.fsolve(f, -3)
array([-2.])
optimize.fsolve(f, 2)
array([ 3.])
La clase interp1d en scipy.interpolate es un método conveniente para crear una función basada en puntos de datos fijos que pueden ser evaluados en cualquier lugar dentro del dominio definido por los datos dados usando interpolación lineal. Una instancia de esta clase se crea pasando los vectores 1-d que comprenden los datos.
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
x = np.linspace(0, 10, num=10, endpoint=True)
y = np.cos(-x**2/9.0)
f = interp1d(x, y) # Interpolación Lineal
f2 = interp1d(x, y, kind='cubic') # Interpolación cúbica
xnew = np.linspace(0, 10, num=41, endpoint=True)
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'Int. Lineal', 'Int. cúbica'], loc='best')
plt.show()
El módulo scipy.stats contiene un gran número de distribuciones estadísticas, funciones estadísticas y tests. Para obtener una documentación completa de sus características, consulte http://docs.scipy.org/doc/scipy/reference/stats.html.
También hay un potente paquete de python para el modelado estadístico llamado statsmodels. Vea http://statsmodels.sourceforge.net para más detalles.
from scipy import stats
# Crear una variable aleatoria (discreta) con distribución poissioniana
X = stats.poisson(3.5) # Poisson distribution for a coherent state with n=3.5 coef
n = np.arange(0,15)
# Para que nos de 3 figuras
fig, axes = plt.subplots(3,1, sharex=True)
# Grafica el PMF (probability mass function)
axes[0].step(n, X.pmf(n))
# Grafica el CDF(commulative distribution function)
axes[1].step(n, X.cdf(n))
# Grafica un histograma de 1000 realizaciones aleatorias de la variable estocástica X
axes[2].hist(X.rvs(size=1000));
# create a (continous) random variable with normal distribution
Y = stats.norm()
x = np.linspace(-5,5,100)
fig, axes = plt.subplots(3,1, sharex=True)
# plot the probability distribution function (PDF)
axes[0].plot(x, Y.pdf(x))
# plot the commulative distributin function (CDF)
axes[1].plot(x, Y.cdf(x));
# plot histogram of 1000 random realizations of the stochastic variable Y
axes[2].hist(Y.rvs(size=1000), bins=50);
#Statistics:
X.mean(), X.std(), X.var() # poission distribution
(3.5, 1.8708286933869707, 3.5)
Y.mean(), Y.std(), Y.var() # normal distribution
(0.0, 1.0, 1.0)
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../styles/StyleCursoPython.css'
HTML(open(css_file, "r").read())