¿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.
#Usaremos Arrays:
#Pintaremos cosas:
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)
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 =
B =
v1 =
v2 =
v3 =
print(A,'\n\n', B, '\n\n',v1, '\n\n',v2, '\n\n',v3)
np.dot(A,B)
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 =
M =
x = np.??????(M, V)
x
Cargar y guardar datos
Numpy tiene dos funciones específicas para leer matrices y guardarlas: np.loadtxt y np.savetxt
np.?????('array.txt', x, fmt='%.4e', header='Nuestro array')
z = np.??????('array.txt')
z
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:
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.
Para hacernos una idea de las posibles soluciones siempre podemos representar gráficamente esa función:
Y utilizando por ejemplo el método de Brent en el intervalo $[0, 3]$:
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:
.
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
:
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$!
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?
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?
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.
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:
for ee in 0.0167, 0.249, 0.432, 0.775, 0.967:
pass
# plt.legend(["Earth", "Pluto", "Comet Holmes", "28P/Neujmin", "Halley's Comet"], loc=2)
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}$.
Vamos a integrar primero una EDO elemental, cuya solución ya conocemos:
$$y' + y = 0$$$$f(y, t) = \frac{dy}{dt} = -y$$
Condiciones iniciales:
Vector de tiempos donde realizamos la integración:
Integramos y representamos la solución:
Tendremos que acordarnos ahora de cómo reducir las ecuaciones de orden. De nuevo, vamos a probar con un ejemplo académico:
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())