Part 1: Root Finding - Where Physics Reaches Equilibrium

Learning Outcomes¶

By the end of this section, you will be able to:

• Implement bisection, Newton-Raphson, and secant methods from scratch

• Analyze convergence rates (linear, superlinear, quadratic) mathematically

• Diagnose method failures and design robust hybrid approaches

• Apply root finding to Kepler’s equation, Lagrange points, and stellar structure

• Predict iteration counts and verify convergence behavior

The Fundamental Problem¶

In physics, equilibrium occurs where the net force vanishes, where energy is minimized, or where competing effects balance. Mathematically, these all reduce to finding roots of equations: solving $f(x) = 0$.

But here’s the challenge: most astrophysical equations are transcendental equations that cannot be solved analytically. Consider finding where a star’s pressure balances gravity:

With realistic equations of state and composition gradients, this becomes a transcendental equation with no closed-form solution. We need numerical methods.

Building Intuition: The Geometry of Root Finding¶

Before diving into algorithms, let’s understand the geometry. Finding roots is about discovering where a curve crosses the x-axis. Different methods use different geometric insights:

1. Bracketing methods (Bisection): Trap the root between two points and squeeze—when you absolutely must find the root, even if slowly

2. Local approximation methods (Newton): Follow the tangent line to the axis—when you need maximum speed and have derivatives

3. Interpolation methods (Secant): Approximate the curve with simpler functions—the practical compromise when derivatives are expensive

Each approach has trade-offs between:

• Reliability: Will it always find a root?

• Speed: How many iterations to converge?

• Requirements: Do we need derivatives? Good initial guesses?

Method 1: Bisection - The Reliable Workhorse¶

The Mathematical Foundation¶

The bisection method is based on the Intermediate Value Theorem:

If $f$ is continuous on $[a,b]$ and $f(a) \cdot f(b) < 0$, then there exists at least one root $r \in (a,b)$ where $f(r) = 0$.

The algorithm repeatedly halves the interval, keeping the half that contains the root.

Algorithm Analysis¶

Starting with interval $[a_0, b_0]$ where $f(a_0) \cdot f(b_0) < 0$:

1. Compute midpoint: $c_n = \frac{a_n + b_n}{2}$

2. Evaluate $f(c_n)$

3. Update interval:

   • If $|f(c_n)| < \epsilon_{tol}$: Found root to desired accuracy

   • If $f(a_n) \cdot f(c_n) < 0$: root in $[a_n, c_n]$, so set $b_{n+1} = c_n$, $a_{n+1} = a_n$

   • Otherwise: root in $[c_n, b_n]$, so set $a_{n+1} = c_n$, $b_{n+1} = b_n$

After $n$ iterations, the interval width is:

The error in our approximation is bounded by half the interval width:

This is linear convergence with rate $\frac{1}{2}$.

Iterations Required¶

To achieve error $< \epsilon$, we need:

Solving for $n$:

For example, to find a root in $[0, 1]$ to 10 decimal places ($\epsilon = 10^{-10}$):

So we need 33 iterations—slow but guaranteed!

Implementation¶

Here’s a complete bisection implementation to understand the algorithm:

def bisection(f, a, b, tol=1e-10, max_iter=100):
    """Find root of f in [a,b] by bisection"""
    # Check initial bracket
    fa, fb = f(a), f(b)
    if fa * fb >= 0:
        raise ValueError("Need f(a) and f(b) with opposite signs")
    if abs(fa) < tol: return a
    if abs(fb) < tol: return b
    
    for i in range(max_iter):
        c = (a + b) / 2
        fc = f(c)
        
        if abs(fc) < tol or abs(b - a) < 2*tol:
            return c
            
        if fa * fc < 0:
            b, fb = c, fc
        else:
            a, fa = c, fc
            
    raise RuntimeError(f"Failed to converge in {max_iter} iterations")

Common Failure Modes and Fixes¶

Method 2: Newton-Raphson - The Speed Demon¶

The Geometric Insight¶

Newton’s method uses calculus to accelerate convergence. The key insight: near any point, a smooth function looks approximately linear. We can follow the tangent line to where it crosses the axis, getting much closer to the root in a single step.

Mathematical Derivation¶

Starting at point $x_n$, the tangent line to $f(x)$ has equation:

This line crosses the x-axis (where $y = 0$) at:

Solving for $x_{n+1}$:

This is the Newton-Raphson iteration formula.

Convergence Analysis¶

To understand the convergence rate, we analyze the error $e_n = x_n - r$ where $r$ is the true root.

Using Taylor series around $r$:

Since $f(r) = 0$:

Similarly:

Through algebraic manipulation of the Newton iteration formula:

This shows quadratic convergence! The error is squared each iteration, meaning the number of correct digits approximately doubles.

When Newton Fails¶

Despite its speed, Newton’s method can fail spectacularly:

1. Zero derivative: If $f'(x_n) = 0$, the tangent is horizontal and never crosses the axis

2. Poor initial guess: May diverge or converge to wrong root

3. Cycles: Can oscillate between points without converging

4. Ill-conditioned: When $|f'(r)|$ is very small, round-off errors dominate

Common Failure Modes and Fixes¶

Method 3: Secant - The Practical Compromise¶

Motivation: No Derivatives Required¶

Newton’s method requires $f'(x)$, but what if:

• The derivative is expensive to compute?

• We only have $f(x)$ as a black box?

• The function comes from experimental data?

The secant method approximates the derivative using a finite difference:

The Iteration Formula¶

Substituting this approximation into Newton’s formula:

Geometrically, we’re replacing the tangent line with a secant line through two points.

Convergence Analysis¶

The convergence of the secant method is more complex. Through detailed analysis (involving difference equations), one can show that the order of convergence is:

This is the golden ratio! The error satisfies:

This superlinear convergence is slower than Newton’s quadratic rate but faster than bisection’s linear rate.

Comparison of Convergence Rates¶

Common Failure Modes and Fixes¶

Hybrid Methods: Best of All Worlds¶

In practice, we often combine methods to leverage their strengths:

Brent’s Method¶

Combines bisection’s reliability with inverse quadratic interpolation’s speed:

• Maintains a bracketing interval like bisection

• Tries fast interpolation when possible

• Falls back to bisection if interpolation fails

• Never worse than bisection, often much faster

Practical Strategy for Robust Root Finding¶

1. Start with bisection to bracket the root
2. Switch to Newton/Secant once close enough  
3. Fall back to bisection if Newton diverges
4. Monitor condition number throughout

Application: Units and Scaling¶

Always non-dimensionalize your equations before root finding. Working with $x \sim 1$ rather than $x \sim 10^{23}$ improves numerical stability:

Bad: Finding radius where $\rho(r) = 10^{-7}$ g/cm³ with $r$ in cm Good: Define $\tilde{r} = r/R_{\odot}$, $\tilde{\rho} = \rho/\rho_0$ and solve $\tilde{\rho}(\tilde{r}) = 1$

This keeps all numbers near unity, avoiding overflow/underflow and improving conditioning.

Physical Example: Stellar Structure Equilibrium¶

In stellar interiors, hydrostatic equilibrium requires:

For a polytrope with $P = K\rho^{1+1/n}$, the Lane-Emden equation describes stellar structure:

where $\xi$ is dimensionless radius and $\theta$ is related to density.

The stellar surface occurs where $\theta(\xi_1) = 0$—a root-finding problem!

For different polytropic indices:

• $n = 0$ (constant density): $\xi_1 = \pi$ exactly

• $n = 1$ (isothermal core): $\xi_1 = \pi/\sqrt{2}$

• $n = 3$ (Eddington standard model): $\xi_1 \approx 6.897$

• $n = 5$ (minimum mass star): $\xi_1 \to \infty$ (infinite radius!)

Newton’s method excels here because we can compute $\theta$ and $d\theta/d\xi$ from the differential equation. The convergence is rapid even for the sensitive $n = 3$ case.

Worked Example: Kepler’s Equation¶

The most famous root-finding problem in astronomy is Kepler’s equation, which relates mean anomaly $M$ (proportional to time) to eccentric anomaly $E$ (related to position):

where $e$ is the orbital eccentricity and both $M$ and $E$ are in radians.

Why This Is Hard¶

• Transcendental equation (no algebraic solution)

• Must be solved millions of times in orbit propagation

• High eccentricity makes convergence difficult

• Conditioning worsens as $e \to 1$ (parabolic orbit)

Newton’s Method Applied¶

Define $f(E) = E - e\sin(E) - M$

Then $f'(E) = 1 - e\cos(E)$

Newton iteration:

Note: $f'(E) = 0$ when $e\cos(E) = 1$, which can happen for high eccentricity!

Initial Guess Strategy¶

Good initial guesses are crucial:

• For $e < 0.8$: Use $E_0 = M$ (nearly circular)

• For $e \geq 0.8$: Use $E_0 = \pi$ (highly eccentric)

Convergence Analysis¶

Practical Algorithm for Finding Initial Brackets¶

To find initial brackets for bisection:

def find_brackets(f, x_min, x_max, n_scan=100):
    """Scan interval for sign changes"""
    dx = (x_max - x_min) / n_scan
    brackets = []
    
    x_prev = x_min
    f_prev = f(x_prev)
    
    for i in range(1, n_scan + 1):
        x = x_min + i * dx
        fx = f(x)
        
        if f_prev * fx < 0:
            brackets.append([x_prev, x])
        
        x_prev, f_prev = x, fx
    
    return brackets

Debugging Root Finding¶

Common root-finding bugs and their diagnosis:

Bridge to Part 2: From Finding Zeros to Measuring Areas¶

You now command three powerful methods for finding where functions equal zero. Each has its place: bisection for reliability, Newton for speed, secant for practicality.

But finding equilibrium points is only half the story. Next, we tackle integration—measuring quantities across space and time. You’ll discover how the same principles of approximation and error analysis apply to computing areas under curves, volumes of regions, and integrals over high-dimensional spaces.

The connection runs deep: integration is the inverse of differentiation, and many integration methods rely on root finding. The error analysis skills you’ve developed here will guide your understanding of quadrature methods.

Next: Part 2 - Quadrature