Cómo resolver ecuaciones algebraicas en Python con SciPy

Por Juan Luis Cano

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$

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 9, 100)
plt.plot(x, exp(x / 3) * cos(x) + 10 * sin(3 * x), x, x ** 2 / 4)

Ecuación 1

Antes que nada debemos definir la función que va a representar la ecuación. Todas las subrutinas de búsqueda de raíces en realidad buscan ceros de funciones, así que debemos escribir

$f(x) \equiv \displaystyle e^{x / 3} \cos{x} + 10 \sin{3 x} - x^2 / 4 = 0$

def f(x):
    """Ecuación no lineal.
    Asume que se ha importado NumPy de la forma
      import numpy as np
    """
    return np.exp(x / 3.0) * np.cos(x) + 10 * np.sin(3 * x) - x ** 2 / 4

Vamos a buscar primero la solución que vemos que está entre 4 y 5. Para ello, utilizamos el método de Brent:

from scipy.optimize import brentq
sol1 = brentq(f, 4.0, 5.0)
print sol1  # 4.40989799172

Sencillo, ¿no? Ahora si queremos buscar la solución que está cerca de 2, podemos probar el método de la secante:

from scipy.optimize import newton
sol2 = newton(f, 2.0)
print sol2  # 2.17349784856

Más sencillo todavía 😉 estas son las dos soluciones:

Soluciones de la ecuación 1

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:

:::python x = np.linspace(0, 3, 151) def A(x): """Área transversal de la tobera. """ def A1(x): return 2.0 * x ** 3 - 3.0 * x ** 2 + 2.0 def A2(x): return -3.0 * x ** 3 / 8.0 + 9.0 * x ** 2 / 4.0 - 27.0 * x / 8.0 + 5.0 / 2.0 return np.piecewise(x, [(0.0 <= x) & (x < 1.0), (1.0 <= x) & (x <= 3.0)], [A1, A2]) plt.fill_between(x, np.sqrt(A(x) / np.pi), -np.sqrt(A(x) / np.pi), facecolor="#eebb22") plt.xlim((0, 3)) plt.title("Tobera") plt.xlabel("x (m)") plt.ylabel("Radio (dm)")

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:

:::python from scipy.optimize import brentq, newton def rel(M, gamma=1.4): """Parte derecha de la relación entre el número de Mach $M$ y la relación de áreas $A / A^*$. """ return (2 * (1 + (gamma - 1) * M ** 2 / 2) / (gamma + 1)) ** ((gamma + 1) / 2 / (gamma - 1)) / M def eq(M, A): """Función implícita entre el número de Mach y la relación de áreas. """ return rel(M) - A # Para cada valor de x resolvemos la ecuación en M M = np.empty_like(x) # El primer paso lo damos con el método de Brent M[0] = brentq(eq, 0.001, 1, args=(A(x[0:1]),)) # Comenzamos a iterar for i in xrange(1, len(x)): # El valor inicial para el método de Newton es la solución del # paso anterior M[i] = newton(eq, M[i - 1], args=(A(x[i - 1:i]),)) # Representamos la solución plt.plot(x, M) plt.plot(x, np.ones_like(x), 'k--') plt.title(u"Distribución de $M$ a lo largo del eje de la tobera") plt.ylabel("M") plt.xlabel("x (m)") plt.annotate(s=u"Garganta", xy=(1.0, 1.0), xytext=(0.5, 1.6), arrowprops=dict(arrowstyle = "->"))

Otra cosa buena de este método es que, en relalidad, 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á:

:::python from scipy.optimize import root def f(x): """Sistema de dos ecuaciones con dos incógnitas. """ return [ x[0] * np.cos(x[1]) - 4, x[1]*x[0] - x[1] - 5 ] sol = root(f, [1, 1], jac=False) # Devuelve objecto Result print sol.x # Result.x contiene la solución

Y esto ha sido todo, no olvides hacernos llegar tus sugerencias y comentarios. ¡Un saludo!

Comentarios