¿Te acuerdas de todos esos esquemas numéricos para integrar ecuaciones diferenciales ordinarias? Es bueno saber que existen y qué peculiaridades tiene cada uno, pero en este curso no queremos implementar esos esquemas: queremos resolver las ecuaciones. Los problemas de evolución están por todas partes en ingeniería y son de los más divertidos de programar.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Como sabemos, las operaciones del álgebra lineal aparecen con mucha frecuencia a la hora de resolver sistemas de ecuaciones en derivadas parciales y en general al linealizar problemas de todo tipo, y suele ser necesario resolver sistemas con un número enorme de ecuaciones e incógnitas. Gracias a los arrays de NumPy podemos abordar este tipo de cálculos en Python, ya que todas las funciones están escritas en C o Fortran y tenemos la opción de usar bibliotecas optimizadas al límite.
El paquete de álgebra lineal en NumPy se llama linalg
, así que importando NumPy con la convención habitual podemos acceder a él escribiendo np.linalg
. Si imprimimos la ayuda del paquete vemos que tenemos funciones para:
help(np.linalg)
Help on package numpy.linalg in numpy: NAME numpy.linalg DESCRIPTION Core Linear Algebra Tools ------------------------- Linear algebra basics: - norm Vector or matrix norm - inv Inverse of a square matrix - solve Solve a linear system of equations - det Determinant of a square matrix - lstsq Solve linear least-squares problem - pinv Pseudo-inverse (Moore-Penrose) calculated using a singular value decomposition - matrix_power Integer power of a square matrix Eigenvalues and decompositions: - eig Eigenvalues and vectors of a square matrix - eigh Eigenvalues and eigenvectors of a Hermitian matrix - eigvals Eigenvalues of a square matrix - eigvalsh Eigenvalues of a Hermitian matrix - qr QR decomposition of a matrix - svd Singular value decomposition of a matrix - cholesky Cholesky decomposition of a matrix Tensor operations: - tensorsolve Solve a linear tensor equation - tensorinv Calculate an inverse of a tensor Exceptions: - LinAlgError Indicates a failed linear algebra operation PACKAGE CONTENTS _umath_linalg info lapack_lite linalg setup DATA absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0... division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192... print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)... FILE /home/siro/anaconda3/lib/python3.4/site-packages/numpy/linalg/__init__.py
El producto matricial usual (no el que se hace elemento a elemento, sino el del álgebra lineal) se calcula con la misma función que el producto matriz-vector y el producto escalar vector-vector: con la función dot
, que no está en el paquete linalg
sino directamente en numpy
.
Una consideración importante a tener en cuenta es que en NumPy no hace falta ser estricto a la hora de manejar vectores como si fueran matrices columna, siempre que la operación sea consistente. Un vector es una matriz con una sola dimensión: por eso si calculamos su traspuesta no funciona.
Juguemos un poco con ella para ver cómo funciona!
Ejercicio:
A = np.array([
[1, 2],
[2, 1],
[1, 5]
])
B = np.array([
[1, 2, 3],
[3, 2, 1]
])
v1 = np.array([1, 0, 1])
v2 = np.array([2, 2, 2])
v3 = np.array([1, 2])
print(A,'\n\n', B, '\n\n',v1, '\n\n',v2, '\n\n',v3)
[[1 2] [2 1] [1 5]] [[1 2 3] [3 2 1]] [1 0 1] [2 2 2] [1 2]
print(np.dot(A,B)) #ETC
[[ 7 6 5] [ 5 6 7] [16 12 8]]
Existe también otra función muy interesante: np.linalg.solve La usaremos para:
Resolver el siguiente sistema:
$$ \begin{pmatrix} 2 & 0 & 1 \\ -1 & 1 & 0 \\ 3 & 2 & -1 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} -1 \\ 3 \\ 0 \end{pmatrix} $$M = np.array([
[2, 0, 1],
[-1, 1, 0],
[3, 2, -1]
])
M = np.array([-1, 3, 0])
x = np.linalg.solve(M, V)
x
array([-1., 2., 1.])
np.dot(M, x)
array([-1., 3., 0.])
Cargar y guardar datos
Numpy tiene dos funciones específicas para leer matrices y guardarlas: np.loadtxt y np.savetxt
np.savetxt('array.txt', x, fmt='%.4e', header='Nuestro array')
z = np.loadtxt('array.txt')
z
array([-1., 2., 1.])
Visto cómo resolver sistemas de ecuaciones lineales, tal vez sea incluso más atractivo resolver ecuaciones no lineales. Para ello, importaremos el paquete optimize
de SciPy:
from scipy import optimize
La ayuda de este paquete es bastante larga (puedes consultarla también en http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html). El paquete optimize
incluye multitud de métodos para optimización, ajuste de curvas y búsqueda de raíces. Vamos a centrarnos ahora en la búsqueda de raíces de funciones escalares. Para más información puedes leer http://pybonacci.org/2012/10/25/como-resolver-ecuaciones-algebraicas-en-python-con-scipy/
Hay básicamente dos tipos de algoritmos para hallar raíces de ecuaciones no lineales:
De los primeros vamos a usar la función brentq
(aunque podríamos usar bisect
) y de los segundos vamos a usar newton
(que en realidad engloba los métodos de Newton y de la secante).
Ejemplo:
$\ln{x} = \sin{x} \Rightarrow F(x) \equiv \ln{x} - \sin{x} = 0$
Lo primero que tengo que hacer es definir la ecuación, que matemáticamente será una función $F(x)$ que quiero igualar a cero.
def F(x):
return np.log(x) - np.sin(x)
Para hacernos una idea de las posibles soluciones siempre podemos representar gráficamente esa función:
x = np.linspace(0, 10, num=100)
plt.plot(x, F(x), 'k', lw=2, label="$F(x)$")
plt.plot(x, np.log(x), label="$\log{x}$")
plt.plot(x, np.sin(x), label="$\sin{x}$")
plt.plot(x, np.zeros_like(x), 'k--')
plt.legend(loc=4)
-c:2: RuntimeWarning: divide by zero encountered in log -c:3: RuntimeWarning: divide by zero encountered in log
<matplotlib.legend.Legend at 0x7f47eb8c6cc0>
Y utilizando por ejemplo el método de Brent en el intervalo $[0, 3]$:
optimize.brentq(F, 0, 3)
2.219107148913746
1 / 0
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) <ipython-input-7-b710d87c980c> in <module>() ----> 1 1 / 0 ZeroDivisionError: division by zero
1 / np.array([0])
-c:1: RuntimeWarning: divide by zero encountered in true_divide
array([ inf])
Obtener por ambos métodos (newton
y brentq
) una solución a la ecuación $\tan{x} = x$ distinta de $x = 0$. Visualizar el resultado.
Nuestras funciones siempre tienen que tomar como primer argumento la incógnita, el valor que la hace cero. Si queremos incluir más, tendremos que usar el argumento args
de la funciones de búsqueda de raíces. Este patrón se usa también en otras partes de SciPy, como ya veremos.
Vamos a resolver ahora una ecuación que depende de un parámetro:
.
def G(x, C):
return C - np.sqrt(x) - np.log(x)
Nuestra incógnita sigue siendo $x$, así que debe ir en primer lugar. El resto de parámetros van a continuación, y sus valores se especifican a la hora de resolver la ecuación usando args
:
optimize.newton(G, 2.0, args=(2,))
1.8773216666875554
Esta es la relación isentrópica entre el número de Mach $M(x)$ en un conducto de área $A(x)$:
Para un conducto convergente:
$$ \frac{A(x)}{A^*} = 3 - 2 x \quad x \in [0, 1]$$Hallar el número de Mach en la sección $x = 0.9$.
def A(x):
return 3 - 2 * x
x = np.linspace(0, 1)
area = A(x)
r = np.sqrt(area / np.pi)
plt.fill_between(x, r, -r, color="#ffcc00")
<matplotlib.collections.PolyCollection at 0x7f47eb80a0b8>
¿Cuál es la función $F$ ahora? Hay dos opciones: definir una función $F_{0.9}(M)$ que me da el número de Mach en la sección $0.9$ o una función $F(M; x)$ con la que puedo hallar el número de Mach en cualquier sección. Bonus points si haces la segunda opción :)
Para resolver la ecuación utiliza el método de Brent (bisección). ¿En qué intervalo se encontrará la solución? ¡Si no te haces una idea es tan fácil como pintar la función $F$!
def F(M, x, g):
return A(x) - (1 / M) * ((2 / (1 + g)) * (1 + (g - 1) / 2 * M ** 2)) ** ((g + 1) / (2 * (g - 1)))
optimize.brentq(F, 0.01, 1, args=(0.9, 1.4))
0.5902487609888621
Representar la ecuación de Kepler
$$M = E - e \sin E$$que relaciona dos parámetros geométricos de las órbitas elípticas, la anomalía media $M$ y la anomalía excéntrica $E$.
para los siguientes valores de excentricidad:
Para reproducir esta gráfica:
from IPython.display import HTML
HTML('<iframe src="http://en.m.wikipedia.org/wiki/Kepler%27s_equation" width="800" height="400"></iframe>')
Para ello utilizaremos el método de Newton (secante).
1- Define la función correspondiente a la ecuación de Kepler, que no solo es una ecuación implícita sino que además depende de un parámetro. ¿Cuál es la incógnita?
def F(E, e, M):
return M - E + e * np.sin(E)
2- Como primer paso, resuélvela para la excentricidad terrerestre y anomalía media $M = 0.3$. ¿Qué valor escogerías como condición inicial?
optimize.newton(F, 0.3, args=(0.0167, 0.3))
0.30501513714875778
3- Como siguiente paso, crea un dominio (linspace
) de anomalías medias entre $0$ y $2 \pi$ y resuelve la ecuación de Kepler con excentricidad terrestre para todos esos valores. Fíjate que necesitarás un array donde almacenar las soluciones. Representa la curva resultante.
N = 500
M = np.linspace(0, 2 * np.pi, N)
sol = np.zeros_like(M)
for ii in range(N):
sol[ii] = optimize.newton(F, sol[ii - 1], args=(0.249, M[ii]))
plt.plot(M, sol)
[<matplotlib.lines.Line2D at 0x7f47eb7812e8>]
4- Como último paso, solo tienes que meter parte del código que ya has escrito en un bucle que cambie el valor de la excentricidad 5 veces. Es aconsejable que tengas todo ese código en una única celda (esta de aquí abajo).
Vamos a introducir aquí un truco muy útil en Python:
M = np.linspace(0, 2 * np.pi, N)
sol = np.zeros_like(M)
plt.figure(figsize=(6, 6))
for ee in 0.0167, 0.249, 0.432, 0.775, 0.967:
# Para cada valor de excentricidad sobreescribimos el array sol
for ii in range(N):
sol[ii] = optimize.newton(F, sol[ii - 1], args=(ee, M[ii]))
plt.plot(M, sol)
plt.xlim(0, 2 * np.pi)
plt.ylim(0, 2 * np.pi)
plt.xlabel("$M$", fontsize=15)
plt.ylabel("$E$", fontsize=15)
plt.gca().set_aspect(1)
plt.grid(True)
plt.legend(["Earth", "Pluto", "Comet Holmes", "28P/Neujmin", "Halley's Comet"], loc=2)
plt.title("Kepler's equation solutions")
<matplotlib.text.Text at 0x7f47eb737d68>
Para integrar EDOs vamos a usar la función odeint
del paquete integrate
, que permite integrar sistemas del tipo:
con condiciones iniciales $\mathbf{y}(\mathbf{0}) = \mathbf{y_0}$.
from scipy.integrate import odeint
Vamos a integrar primero una EDO elemental, cuya solución ya conocemos:
$$y' + y = 0$$$$f(y, t) = \frac{dy}{dt} = -y$$def f(y, t):
return -y
Condiciones iniciales:
y0 = 1
Vector de tiempos donde realizamos la integración:
t = np.linspace(0, 3)
Integramos y representamos la solución:
sol = odeint(f, y0, t)
plt.plot(t, sol)
[<matplotlib.lines.Line2D at 0x7f47eab07128>]
Tendremos que acordarnos ahora de cómo reducir las ecuaciones de orden. De nuevo, vamos a probar con un ejemplo académico:
def f(y, t):
return np.array([y[1], -y[0]])
t = np.linspace(0, 10)
y0 = np.array([1.0, 0.0])
sol = odeint(f, y0, t)
plt.plot(t, sol[:, 0], label='$y$')
plt.plot(t, sol[:, 1], '--k', label='$\dot{y}$')
plt.legend()
<matplotlib.legend.Legend at 0x7f47eaae7828>
En nuestra edición anterior del curso de AeroPython puedes ver una aplicación muy interesante de lo que hemos visto hasta ahora al salto de Felix Baumgartner. ¡Aquí lo tienes!
$$\displaystyle m \frac{d^2 y}{d t^2} = -m g + D$$Clase en vídeo, parte del Curso de Python para científicos e ingenieros grabado en la Escuela Politécnica Superior de la Universidad de Alicante.
from IPython.display import YouTubeVideo
YouTubeVideo("1R_JnTajrRY", width=560, height=315, list="PLGBbVX_WvN7bMwYe7wWV5TZt1a58jTggB")
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
%%html
<a href="https://twitter.com/AeroPython" class="twitter-follow-button" data-show-count="false">Follow @AeroPython</a>
<script>!function(d,s,id){var js,fjs=d.getElementsByTagName(s)[0],p=/^http:/.test(d.location)?'http':'https';if(!d.getElementById(id)){js=d.createElement(s);js.id=id;js.src=p+'://platform.twitter.com/widgets.js';fjs.parentNode.insertBefore(js,fjs);}}(document, 'script', 'twitter-wjs');</script>
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../static/styles/style.css'
HTML(open(css_file, "r").read())