Why can we solve some equations easily, while others seem impossible? And another thing: why is this knowledge hidden from us?
As data scientists, applied scientists, and engineers, we often create mathematical models. For example, consider the model: y = x². Given a value for x, we can apply it forward to compute y. For instance, if x = 3, then y = 9.
We can also apply this model backward. Starting with y = x², we rearrange to solve for x: x = ±√y. If y = 9, then x = ±3. The expression x = ±√y is an example of a closed-form solution — an expression that uses a finite combination of standard operations and functions.
However, not all models are so straightforward. Sometimes, we encounter equations where we can’t simply “solve for x” and get a closed-form expression. In such cases, we might hear, “That’s not solvable — you need numerical methods.” Numerical methods are powerful. They can provide precise approximations. Still, it frustrates me (and perhaps you) that no one ever seems to explain when closed-form solutions are possible and when they aren’t.
The great Johannes Kepler shared our frustration. When studying planetary motion, he created this model:
This equation converts a body’s position along its orbit (x) into its time along the orbit (y). Kepler sought a closed-form solution for x to turn time into a position. However, even 400 years later, the best we have are numerical methods.
In this article, we’ll build intuition about when to expect a closed-form solution. The only way to determine this rigorously is by using advanced mathematics — such as Galois theory, transcendental number theory, and algebraic geometry. These topics go far beyond what we, as applied scientists and engineers, typically learn in our training.
Instead of diving into these advanced fields, we’ll cheat. Using SymPy, a Python-based computer algebra system, we’ll explore different classes of equations to see which it can solve with a closed-form expression. For completeness, we’ll also apply numerical methods.
We’ll explore equations that combine polynomials, exponentials, logarithms, and trigonometric functions. Along the way, we’ll discover specific combinations that often resist closed-form solutions. We’ll see that if you want to create an equation with (or without) a closed-form solution, you should avoid (or try) the following:
- Fifth degree and higher polynomials
- Mixing x with exp(x) or log(x) — if Lambert’s W function is off-limits
- Combining exp(x) and log(x) within the same equation
- Some pairs of trigonometric functions with commensurate frequencies
- Many pairs of trigonometric functions with non-commensurate frequencies
- Mixing trigonometric functions with x, exp(x), or log(x)
Aside 1: I’m not a mathematician, and my SymPy scripts are not higher mathematics. If you find any mistakes or overlooked resources, forgive my oversight. Please share them with me, and I’ll gladly add a note.
Aside 2: Welch Lab’s recent video, Kepler’s Impossible Equation, reminded me of my frustration about not knowing when an equation can be solved in a closed form. The video sparked the investigation that follows and provides our first example.
Imagine you are Johannes Kepler’s research programmer. He has created the following model of orbital motion:
y = x −c sin(x)
where:
- x is the body’s position along its orbit. We measure this position as an angle (in radians). The angle starts at 0 radians when the body is closest to the Sun. When the body has covered ¼ of its orbit’s distance, the angle is π/2 radians (90°). When it has covered half of its orbit’s distance, the angle is π (180°), and so on. Recall that radians measure angles from 0 to 2π rather than from 0 to 360°.
- c is the orbit’s eccentricity, ranging from 0 (a perfect circle) to just under 1 (a highly elongated ellipse). Suppose Kepler has observed a comet with c = 0.967.
- y is the body’s time along its orbit. We measure this time as an angle (in radians). For instance, if the comet has an orbital period of 76 Earth years, then π/2 (90°) corresponds to ¼ of 76 years, or 19 years. A time of π (180°) corresponds to ½ of 76 years, or 38 years. A time of 2π (360°) is the full 76-year orbital period.
This diagram shows the comet’s position at π/2 radians (90°), which is ¼ of the way along its orbit:
Kepler asks for the time when the comet reaches position π/2 radians (90°). You create and run this Python code:
import numpy as npdef 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.")
You report back to Kepler:
It takes approximately 7.30 Earth years for the comet to move from 0 to π/2 radians.
Aside: The comet covers 25% of its orbit distance in under 10% of its orbital period because it speeds up when closer to the Sun.
No good deed goes unpunished. Kepler, fascinated by the result, assigns you a new task: “Can you tell me how far along its orbit the comet is after 20 Earth years? I want to know the position in radians.”
“No problem,” you think. “I’ll just use a bit of high school algebra.”
First, you convert 20 Earth years into radians:
- time_radians = (20 / 76) × 2π = (10 / 19)π
Next, you rearrange Kepler’s equation, setting it equal to 0.
- x − 0.967 sin(x) − (10 / 19)π = 0
Now you want to find the value of x that makes this equation true. You decide to graph the equation to see where it crosses zero:
import numpy as np
import matplotlib.pyplot as pltc = 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()
So far, so good. The graph shows that a solution for x exists. But when you try to rearrange the equation to solve for x using algebra, you hit a wall. How do you isolate x when you have a combination of x and sin(x)?
“That’s okay,” you think. “We’ve got Python, and Python has the SymPy package,” a powerful and free computer algebra system.
You pose the problem to SymPy:
# Warning: This code will fail.
import sympy as sym
from sympy import pi, sin
from sympy.abc import xc = 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)
Unfortunately, it replies with an error:
NotImplementedError: multiple generators [x, sin(x)]
No algorithms are implemented to solve equation x - 967*sin(x)/1000 - 10*pi/19
SymPy is quite good at solving equations, but not all equations can be solved in what’s called closed form — a solution expressed in a finite number of elementary functions such as addition, multiplication, roots, exponentials, logarithms, and trigonometric functions. When we combine a term such as x with a trigonometric term like sin(x), isolating x can become fundamentally impossible. In other words, these types of mixed equations often lack a closed-form solution.
That’s okay. From the graph, we know a solution exists. SymPy can get us arbitrarily close to that solution using numerical methods. We use SymPy’s nsolve():
import sympy as sym
from sympy import pi, sin
from sympy.abc import xc = 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.")
Which reports:
After 20 Earth years, the comet will travel 2.3449 radians (134.35°) along its orbit.
We can summarize the results in a table:
Are we sure there is not a closed-form solution? We add a question mark to our “No” answer. This reminds us that SymPy’s failure is not a mathematical proof that no closed-form solution exists. We label the last column “A Numeric” to remind ourselves that it represents one numerical solution. There could be more.