![\"AeroPython\"](\"../static/aeropython_name_mini.png\")
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Clase 4: Ecuaciones algebraicas, 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": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##Ecuaciones algebraicas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Help on package numpy.linalg in numpy:\n",
"\n",
"NAME\n",
" numpy.linalg\n",
"\n",
"DESCRIPTION\n",
" Core Linear Algebra Tools\n",
" -------------------------\n",
" Linear algebra basics:\n",
" \n",
" - norm Vector or matrix norm\n",
" - inv Inverse of a square matrix\n",
" - solve Solve a linear system of equations\n",
" - det Determinant of a square matrix\n",
" - lstsq Solve linear least-squares problem\n",
" - pinv Pseudo-inverse (Moore-Penrose) calculated using a singular\n",
" value decomposition\n",
" - matrix_power Integer power of a square matrix\n",
" \n",
" Eigenvalues and decompositions:\n",
" \n",
" - eig Eigenvalues and vectors of a square matrix\n",
" - eigh Eigenvalues and eigenvectors of a Hermitian matrix\n",
" - eigvals Eigenvalues of a square matrix\n",
" - eigvalsh Eigenvalues of a Hermitian matrix\n",
" - qr QR decomposition of a matrix\n",
" - svd Singular value decomposition of a matrix\n",
" - cholesky Cholesky decomposition of a matrix\n",
" \n",
" Tensor operations:\n",
" \n",
" - tensorsolve Solve a linear tensor equation\n",
" - tensorinv Calculate an inverse of a tensor\n",
" \n",
" Exceptions:\n",
" \n",
" - LinAlgError Indicates a failed linear algebra operation\n",
"\n",
"PACKAGE CONTENTS\n",
" _umath_linalg\n",
" info\n",
" lapack_lite\n",
" linalg\n",
" setup\n",
"\n",
"DATA\n",
" absolute_import = _Feature((2, 5, 0, 'alpha', 1), (3, 0, 0, 'alpha', 0...\n",
" division = _Feature((2, 2, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0), 8192...\n",
" print_function = _Feature((2, 6, 0, 'alpha', 2), (3, 0, 0, 'alpha', 0)...\n",
"\n",
"FILE\n",
" /home/siro/anaconda3/lib/python3.4/site-packages/numpy/linalg/__init__.py\n",
"\n",
"\n"
]
}
],
"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!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**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": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1 2]\n",
" [2 1]\n",
" [1 5]] \n",
"\n",
" [[1 2 3]\n",
" [3 2 1]] \n",
"\n",
" [1 0 1] \n",
"\n",
" [2 2 2] \n",
"\n",
" [1 2]\n"
]
}
],
"source": [
"A = np.array([\n",
" [1, 2],\n",
" [2, 1],\n",
" [1, 5]\n",
" ])\n",
"B = np.array([\n",
" [1, 2, 3],\n",
" [3, 2, 1]\n",
" ])\n",
"\n",
"v1 = np.array([1, 0, 1])\n",
"v2 = np.array([2, 2, 2])\n",
"v3 = np.array([1, 2])\n",
"\n",
"print(A,'\\n\\n', B, '\\n\\n',v1, '\\n\\n',v2, '\\n\\n',v3)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 7 6 5]\n",
" [ 5 6 7]\n",
" [16 12 8]]\n"
]
}
],
"source": [
"print(np.dot(A,B)) #ETC"
]
},
{
"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": 21,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"M = np.array([\n",
" [2, 0, 1],\n",
" [-1, 1, 0],\n",
" [3, 2, -1]\n",
"])\n",
"M = np.array([-1, 3, 0])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([-1., 2., 1.])"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = np.linalg.solve(M, V)\n",
"x"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([-1., 3., 0.])"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.dot(M, 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": 24,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"np.savetxt('array.txt', x, fmt='%.4e', header='Nuestro array')"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"array([-1., 2., 1.])"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z = np.loadtxt('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": 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\")