Resolver ecuaciones en Python: forma cerrada versus numérica

Encuentre soluciones de forma cerrada cuando sea posible; utilice métodos numéricos cuando sea necesario

¿Por qué podemos resolver algunas ecuaciones fácilmente, mientras que otras parecen imposibles? Y otra cosa: ¿por qué se nos oculta este conocimiento?

Como científicos de datos, científicos aplicados e ingenieros, a menudo creamos modelos matemáticos. Por ejemplo, considere el modelo: y = x². Dado un valor para x,wPodemos aplicarlo hacia adelante para calcular y. Por ejemplo, si x = 3entonces y = 9.

También podemos aplicar este modelo al revés. comenzando con y = x²reorganizamos para resolver incógnita: x = ±√y. Si y = 9entonces x = ±3. la expresión x = ±√y es un ejemplo de un solución de forma cerrada — una expresión que utiliza una combinación finita de operaciones y funciones estándar.

Sin embargo, no todos los modelos son tan sencillos. A veces nos encontramos con ecuaciones en las que no podemos simplemente “resolver” incógnita”y obtener una expresión de forma cerrada. En tales casos, podríamos escuchar: “Eso no tiene solución; se necesitan métodos numéricos”. Los métodos numéricos son poderosos. Pueden proporcionar aproximaciones precisas. Aún así, me frustra a mí (y quizás a usted) que nadie parezca explicar cuándo las soluciones de forma cerrada son posibles y cuándo no.

El gran Johannes Kepler compartió nuestra frustración. Al estudiar el movimiento planetario, creó este modelo:

Esta ecuación convierte la posición de un cuerpo a lo largo de su órbita (incógnita) en su tiempo a lo largo de la órbita (y). Kepler buscó una solución de forma cerrada para incógnita convertir el tiempo en una posición. Sin embargo, incluso 400 años después, lo mejor que tenemos son los métodos numéricos.

En este artículo, desarrollaremos una intuición sobre cuándo esperar una solución de formato cerrado. La única manera de determinar esto rigurosamente es mediante el uso de matemáticas avanzadas, como la teoría de Galois, la teoría de números trascendentales y la geometría algebraica. Estos temas van mucho más allá de lo que nosotros, como científicos e ingenieros aplicados, normalmente aprendemos en nuestra formación.

En lugar de sumergirnos en estos campos avanzados, haremos trampa. Usando SymPyun sistema de álgebra computacional basado en Python, exploraremos diferentes clases de ecuaciones para ver cuáles puede resolver con una expresión de forma cerrada. Para completar, también aplicaremos métodos numéricos.

Exploraremos ecuaciones que combinan polinomios, exponenciales, logaritmos y funciones trigonométricas. En el camino, descubriremos combinaciones específicas que a menudo se resisten a las soluciones de forma cerrada. Veremos que si quieres crear una ecuación. con (o sin) una solución de forma cerrada, deberías evitar (o intentar) la siguiente:

• Polinomios de quinto grado y superiores
• Mezclando incógnita con exp(incógnita) o iniciar sesión (incógnita) – si Lambert W. la función está fuera de los límites
• Combinando exp(incógnita) y iniciar sesión (incógnita) dentro de la misma ecuación
• Algunos pares de funciones trigonométricas con frecuencias proporcionales
• Muchos pares de funciones trigonométricas con frecuencias no proporcionales.
• Mezclar funciones trigonométricas con incógnitaexp(x) o log(x)

Aparte 1: no soy matemático y mis scripts SymPy no son matemáticas superiores. Si encuentra algún error o recurso pasado por alto, perdone mi descuido. Por favor compártelos conmigo y con gusto agregaré una nota.

Aparte 2: vídeo reciente de Welch Lab, La ecuación imposible de Keplerme recordó mi frustración por no saber cuándo se puede resolver una ecuación de forma cerrada. El vídeo desató la investigación que sigue y proporciona nuestro primer ejemplo.

Imagine que es el programador de investigación de Johannes Kepler. Ha creado el siguiente modelo de movimiento orbital:

y = incógnita −do pecado(incógnita)

dónde:

• incógnita es la posición del cuerpo a lo largo de su órbita. Medimos esta posición como un ángulo (en radianes). El ángulo comienza en 0 radianes cuando el cuerpo está más cerca del Sol. Cuando el cuerpo ha recorrido ¼ de la distancia de su órbita, el ángulo es π/2 radianes (90°). Cuando ha recorrido la mitad de la distancia de su órbita, el ángulo es π (180°), y así sucesivamente. Recuerde que los radianes miden ángulos de 0 a 2π en lugar de 0 a 360°.
• do es la excentricidad de la órbita, que va desde 0 (un círculo perfecto) hasta poco menos de 1 (una elipse muy alargada). Supongamos que Kepler ha observado un cometa con c = 0,967.
• y es el tiempo que pasa el cuerpo a lo largo de su órbita. Medimos este tiempo como un ángulo (en radianes). Por ejemplo, si el cometa tiene un período orbital de 76 años terrestres, entonces π/2 (90°) corresponde a ¼ de 76 años, o 19 años. Un tiempo de π (180°) corresponde a la mitad de 76 años, o 38 años. Un tiempo de 2π (360°) es el período orbital completo de 76 años.

Este diagrama muestra la posición del cometa en π/2 radianes (90°), que es ¼ del camino a lo largo de su órbita:

Kepler pregunta la hora en que el cometa alcanza la posición π/2 radianes (90°). Creas y ejecutas este código Python:

import numpy as np
def kepler_equation(x):
return x - c * np.sin(x)
c = 0.967
position_radians = np.pi / 2 # aka 90 degrees
time_radians = kepler_equation(position_radians)
orbital_period_earth_years = 76
t_earth_years = (time_radians / (2 * np.pi)) * orbital_period_earth_years
print(f"It takes approximately {t_earth_years:.2f} Earth years for the comet to move from 0 to π/2 radians.")

Le informa a Kepler:

It takes approximately 7.30 Earth years for the comet to move from 0 to π/2 radians.

Aparte: el cometa cubre el 25% de su distancia orbital en menos del 10% de su período orbital porque acelera cuando está más cerca del Sol.

Ninguna buena acción queda impune. Kepler, fascinado por el resultado, te asigna una nueva tarea: “¿Puedes decirme qué tan lejos se encuentra el cometa en su órbita después de 20 años terrestres? Quiero saber la posición en radianes”.

“No hay problema”, piensas. “Solo usaré un poco de álgebra de la escuela secundaria”.

Primero, conviertes 20 años terrestres a radianes:

• tiempo_radianes = (20/76) × 2π = (10/19)π

A continuación, reordenas la ecuación de Kepler, igualándola a 0.

• incógnita − 0,967 pecado(incógnita) − (10/19)π = 0

Ahora quieres encontrar el valor de incógnita eso hace que esta ecuación sea verdadera. Decides graficar la ecuación para ver dónde cruza cero:

import numpy as np
import matplotlib.pyplot as plt
c = 0.967
time_earth_years = 20
orbital_period_earth_years = 76
time_radians = (time_earth_years / orbital_period_earth_years) * 2 * np.pi
def function_to_plot(x):
return x - c * np.sin(x) - time_radians
x_vals = np.linspace(0, 2 * np.pi, 1000)
function_values = function_to_plot(x_vals)
plt.figure(figsize=(10, 6))
plt.axhline(0, color='black', linestyle='--') # dashed horizontal line at y=0
plt.xlabel("Position (radians)")
plt.ylabel("Function Value")
plt.title("Graph of x - c sin(x) - y to Find the Root")
plt.grid(True)
plt.plot(x_vals, function_values)
plt.show()

Hasta ahora, todo bien. El gráfico muestra que una solución para incógnita existe. Pero cuando intentas reorganizar la ecuación para resolver incógnita Usando álgebra, chocas contra una pared. ¿Cómo se aísla? incógnita cuando tienes una combinación de incógnita y el pecado (incógnita)?

“Está bien”, piensas. “Tenemos Python, y Python tiene la Paquete SymPy”, un sistema de álgebra informático potente y gratuito.

Le planteas el problema a SymPy:

# Warning: This code will fail.
import sympy as sym
from sympy import pi, sin
from sympy.abc import x
c = 0.967
time_earth_years = 20
orbital_period_earth_years = 76
time_radians = (time_earth_years / orbital_period_earth_years) * 2 * pi
equation = x - c * sin(x) - time_radians
solution = sym.solve(equation, x)
#^^^^^^^^^^^^^error^^^^^^^^^^^^^^
print(solution)

Desafortunadamente, responde con un error:

NotImplementedError: multiple generators [x, sin(x)]
No algorithms are implemented to solve equation x - 967*sin(x)/1000 - 10*pi/19

SymPy es bastante bueno resolviendo ecuaciones, pero no todas las ecuaciones se pueden resolver en lo que se llama forma cerrada — una solución expresada en un número finito de funciones elementales como suma, multiplicación, raíces, exponenciales, logaritmos y funciones trigonométricas. Cuando combinamos un término como incógnita con un término trigonométrico como sin⁡(incógnita), aislando incógnita puede volverse fundamentalmente imposible. En otras palabras, este tipo de ecuaciones mixtas a menudo carecen de una solución en forma cerrada.

Está bien. Del gráfico sabemos que existe una solución. SymPy puede acercarnos arbitrariamente a esa solución utilizando métodos numéricos. Usamos SymPy nsolve():

import sympy as sym
from sympy import pi, sin
from sympy.abc import x
c = 0.967
time_earth_years = 20
orbital_period_earth_years = 76
time_radians = (time_earth_years / orbital_period_earth_years) * 2 * pi
equation = x - c * sin(x) - time_radians
initial_guess = 1.0   # Initial guess for the numerical solver
position_radians = sym.nsolve(equation, x, initial_guess)
print(f"After {time_earth_years} Earth years, the comet will travel {position_radians:.4f} radians ({position_radians * 180 / pi:.2f}°) along its orbit.")

Que informes:

After 20 Earth years, the comet will travel 2.3449 radians (134.35°) along its orbit.

Podemos resumir los resultados en una tabla:

¿Estamos seguros de que no existe una solución de forma cerrada? Agregamos un signo de interrogación a nuestra respuesta “No”. Esto nos recuerda que el fracaso de SymPy no es una prueba matemática de que no existe una solución de forma cerrada. Etiquetamos la última columna “A Numérico” para recordarnos que representa una solución numérica. Podría haber más.