¿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
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$$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/Pybonacci" class="twitter-follow-button" data-show-count="false">Follow @Pybonacci</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())