{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"AeroPython\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Clase 4: 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": [ "#Usaremos Arrays:\n" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#Pintaremos cosas:\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##Ecuaciones algebraicas\n", "\n", "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.\n", "\n", "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:\n", "\n", "* Funciones básicas (norma de un vector, inversa de una matriz, determinante, traza)\n", "* Resolución de sistemas\n", "* Autovalores y autovectores\n", "* Descomposiciones matriciales (QR, SVD)\n", "* Pseudoinversas" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "help(np.linalg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`.\n", "\n", "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.\n", "\n", "Juguemos un poco con ella para ver cómo funciona!\n", "\n", "**Ejercicio:**\n", "* Crear una matriz de 2x3, y otra matriz de 3x2\n", "* Crear dos arrays de 3 elementos y otro de 2\n", "* Multiplicar los arrays entre sí y observar el resultado" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "A = \n", "B = \n", "\n", "v1 = \n", "v2 = \n", "v3 = \n", "\n", "print(A,'\\n\\n', B, '\\n\\n',v1, '\\n\\n',v2, '\\n\\n',v3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "np.dot(A,B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Existe también otra función muy interesante: np.linalg.solve\n", "La usaremos para: \n", "\n", "**Resolver el siguiente sistema:**\n", "\n", "$$ \\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} $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "M =\n", "\n", "M =" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x = np.??????(M, V)\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Cargar y guardar datos**\n", "\n", "Numpy tiene dos funciones específicas para leer matrices y guardarlas: np.loadtxt y np.savetxt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "np.?????('array.txt', x, fmt='%.4e', header='Nuestro array')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "z = np.??????('array.txt')\n", "z" ] }, { "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": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "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": [ "
**Nota**: La función `root` se utiliza para hallar soluciones de *sistemas* de ecuaciones no lineales así que obviamente también funciona para ecuaciones escalares. No obstante, vamos a utilizar las funciones `brentq` y `newton` para que el método utilizado quede más claro.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hay básicamente dos tipos de algoritmos para hallar raíces de ecuaciones no lineales:\n", "\n", "* Aquellos que operan en un intervalo $[a, b]$ tal que $f(a) \\cdot f(b) < 0$. Más lentos, convergencia asegurada.\n", "* Aquellos que operan dando una condición inicial $x_0$ más o menos cerca de la solución. Más rápidos, convergencia condicionada.\n", "\n", "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)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Ejemplo**:\n", "\n", "$\\ln{x} = \\sin{x} \\Rightarrow F(x) \\equiv \\ln{x} - \\sin{x} = 0$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para hacernos una idea de las posibles soluciones siempre podemos representar gráficamente esa función:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Y utilizando por ejemplo el método de Brent en el intervalo $[0, 3]$:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "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": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Si 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": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "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": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "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": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFmdJREFUeJzt3VlsXNd9x/Hff/YhZVmRJVMWSZGURFKbZcuSZTlJEbZN\nC8UPTpMuabokTYvGDeCiDwWapUWtPnTJQ4EgCGAIzgI/xQHaIlUTN6nbRq1RJEaUxI6bWKjVxoDs\npEKzOE1jbRT/fThnOCOapIZzyTvDOd8PMNBcziXn+Jr8nXPPNubuAgCko9DtAgAA8kXwA0BiCH4A\nSAzBDwCJIfgBIDEEPwAkJnPwm9nHzeyCmT27zDkfNrPnzewZMzuU9T0BAJ1bjRb/JyQdX+pFM7tP\n0m53n5T0bkkPr8J7AgA6lDn43f1JST9Y5pT7JT0az31K0iYzG8r6vgCAzuTRxz8s6XzL8YuSRnJ4\nXwDAIvIa3LUFx+wTAQBdUsrhPV6SNNpyPBK/dh0zozIAgA64+8LG9bLyaPGfkvQOSTKzY5JedvcL\ni53or3xDPvuyfG5O7p7s46GHHup6GXrlwbXgWnAtln90InOL38w+KekNkraY2XlJD0kqS5K7n3T3\nx83sPjM7J+nHkt615A87e0zyS5IKUukWqbxdqoxLtUmpMiZVRqTySPi3uFmyFVVyAACtQvC7+9vb\nOOfBtn7Y3I+az69+OzxeOROOrS5ZOf7Ay5Jfi5XDbVJ1XKpOhkqiMiJVRkMFUdpC5QAAC+TRx786\n/GJ4tJq9EB4Xnw7HVpOsEs9vVA6bF9w5jIeKoVFBFG/pucphZmam20XoGVyLJq5FE9ciG+u0j2i1\nmZn7mTzeaLHKIXYrVSfCnUN1PN41jMbKYVPPVQ4AIElmJl/h4G56wd+O+W4lD5WDFLqNysOxcpiW\nqmMLKocNXS0ygDQR/LmxlsphTpq7KBUqUunWUAlUd8Uxh1g5VHZIleHmGAUArBKCv6cUpUI9/KvZ\nUDkUN0qlbeFuoToVKojKjvgYlUpb6VICsCIE/3pj1fCQWsYbtkiV7fGuYTqON+xo3j0U6l0tMoDe\nQvD3HZMKA5KVQqUwdzEEf3lbqAiqU1Jtd6wU4p1D6VbuGoCEEPwpskqYqSRJfkXy2XjXMBLuGmp7\nwoB0466hPBLGIwD0BYIfi4h3DSrpurGG8napMiHVpmO30li8cxiTijd1u9AA2kTwozPzaxu8ZYbS\nUKgEanuk2lSzUqiOs10G0EMIfqyN67qTLkvyOHV1LFQKjUHo6nj4WmmIigHICcGPLinF2UbWHGco\nN9Y0TLWMM4yHyqE0JFleHwUB9DeCHz2qGMcZCnHa6qxU3hoGnKvTCyqGCWYmAStA8GOdalQM8Y5B\n18Jitsp47EqKFUPj0YMb6wHdQvCjT7V2JTX2ThoK3Ua1vS13DI2KgVlJSAfBjzTNDz43ZiVV426r\nO6XagbAdd3VCquwM3UusY0AfIfiBxVg9TledleYuhW22KyNhDUNtn1TdHSqJ6gQzkrDuEPzAihXi\n+EIxdiNdi91IE1Jtv1TfG+4UGhUDeyWhxxD8wGqb30jPpblXwuculEfDoHP9QPzgnlgxcLeALiD4\ngVyZVBiUrCjNXZY0F7fd3hkqhdqeuB3GrrCwjbEFrAGCH+glVot3C3Fn1dJr4mrnvXHQeXezYmAm\nEjpE8APrRqk5XtDYbrsyGgac6wdDF1Jtdxh4Zm8kLIPgB/pCoWWl8yVJxbjN9lSoFOrTcSbSblY5\ng+AH+p+1zEK6Isml8nBYq1A/GLqRqrvD3UJpG5VCAgh+IHWFQUml5i6q5e2h22jgjji2MBmOuVPo\nGwQ/gCXYgvUKirunTi+oFKbCIDTWDYIfQAfitNTGmIKVm7OP6nfGMYWpONA82O3CYgGCH8Aqa+yc\nqvixnTfFdQq3xzGFRqUwLlmpqyVNFcEPID/zm+NdCwvYykNxPOGQVN8fK4VpqbSF8YQ1RPAD6A3z\ng8yXwsrmylhYzVw/FD/HOU5JZTVzZgQ/gB5Xil1Hc3E189Y4FfUuaeD2uGPqnnCXgLYQ/ADWr/m7\nhIthgHl+LOFQ2CW1tjfcOVix2yXtKQQ/gD5UkQpxLMGvxAVr+6TBI82ts6uT8Zz0EPwAElKIdwmK\n3Ua3hK6igbtjt1G8S+jzDfAIfgCQpMIGSYX4GQo3hc9PGDgcZhzV9oUKoU8WqhH8ALAcGwjrDfxi\nGGSuTkoDd4VHfX/oOlpnFQLBDwCdsHoYUL6uQjgSK4QDUn2fVLy526VcFMEPAKvpujuEm1rGEA41\nK4TGyuZuFZHgB4AcFAYlFcMYQumW0E00cE/Y8K5+e1ibYOVcikLwA0DXWBxUdmnuUvjwnPqd0uAx\naeBgqBDKw6u+fQXBDwA9pxw+WtOvhtCvTkmDR6XBu+NGd/sz7XpK8APAetEYP5jvLrpd2vC6OH5w\nUKqMt3V3QPADwLpWjOMHs5LPhbGCgbtjd9GdYUC5UL/uOwh+AOhHrXcH5W2hEhh8vbThXtnGNxD8\nAJCGsqRZ2RFfcfAX1qhEAIA1dbXj7yT4ASAxmYPfzI6b2Vkze97M3rvI6zNm9kMz+1p8/FHW9wQA\ndC7TpyObWVHSRyS9UdJLkr5sZqfc/bkFp/6Lu9+f5b0AAKsja4v/qKRz7v6Cu1+V9JikNy9yHp+0\nDAA9ImvwD0s633L8YvxaK5f0WjN7xsweN7N9Gd8TAJBBpq4ehVC/ka9KGnX3V8zsTZI+LWlqsRNP\nnGw+nzkszRzJWDoA6DOnz0inv9I46mw6fqZ5/GZ2TNIJdz8ej98vac7dP7jM93xL0mF3//6CrzOP\nHwBWxLoyj/+MpEkzGzeziqS3STp1XbHMhszChhNmdlShsvn+q38UACAPmbp63H3WzB6U9HlJRUkf\nc/fnzOyB+PpJSb8g6T1mNivpFUm/nLHMAIAM2LIBANat7nT1AADWGYIfABJD8ANAYgh+AEgMwQ8A\niSH4ASAxBD8AJIbgB4DEEPwAkBiCHwASQ/ADQGIIfgBIDMEPAIkh+AEgMQQ/ACSG4AeAxBD8AJAY\ngh8AEkPwA0BiCH4ASAzBDwCJIfgBIDEEPwAkhuAHgMQQ/ACQGIIfABJD8ANAYgh+AEgMwQ8AiSH4\nASAxBD8AJIbgB4DEEPwAkBiCHwASQ/ADQGIIfgBIDMEPAIkh+AEgMQQ/ACSG4AeAxBD8AJAYgh8A\nEkPwA0BiCH4ASAzBDwCJIfgBIDEEPwAkJnPwm9lxMztrZs+b2XuXOOfD8fVnzOxQ1vcEAHQuU/Cb\nWVHSRyQdl7RP0tvNbO+Cc+6TtNvdJyW9W9LDWd4TAJBN1hb/UUnn3P0Fd78q6TFJb15wzv2SHpUk\nd39K0iYzG8r4vgCADmUN/mFJ51uOX4xfu9E5IxnfFwDQoVLG7/c2z7N2vu/EyebzmcPSzJHOCgUA\n/er0Gen0VxpH7Ubw9bIG/0uSRluORxVa9MudMxK/9ionHshYGgDoczNHWhvFpj95ZOXhn7Wr54yk\nSTMbN7OKpLdJOrXgnFOS3iFJZnZM0svufiHj+wIAOpSpxe/us2b2oKTPSypK+pi7P2dmD8TXT7r7\n42Z2n5mdk/RjSe/KXGoAQMfMvbM+otVmZu5nul0KAFhPTHbE5e4Lx1GXxcpdAEgMwQ8AiSH4ASAx\nBD8AJIbgB4DEEPwAkBiCHwASQ/ADQGIIfgBIDMEPAIkh+AEgMQQ/ACSG4AeAxBD8AJAYgh8AEkPw\nA0BiCH4ASAzBDwCJIfgBIDEEPwAkhuAHgMQQ/ACQGIIfABJD8ANAYgh+AEgMwQ8AiSH4ASAxBD8A\nJIbgB4DEEPwAkBiCHwASQ/ADQGIIfgBIDMEPAIkh+AEgMQQ/ACSG4AeAxBD8AJAYgh8AEkPwA0Bi\nCH4ASAzBDwCJIfgBIDEEPwAkhuAHgMQQ/ACQGIIfABJT6vQbzWyzpE9JGpP0gqRfcveXFznvBUn/\nK+mapKvufrTT9wQAZJelxf8+SU+4+5Skf4rHi3FJM+5+iNAHgO7LEvz3S3o0Pn9U0s8tc65leB8A\nwCrKEvxD7n4hPr8gaWiJ81zSP5rZGTP77QzvBwBYBcv28ZvZE5K2LfLSH7YeuLubmS/xY17n7t8x\ns62SnjCzs+7+5GInnjjZfD5zWJo5slzpACA9p89Ip7/SOFoqdpdn7h1+o9lZhb77/zaz2yR9wd33\n3OB7HpL0f+7+l4u85n6mo6IAQKJMdsTl7ivqTs/S1XNK0jvj83dK+vSrimQ2YGY3xeeDkn5W0rMZ\n3hMAkFGW4P8LST9jZv8h6afiscxsu5l9Np6zTdKTZva0pKckfcbd/yFLgQEA2XTc1bPa6OoBgJXK\nv6sHALAOEfwAkBiCHwDWC6tJhY2SlaXSkLTprR39mI736gEArBWTChskueSXpcqENHhEGrxXGjgk\n1W+Xihub564QwQ8A3WTV8PBLUmFQqu2TNtwrDRyR6gel2qRkqxvVBD8A5KIQgl1zkl8JrfiBw9Lg\nPdLAwRDypVtyKQnBDwCrrTAgqSTNvRLCvH5AGjwmDdwZAr66S7Ji14pH8ANAp6wmFarS3MXwvDYV\numgGG900+6TiYLdL+SoEPwDcUDm04v2KJJeqO0M3zcDd0sDtUm2/VN7a7UK2jeAHgHml2E0zK/ms\nVBkL3TMDR0N3Tf2AVB6WbH1/xAjBDyBBjYC/JvlVqbIjdM0MHJXq+8OjMiZZfy51IvgB9C+rSFYP\n4a5rUmU8zKAZOBLCvbY3fK1PA34pBD+A9c/qIeT9cgjxyoQ0cIdUPyzV94VHeXTdd9GsFoIfwPpR\nGJRUDNMkixuk6m6pfmdYzVrbE1rxpVsJ+Bsg+AH0mGJzodPcRak8JFWnpcHDYfZMfW8I+eLN3S7o\nukXwA+gOq0uFijR3ORxXxkKLfeCuMP+9vleq7AznYFUR/ADWUDHOnlFcxbo5ds/cEfrga3vCozRE\n90yOCH4A2TX63v2ypEJsve+Lrfc9Um06BH6h1u2SQgQ/gHZZNWxLoNnQPVMekqqTcf+ZAzHcp6TS\nVlrvPY7gB9CisbDJ48yZTWF7gvqBEPDVqbAfTWXHqm8VjPzwfw5ITqPf3UK4Fwak6lgcUL0ztNxr\nU2EHSbpm+hLBD/Slxt7vJvnF0E1T2RFCfT7cJ0O/+/wnOSEVBD+wbi0W7qNhzvvAHbHPPYZ7aVO3\nC4seQvADPa2xmEmxW6Yewr02HaZENmbLEO5YAYIf6DarxNkyjQHVjbFbZk8M99hqr+6Sijd1u7To\nAwQ/kIfGKlWfleYuSaUtUnUi7A5Z298S7hOhVQ+sIYIfWBWFOFOmEBcxSSrfFlep7g+t9+qucFwZ\nZSokuorfPqBdVg2P+S6Zm+Jg6lSY597ojqnuYhETehrBD8yLrXYrxo3DrkmlbbFLprFp2EQMd7pk\nsH4R/EhL4wM7dC1s+Vu8WaqMxHnt+5v97NWdbByGvkXwo8+Umy1xvyipJJW3x20H9sftBiZCuFfG\npUK1m4UFuoLgxzrTst2AX1HojhkKId7Y4re6U6qOh4BnbjvwKgQ/ekxjdkwxBLtfCVMfG4uWavtC\nsFfGQ6udj9kDVozgR84srES1ouRXwyBqabNUHglz2RvB3uiKKd8WzgWwagh+rLLWFvvVMKf9umCP\n89kr46E7prydOe1AzviLwwoV4+Bpoyvmagj2ymiYEVPbG1vrEy3BTosd6CUEPxaoNPdg98uS5kI/\nemVH2OmxNh0HTsfCo7xdskI3CwxghQj+1FhdsrLC6tOLYf+Y0lAI8cYHcDRCvTIWBlYZPAX6CsHf\nV1pWnvpsc4FS+bbQ9VLbI9V2tQT7DnZ7BBJE8K8n83vFqKUbZqtUHpVqu8MHcFTHQ6BXxqTKcGzd\nA0ATwd8zFrbWL4XWeGlbCPPqZAj3+VDfIRU30w0DYMUI/rzM7xEzJ/klSYXQWq+Mxk2/psMHXld2\nhK+Vh9lOAMCaIPhXg1ViF4zFuetX4tz12+JWAtNxiuOO8CiPhr53WusAuoDgv6FSnN5YlHQ1dsFs\nbJkJMxkXJI02W+ulIaY4AuhZiQd/y2IkxX71Qj2G+mjzQ6yrO8LK08ooA6YA1r0+Dv6WrQMae68X\n6rFffST2q0/GVvpIS796rdsFB4A11XHwm9kvSjohaY+ku939q0ucd1zShxQS+KPu/sFO37OptaW+\nRKhXdoWWemU0ttZH+MQkAFC2Fv+zkt4i6eRSJ5hZUdJHJL1R0kuSvmxmp9z9uWWLVKhLKqjZ/TIo\nlW8NAV7d2fzA6spwGCitDPdVqJ8+fVozMzPdLkZP4Fo0cS2auBbZdBz87n5Wkmz5mSlHJZ1z9xfi\nuY9JerOkxYN/y+/EvdZHQwu9PBJmxiQ2rZFf6iauRRPXoolrkc1a9/EPSzrfcvyipHuWPHvs4TUu\nDgBg2eA3syckbVvkpQ+4+9+18fO9o1IBANaMuWfLZjP7gqTfX2xw18yOSTrh7sfj8fslzS02wGtm\nVBIA0AF3X9Fq0NXq6lnqTc9ImjSzcUnflvQ2SW9f7MSVFhwA0JmOl5ea2VvM7LykY5I+a2Z/H7++\n3cw+K0nuPivpQUmfl/RNSZ9afkYPAGCtZe7qAQCsL7luKGNmx83srJk9b2bvXeKcD8fXnzGzQ3mW\nL083uhZm9qvxGnzdzP7NzA52o5x5aOf3Ip53t5nNmtlb8yxfntr8G5kxs6+Z2b+b2emci5ibNv5G\ntpjZ58zs6XgtfqMLxVxzZvZxM7tgZs8uc87KctPdc3koLLM9J2lcUlnS05L2LjjnPkmPx+f3SPpS\nXuXL89HmtbhX0s3x+fGUr0XLef8s6TOSfr7b5e7i78UmSd+QNBKPt3S73F28Fick/XnjOkj6nqRS\nt8u+BtfiJyQdkvTsEq+vODfzbPHPL+Zy96uSGou5Wt0v6VFJcvenJG0ys6Ecy5iXG14Ld/+iu/8w\nHj4laSTnMualnd8LSfpdSX8l6X/yLFzO2rkWvyLpr939RUly9+/mXMa8tHMtviNpY3y+UdL3PIwr\n9hV3f1LSD5Y5ZcW5mWfwL7aYa7iNc/ox8Nq5Fq1+S9Lja1qi7rnhtTCzYYU/+sYKv34dmGrn92JS\n0mYz+4KZnTGzX8+tdPlq51o8Imm/mX1b0jOSfi+nsvWaFedmnrtztvvHunBaZz/+kbf932RmPynp\nNyW9bu2K01XtXIsPSXqfu7uFPUL6depvO9eiLOkuST8taUDSF83sS+7+/JqWLH/tXIsPSHra3WfM\nbJekJ8zsDnf/0RqXrRetKDfzDP6XJI22HI8q1EzLnTMSv9Zv2rkWigO6j0g67u7L3eqtZ+1ci8OS\nHov7Qm2R9CYzu+rup/IpYm7auRbnJX3X3S9Kumhm/yrpDkn9FvztXIvXSvpTSXL3/zSzb0maVlg/\nlJIV52aeXT3zi7nMrKKwmGvhH+4pSe+Q5lf9vuzuF3IsY15ueC3MbIekv5H0a+5+rgtlzMsNr4W7\n73T3CXefUOjnf08fhr7U3t/I30p6vZkVzWxAYTDvmzmXMw/tXIuzCjv/KvZpT0v6r1xL2RtWnJu5\ntfjdfdbMGou5ipI+5u7PmdkD8fWT7v64md1nZuck/VjSu/IqX57auRaS/ljSayQ9HFu6V939aLfK\nvFbavBZJaPNv5KyZfU7S1yXNSXrE3fsu+Nv8vfgzSZ8ws2cUGrF/4O7f71qh14iZfVLSGyRtiYtm\nH1Lo8us4N1nABQCJ4RPBASAxBD8AJIbgB4DEEPwAkBiCHwASQ/ADQGIIfgBIDMEPAIn5f/qGAQ7x\nBpCnAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def A(x):\n", " return 3 - 2 * x\n", "\n", "x = np.linspace(0, 1)\n", "area = A(x)\n", "r = np.sqrt(area / np.pi)\n", "plt.fill_between(x, r, -r, color=\"#ffcc00\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "¿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 :)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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$!" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Ecuación de Kepler" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Representar la ecuación de Kepler\n", "\n", "$$M = E - e \\sin E$$\n", "\n", "que relaciona dos parámetros geométricos de las órbitas elípticas, la anomalía media $M$ y la anomalía excéntrica $E$.\n", "\n", "![Anomalías excéntrica y media](http://upload.wikimedia.org/wikipedia/commons/thumb/f/f8/Kepler%27s-equation-scheme.png/250px-Kepler%27s-equation-scheme.png)\n", "\n", "para los siguientes valores de excentricidad:\n", "\n", "* Tierra: $0.0167$\n", "* Plutón: $0.249$\n", "* Cometa Holmes: $0.432$\n", "* 28P/Neujmin: $0.775$\n", "* Cometa Halley: $0.967$\n", "\n", "Para reproducir esta gráfica:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import HTML\n", "HTML('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para ello utilizaremos el método de Newton (secante).\n", "\n", "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?" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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?" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vamos a introducir aquí un truco muy útil en Python:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for ee in 0.0167, 0.249, 0.432, 0.775, 0.967:\n", " pass\n", "\n", "# plt.legend([\"Earth\", \"Pluto\", \"Comet Holmes\", \"28P/Neujmin\", \"Halley's Comet\"], loc=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ecuaciones diferenciales ordinarias" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para integrar EDOs vamos a usar la función `odeint` del paquete `integrate`, que permite integrar sistemas del tipo:\n", "\n", "$$ \\frac{d\\mathbf{y}}{dt}=\\mathbf{f}\\left(\\mathbf{y},t\\right)$$\n", "\n", "con condiciones iniciales $\\mathbf{y}(\\mathbf{0}) = \\mathbf{y_0}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
**¡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": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "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": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Condiciones iniciales:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vector de tiempos donde realizamos la integración:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Integramos y representamos la solución:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### EDOs de orden superior" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tendremos que acordarnos ahora de cómo reducir las ecuaciones de orden. De nuevo, vamos a probar con un ejemplo académico:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$y + y'' = 0$$\n", "\n", "$$\\mathbf{y} \\leftarrow \\pmatrix{y \\\\ y'}$$\n", "\n", "$$\\mathbf{f}(\\mathbf{y}) = \\frac{d\\mathbf{y}}{dt} = \\pmatrix{y \\\\ y'}' = \\pmatrix{y' \\\\ y''} = \\pmatrix{y' \\\\ -y}$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "Clase en vídeo, parte del [Curso de Python para científicos e ingenieros](http://cacheme.org/curso-online-python-cientifico-ingenieros/) grabado en la Escuela Politécnica Superior de la Universidad de Alicante." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import YouTubeVideo\n", "\n", "YouTubeVideo(\"1R_JnTajrRY\", width=560, height=315, list=\"PLGBbVX_WvN7bMwYe7wWV5TZt1a58jTggB\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si te ha gustado esta clase:\n", "\n", "Tweet\n", "\n", "\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "####

¡Síguenos en Twitter!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###### Follow @AeroPython " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### \"Licencia
Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "_Las siguientes celdas contienen configuración del Notebook_\n", "\n", "_Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como [seguro](http://ipython.org/ipython-doc/dev/notebook/security.html)_\n", "\n", " File > Trusted Notebook" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "Follow @AeroPython\n", "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "Follow @AeroPython\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "El estilo se ha aplicado =)\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Esta celda da el estilo al notebook\n", "from IPython.core.display import HTML\n", "css_file = '../static/styles/style.css'\n", "HTML(open(css_file, \"r\").read())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }