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\)

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 9, 100)
plt.plot(
    x, np.exp(x / 3) * np.cos(x) + 10 * np.sin(3 * x), 
    x, x ** 2 / 4
)
plt.show()

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:

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)")
plt.show()

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:

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 range(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("Distribución de $M$ a lo largo del eje de la tobera")
plt.ylabel("M")
plt.xlabel("x (m)")
plt.annotate(
    s="Garganta", 
    xy=(1.0, 1.0), 
    xytext=(0.5, 1.6), 
    arrowprops=dict(arrowstyle = "->")
)
plt.show()

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á:

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!