![\"AeroPython\"](\"../static/aeropython_name_mini.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Clase 4b: Ecuaciones no lineales y EDOs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_¿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._"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Ecuaciones no lineales"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy import optimize"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"¿No habíamos dicho que en Python no se puede dividir por cero? Observa esto:
"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"ename": "ZeroDivisionError",
"evalue": "division by zero",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mZeroDivisionError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32mSi manejamos arrays de NumPy las operaciones siguen las reglas dadas en el estándar de punto flotante (IEEE 754). Las divisiones por cero resultan en infinito, 0 / 0 es NaN, etc. Podemos controlar si queremos warnings o errores con la función `np.seterr`.
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ejercicio"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Argumentos extra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vamos a resolver ahora una ecuación que depende de un parámetro:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\\sqrt{x} + \\log{x} = C$$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def G(x, C):\n",
" return C - np.sqrt(x) - np.log(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**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`:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.8773216666875554"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"optimize.newton(G, 2.0, args=(2,))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Flujo compresible"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Esta es la relación isentrópica entre el número de Mach $M(x)$ en un conducto de área $A(x)$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ \\frac{A(x)}{A^*} = \\frac{1}{M(x)} \\left( \\frac{2}{1 + \\gamma} \\left( 1 + \\frac{\\gamma - 1}{2} M(x)^2 \\right) \\right)^{\\frac{\\gamma + 1}{2 (\\gamma - 1)}}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Para un conducto convergente:\n",
"\n",
"$$ \\frac{A(x)}{A^*} = 3 - 2 x \\quad x \\in [0, 1]$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hallar el número de Mach en la sección $x = 0.9$."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"**¡Importante!**: La función del sistema recibe como primer argumento $\\mathbf{y}$ (un array) y como segundo argumento el instante $t$ (un escalar). Esta convención va exactamente al revés que en MATLAB y si se hace al revés obtendremos errores o, lo que es peor, resultados incorrectos.
"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy.integrate import odeint"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vamos a integrar primero una EDO elemental, cuya solución ya conocemos:\n",
"\n",
"$$y' + y = 0$$\n",
"\n",
"$$f(y, t) = \\frac{dy}{dt} = -y$$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def f(y, t):\n",
" return -y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Condiciones iniciales:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"y0 = 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vector de tiempos donde realizamos la integración:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"t = np.linspace(0, 3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Integramos y representamos la solución:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"[![\"Licencia](\"http://i.creativecommons.org/l/by/4.0/88x31.png\")