Cómo resolver ecuaciones algebraicas en Python con SciPy

Introducción

En este artículo vamos a utilizar las rutinas de búsqueda de raíces ya disponibles en el módulo scipy.optimize para resolver ecuaciones algebraicas con Python. Ya vimos hace tiempo cómo encontrar el mínimo de una función con SciPy, y también cómo implementar los métodos de la bisección y de Newton en Python. Ahora, además, exploraremos el caso de sistemas de ecuaciones.
En esta entrada se ha usado python 2.7.3, numpy 1.6.2 y scipy 0.11.

Ejemplo básico: métodos de Brent y de Newton

Aunque el método de la bisección es muy conocido porque conceptualmente es muy sencillo de entender, nos vamos a olvidar de él y vamos a utilizar en su lugar el método de Brent, que viene implementado en SciPy en la función scipy.optimize.brentq, ya que según la documentación es la mejor opción. El método de Brent es un algoritmo complicado que combina el método de la bisección, el de la secante e interpolación cuadrática inversa: como el método de la bisección, tiene convergencia garantizada siempre que demos un intervalo \([a, b]\) donde \(f(a) f(b) < 0\); es decir, que haya un cambio de signo en el intervalo.

En el caso en que conozcamos un valor próximo a la solución, seguiremos utilizando el método de Newton o el de la secante a través de la función scipy.optimize.newton, que selecciona uno u otro en el caso de que le pasemos la función derivada o no. Vamos a resolver este ejemplo que me acabo de inventar:

\( \displaystyle e^{x / 3} \cos{x} + 10 \sin{3 x} = x^2 / 4\)

Soluciones de la ecuación

Aplicación: funciones implícitas

Podemos utilizar estos métodos de búsqueda de raíces de ecuaciones algebraicas para representar funciones implícitas. Mientras que las funciones explícitas vienen dadas de la forma \( y = f(x)\), las implícitas están expresadas como \( F(x, y) = 0\), de tal forma que no es posible «despejar» una variable en función de la otra. Sin embargo, podemos obtener una representación siguiendo el siguiente método:

  1. Discretizamos el intervalo \( [x_0, x_{N – 1}]\) en \( N\) elementos.
  2. Para cada valor \( \tilde{x}_i\) con \( i \in 0, …, N – 1\) resolvemos la ecuación \( F(\tilde{x}_i, y) = 0\) y almacenamos la solución en \(y_i\).
  3. Proseguimos hasta que tengamos la lista de valores \( y_0, y_1, …, y_{N – 1}\).

Lo bueno de este algoritmo es que, una vez dado el primer paso, para el siguiente tenemos un valor muy cerca de la solución, por lo que podemos usar métodos de convergencia rápida como el de Newton.

Vamos a tomar un ejemplo con el que me encontré el otro día mientras estudiaba Aerotermodinámica. Tenemos una tobera cuya área transversal viene definida, en función de la coordenada (x) sobre el eje, por
\( A(x) = \begin{cases} 2 x^3 – 3 x^2 + 2 & 0 \leq x \le 1 \\ -\frac{3}{8} x^3 + \frac{9}{4} x^2 – \frac{27}{8} x + \frac{5}{2} & 1 \leq x \le 3 \end{cases}\)
y que tiene esta pinta:


Nota: Recuerda que puedes leer en Pybonacci cómo definir funciones definidas a trozos en NumPy.

Y quiero conocer la distribución del número de Mach \( M\) a lo largo de la misma, utilizando la ecuación
\( \displaystyle \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)}}\)
donde \( A^*\) es el área crítica y \( \gamma = 1.4\) para el caso del aire. Date cuenta de que es imposible despejar \( x\) en función de \( M\) o viceversa, de tal forma que tenemos que resolver la relación implícita que existe entre estas variables. Este es el código:

Otra cosa buena de este método es que, en realidad, este tipo de funciones implícitas definen curvas con más de una «rama» pero, físicamente, para la solución de nuestro problema, solo una suele ser la correcta, como es este caso. Si resolvemos la ecuación de esta manera podemos solucionar la rama que nos interesa.

Sistemas de ecuaciones no lineales

Si estamos en el caso de tener un sistema de ecuaciones no lineales, SciPy también tiene métodos para resolverlos, a través del paquete MINPACK. Utilizaremos la función root, que además soporta varios métodos de resolución en función de si proporcionamos el jacobiano del sistema o no.

Vamos a resolver este ejemplo extraído de la documentación de SciPy:
\( \begin{cases} x \cos{y} & = 4 \\ x y – y & = 5 \end{cases}\)
Una vez más, tenemos que transformar este sistema en \( \vec{f}(\vec{x}) = \vec{0}\), donde las flechas indican vectores. En este caso,
\( \vec{x} = \begin{pmatrix} x \\ y \end{pmatrix}\)
\( \vec{f}(\vec{x}) \equiv \begin{pmatrix} x \cos{y} – 4 \\ x y – y – 5 \end{pmatrix} = \vec{0}\)

El código en Python será:

Y esto ha sido todo, no olvides hacernos llegar tus sugerencias y comentarios.

¡Un saludo!

4 pensamientos sobre “Cómo resolver ecuaciones algebraicas en Python con SciPy”

  1. Viendo la tobera me acabo de acordar de una duda que tengo desde hace tiempo, y seguro que me puedes responder u orientarme hacia las lecturas correctas 😉
    Tengo asumido que una tobera divergente acelera fluidos supersónicos (haces un par de hipótesis, pones las ecuaciones y sale) pero, ¿cuál es aquí la diferencia fundamental (si es que hay alguna relevante e “intuitiva”…) entre el régimen súper y subsónico para entender “intuitivamente” por qué es así? Gracias!
    Como siempre, un gran post.

    1. ¡Hola Jorge! Muchas gracias por el comentario; tengo que reconocer que, aunque me he hecho esa misma pregunta más de una vez y he debido de escuchar la respuesta en clase más de una vez también, he tenido que acudir a Internet a refrescar mi memoria (shame on me). Y he llegado a donde llega uno siempre cuando le surgen dudas sobre estos temas: The Beginner’s Guide to Aeronautics 🙂
      http://www.grc.nasa.gov/WWW/k-12/rocket/nozzle.html
      La relación entre velocidad, área y número de Mach es esta:
      $latex displaystyle (1 – M^2) frac{dV}{V} = – frac{dA}{A}$
      Puedes leer la información ahí, pero lo que sucede cualitativamente es que, cuando el fluido es subsónico, la densidad apenas cambia, mientras que en el caso supersónico los efectos de compresibilidad son notables, y para cumplir la conservación de la masa el fluido se tiene que acelerar.
      Espero haberte aclarado la duda, ¡un saludo!

      1. ¡Ah, la compresibilidad! La página de la nasa parece una buena introducción para muchos temas…
        Muchas gracias por la explicación (y por el correo de aviso, ha sido toda una sorpresa).
        Un saludo!

  2. Pingback: ¿Cómo borrar por encima de una línea en matplotlib? | Pybonacci

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

+ sixty one = sixty eight