Part 1: Root Finding

Learning Objectives¶

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 a differentiable potential energy is stationary, or where competing effects balance. Mathematically, these all reduce to finding roots of equations: solving $f(x) = 0$.

Many equations in computational physics have no closed-form solution. For example, a stellar model must satisfy hydrostatic equilibrium:

This is a differential equation, not itself a scalar root equation. Once we discretize a stellar model or impose a boundary condition, root finding can determine an unknown central value, radius, or model parameter. We need to distinguish solving the ODE from solving that resulting residual equation.

Building Intuition: The Geometry of Root Finding¶

Before diving into algorithms, let’s understand the geometry. Finding roots is about locating zeros; a curve can cross the x-axis or touch it without crossing. Different methods use different geometric insights:

1. Bracketing methods (bisection) maintain an interval with a sign change.

2. Local methods (Newton) use a tangent-line model near an iterate.

3. Interpolation methods (secant) replace the analytic derivative with a slope from two evaluations.

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 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)=0$: return the sampled zero; otherwise stop only when the half-width meets the requested position tolerance

   • 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 33 interval halvings suffice, followed by the new midpoint, provided continuity and reliable function signs hold. Count midpoint evaluations separately from interval updates.

Implementation¶

This reading uses a compact reference interface to illustrate bracket reduction. It differs from the Project 3 contract:

Use the project brief for submitted interfaces, history recording, and failure behavior. The example below is a reading reference, not a replacement module.

import math
import numpy as np


def bisection(f, a, b, tol=1e-10, max_iter=100):
    """Return a bracket midpoint with absolute position tolerance tol.

    f must be continuous, finite, scalar-valued, and reliably signed.
    tol has x's units; it is not a tolerance on the residual f(x).
    """
    if not (math.isfinite(a) and math.isfinite(b) and a <= b):
        raise ValueError("Need finite ordered endpoints")
    if not math.isfinite(tol) or tol <= 0:
        raise ValueError("tol must be positive and finite")
    if isinstance(max_iter, (bool, np.bool_)) or not isinstance(max_iter, (int, np.integer)) or max_iter <= 0:
        raise ValueError("max_iter must be a positive integer")

    def evaluate(x):
        value = float(f(x))
        if not math.isfinite(value):
            raise ValueError("f must return finite values")
        return value

    fa, fb = evaluate(a), evaluate(b)
    if fa == 0.0:
        return a
    if fb == 0.0:
        return b
    # Comparing signs avoids overflow/underflow from multiplying values.
    if (fa < 0) == (fb < 0):
        raise ValueError("Need opposite signs at the endpoints")

    for _ in range(max_iter):
        c = a / 2 + b / 2   # Avoid overflow in a + b.
        if max(c - a, b - c) <= tol:
            return c
        if c == a or c == b:
            raise RuntimeError("No representable interior point; tolerance unmet")
        fc = evaluate(c)
        if fc == 0.0:
            return c
        if (fa < 0) != (fc < 0):
            b, fb = c, fc
        else:
            a, fa = c, fc
    raise RuntimeError(f"Failed to converge in {max_iter} iterations")

An exactly zero computed function value is accepted here. If evaluating the function can underflow to zero or produce unreliable signs, the bracket is no longer certified by these values. Rescale the evaluation or use higher precision; changing the stopping criterion cannot repair the function evaluation.

Predict → Run → Explain: small residual or well-located root?¶

Multiplying a function by a positive constant does not move its root. Predict whether the same position tolerance should give the same answer for both rows.

for scale in (1.0, 1e-12, 1e-200):
    f_scaled = lambda x: scale * (x - 0.9)
    root = bisection(f_scaled, 0.0, 1.0, tol=1e-10)
    print(f"scale={scale:.0e}: root={root:.12f}, "
          f"position error={abs(root-0.9):.3e}, residual={abs(f_scaled(root)):.3e}")

scale=1e+00: root=0.899999999965, position error=3.492e-11, residual=3.492e-11
scale=1e-12: root=0.899999999965, position error=3.492e-11, residual=3.492e-23
scale=1e-200: root=0.899999999965, position error=3.492e-11, residual=3.492e-211

The position errors agree although the residuals differ by many orders of magnitude. Near a simple root, $f(x)\approx f'(r)(x-r)$: a residual measures position error only after accounting for the slope and the function’s units. If $|f'|\ge m>0$ on the relevant interval, the mean value theorem gives $|x-r|\le |f(x)|/m$. A bracket gives a position bound without estimating that slope.

Common Failure Modes and Fixes¶

Method 2: Newton-Raphson — using a local slope¶

The Geometric Insight¶

Newton’s method uses calculus to accelerate convergence. The key insight: near any point, a smooth function looks approximately linear. We follow the tangent line to where it crosses the axis. Sufficiently near a simple root, under suitable smoothness conditions, this can provide a much more accurate estimate. A tangent step from an arbitrary starting point need not improve the estimate.

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:

For a sufficiently smooth function, a simple root with $f'(r)\ne0$, and a sufficiently close initial guess, this shows quadratic convergence! The error is squared each iteration, meaning the number of correct digits approximately doubles.

When Newton Fails¶

The tangent model can fail for several distinct reasons:

1. Zero derivative away from a root: If $f'(x_n)=0$ and $f(x_n)\ne0$, the horizontal tangent has no zero

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

3. Cycles: Can oscillate between points without converging

4. Sensitive root: A small $|f'(r)|$ amplifies a given absolute function error; interpret that error together with units and function scaling

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¶

For a sufficiently smooth function with a simple root and sufficiently close starting points, 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.

Where the golden ratio comes from¶

The secant slope contains information from two previous errors. Expanding $f(r+e)=a e+b e^2+O(e^3)$, with $a=f'(r)\ne0$ and $b=f''(r)/2$, gives

Substituting in the secant update and cancelling the leading linear term gives

Suppose this nonzero leading term controls the iteration and the errors obey $|e_k|\sim K|e_{k-1}|^p$. Then $|e_{k-1}|$ scales like $|e_k|^{1/p}$, so

Thus $p^2-p-1=0$, and the positive root is $(1+\sqrt5)/2$. The result comes from reusing two iterates, not from a universal empirical ranking of solvers. If the leading coefficient vanishes, the root is repeated, or rounding dominates, this particular asymptotic calculation does not describe the measured sequence.

Comparison of Convergence Rates¶

The bracket bound can contract even when successive midpoint errors do not. The computed example below uses $f(x)=x^2-2$; it compares actual errors with the bisection bound, not with residuals or update sizes.

There is no universal absolute cutoff for a “small” secant denominator: multiplying the equation by a constant scales that denominator while leaving the exact secant step unchanged. Compare evaluation uncertainty with the local function-value scale and inspect the proposed step. If function error dominates the difference, report an unreliable slope or use a justified safeguard. Project 3 requires explicit handling of zero denominators and nonfinite values; it does not prescribe an extra arbitrary near-zero threshold.

Common Failure Modes and Fixes¶

Backtracking: choose how far to follow Newton¶

Newton proposes a direction $p=-f(x)/f'(x)$. A line search chooses a step length $\alpha$ along that direction:

The full Newton step uses $\alpha=1$. Backtracking starts there and shortens the step when it fails an acceptance test. This controls a local proposal; it does not preserve a root bracket or select a physically meaningful root.

Define progress before testing a step¶

For a scalar equation, use the merit function $\phi(x)=\tfrac12 f(x)^2$. With a differentiable $f$ and finite nonzero $f'(x)$, an exact Newton direction satisfies

at a nonroot. It is a descent direction locally, although a full step can increase the merit function. The Armijo sufficient-decrease test is

Choose $0<c<1/2$ here; the example uses $c=10^{-4}$ and halves $\alpha$ after each rejection. These are illustrative algorithm settings, not accuracy tolerances. Under the differentiability and descent assumptions, sufficiently small positive steps satisfy this test in exact arithmetic. A finite search budget, floating-point stagnation, or unreliable evaluations can still prevent an implementation from accepting a step.

Predict, compare, explain: one accepted step¶

For $f(x)=\arctan x$, start at $x=2$. Predict whether the full Newton step reduces $|f|$, then compare these candidate trials (including an extra quarter-step for comparison). This is one step, not a complete solve.

Backtracking accepts the first passing trial, not necessarily the best point on the line. Explain why an accepted step neither identifies a root nor preserves a bracket. No line-search implementation is assigned.

The full step reaches approximately −3.536 and increases the residual magnitude. One halving reaches approximately −0.768 and satisfies sufficient decrease. This is progress, not a converged root. A complete solver still needs an outer iteration budget, stopping rules, diagnostics, and a physical-domain check. Domain-restricted functions also need an explicit feasible-trial policy; this example does not turn arbitrary domain errors into successful steps.

Three different ways to control a proposal¶

FNC’s quasi-Newton chapter develops Levenberg’s method for nonlinear systems: regularization changes the proposed step, and residual reduction governs acceptance. It is related to step control but is not the scalar backtracking algorithm above. For a production Newton/line-search interface, see PETSc. Levenberg’s method is a preview for later nonlinear least squares, not another algorithm to implement this week.

Jacobians and local fitting · Week 5 companion¶

For observations $y_i$ and a model $g(t_i;\theta)$, define residuals $r_i=g(t_i;\theta)-y_i$. With $m$ observations and $n$ parameters, the Jacobian $J$ has shape $m\times n$ and entries $J_{ij}=\partial r_i/\partial\theta_j$. Each row describes one observation; each column describes sensitivity to one parameter.

For a small parameter update $p$, $r(\theta+p)\approx r(\theta)+Jp$. The product predicts a change, not the new residual itself.

For $g(t;A,k)=Ae^{-kt}$ at $A=2$, $k=\ln2$, and $t=0,1,2$ (scaled, dimensionless coordinates),

The update $p=(0.10,0.02)^{\mathsf T}$ predicts changes $(0.10,0.03,0.005)^{\mathsf T}$. Actual changes are approximately $(0.10,0.029209,0.004414)^{\mathsf T}$. Reducing the step by ten reduces the largest discrepancy by about 99: the linear approximation improves locally.

The least-squares loss is $\Phi=\tfrac12\lVert r\rVert^2$, where $\lVert r\rVert^2=\sum_i r_i^2$. Transpose, written $\mathsf T$, exchanges rows and columns. The loss gradient (its parameter partial derivatives) is $J^{\mathsf T}r$.

Gauss–Newton minimizes $\tfrac12\lVert r+Jp\rVert^2$. A Levenberg step adds a penalty on the current update:

Setting $\nabla_p q_\lambda=J^{\mathsf T}(r+Jp)+\lambda p=0$ gives

Here $I$ is the identity matrix and $\lambda>0$ controls damping. The parameters must be scaled sensibly; Marquardt variants incorporate scaling. This is a penalty on an update, not a fixed prior on the parameters. Check actual loss reduction before accepting a trial. The displayed normal equations explain the method; production solvers can use QR or SVD to avoid the conditioning penalty of explicitly forming $J^{\mathsf T}J$.

In regression and calibration, noisy data generally cannot all have zero residual. This differs from solving a square system $F(u)=0$, where a nonsingular (invertible) Jacobian gives the Newton proposal $J_Fp=-F$. Damping can control a poor step but cannot create information absent from the observations. For example, early observations of a decay curve weakly constrain $k$ because $\partial g/\partial k=-At e^{-kt}$ is small near $t=0$.

Explain: why can a tiny damped step coexist with a large residual? No implementation of these advanced fitting methods is required for Project 3.

Hybrid Methods: preserve a bracket while using local models¶

A safeguarded method uses interpolation or a derivative when the proposed step is acceptable and falls back to bracket reduction when it is not.

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

• Safeguards maintain convergence from a valid continuous bracket; often needs fewer evaluations than bisection

Practical Strategy for Robust Root Finding¶

1. Find a continuous sign-change bracket by physics or a scan, then use bisection
2. Propose a local Newton or interpolation step while retaining the bracket
3. Reject unsafe proposals and enforce regular bracket reduction
4. Monitor bracket width, residual scale, stagnation, and evaluation cost

All three hard-example runs reach the stated 20-iteration cap. Their final position errors are different; none should be reported as satisfying its stopping contract merely because one plotted error is small.

Investigation: which root estimate is justified?¶

Suppose three methods report nearly the same root. One stopped because the residual was small, one because the last step was small, and one because a bracket was narrow. These are different observations. Which one controls the position of the root, and what additional information does each require?

Return to the cubic $f(x)=x^3-2x-5$ from the opening figure. It is continuous, and $f'(x)=3x^2-2\ge10$ on $[2,2.5]$, so a sign-changing bracket there contains one root. Before inspecting the recorded histories, predict:

1. Must a bisection midpoint become more accurate at every iteration?

2. Can Newton make a small step without a small position error?

3. Does one Newton iteration cost the same as one secant iteration?

4. Which conclusion survives multiplying the whole equation by 10^-12?

For this problem, a bracket midpoint has position error at most the half-width. At a candidate in $[2,2.5]$, the derivative bound also gives $|x-r|\le |f(x)|/10$. The last step $|x_k-x_{k-1}|$ is neither of these bounds. It measures what the algorithm just did, not its distance from the answer. Stagnation can make it zero even when a target remains unmet.

The final recorded rows from the cubic comparison are:

A printed zero position error means agreement with the binary64 reference, not an exact representation of the mathematical root. A sampled zero residual also requires the function-evaluation caveat discussed with bisection.

The comparison figure reports real errors against an independent decimal reference. Those errors are available for this teaching example; the method cannot normally inspect them. Its stopping condition must use information available during the calculation. Count calls to $f$ and $f'$ separately: a derivative may be cheaper, more expensive, or computed together with $f$. An iteration count by itself does not settle the cost comparison.

Transfer: an accurate solve of an uncertain function¶

Now suppose the function values themselves have a bounded additive uncertainty $\eta$. This is a different experiment from multiplying an exact equation by a constant. Keep $x$ and function values dimensionless, and consider

The unperturbed root is 1. Any perturbed zero obeys $|x-1|\le\eta/|s|$. A smaller slope therefore permits a wider range of root positions for the same absolute uncertainty. A tighter iteration tolerance cannot distinguish among those possible input functions.

# These are possible-input bounds, not errors measured from a root solver.
slopes = np.array([1.0, 0.05])
function_uncertainty = 0.002
root_half_ranges = function_uncertainty / np.abs(slopes)
for slope, half_range in zip(slopes, root_half_ranges):
    print(f"slope={slope:g}: possible roots in [{1-half_range:.3f}, {1+half_range:.3f}]")

slope=1: possible roots in [0.998, 1.002]
slope=0.05: possible roots in [0.960, 1.040]

Predict what happens to this range if we halve the solver tolerance, halve $\eta$, or multiply both $f_s$ and its uncertainty by 100. Only one of those changes improves the information supplied about the root.

A repeated root makes the issue sharper. For $f(x)=(x-1)^2$, a downward constant perturbation of size $\eta>0$ produces two zeros, $1\pm\sqrt{\eta}$, while an upward perturbation removes the real zero. There is no nonzero derivative lower bound at the original root. For the unperturbed function, Newton’s exact recurrence is

A tiny residual is especially easy to obtain here because it is the square of the position error. For example, residual 10^-12 corresponds to position error 10^-6, not 10^-12.

Further reading: conditioning, residuals, and multiple roots.

Application: Units and Scaling¶

Choose scales that make the variables and tolerances interpretable. Non-dimensionalization can avoid extreme arithmetic and improve numerical balancing, but it does not remove intrinsic sensitivity:

Dimensional form: Finding radius where $\rho(r) = 10^{-7}$ g/cm$^3$ with $r$ in cm Scaled form: Define $\tilde{r} = r/R_{\odot}$, $\tilde{\rho} = \rho/\rho_0$ and solve $\tilde{\rho}(\tilde{r}) = 10^{-7}\,\mathrm{g\,cm^{-3}}/\rho_0$ (equal to 1 only if that target is the chosen $\rho_0$)

Track the units of the position tolerance and transform it with the variable. Rescaling an equation changes its residual values; it does not by itself make a root better determined by the physical data.

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$ is the constant-density Lane–Emden limit: $\theta=1-\xi^2/6$, so $\xi_1=\sqrt6$.

• $n=1$ has $\theta=\sin\xi/\xi$, with its central value defined by the limit 1; $\xi_1=\pi$. This is a $P=K\rho^2$ polytrope, not an isothermal core.

• $n=3$ has $\xi_1\approx6.897$ and appears in idealized radiation-supported and relativistically degenerate models.

• $n=5$ has $\theta=(1+\xi^2/3)^{-1/2}$: it approaches zero only at infinite radius, although its total mass is finite. It is not a minimum-mass-star model.

These solutions use $\theta(0)=1$ and $\theta'(0)=0$. Reference: analytic polytropes.

A computed Lane–Emden trajectory can supply both $\theta$ and its derivative near the surface. Root refinement then depends on the accuracy of that trajectory as well as on the root solver. Tightening a root tolerance cannot remove the ODE or interpolation error. Preserve a credible surface bracket and refine the underlying trajectory before claiming additional digits.

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:

For an ellipse, $0\le e<1$, so $f'(E)=1-e\cos E\ge1-e>0$. The derivative can become small near periapsis as $e$ approaches 1, but it cannot vanish in this domain.

Initial Guess Strategy¶

One heuristic on a reduced mean-anomaly interval is to start with $E_0=M$ for $e<0.8$ and $E_0=\pi$ otherwise. The value 0.8 is not a convergence threshold, and the first regime is not uniformly nearly circular. Compare both guesses over the project’s sweep and report failures and iteration counts. If reducing $M$ modulo $2\pi$, restore its branch in the returned $E$.

Convergence Analysis¶

The absolute sensitivity to mean anomaly is

For $e=0$, $E=M$ exactly. For $e=0.9$ the upper bound is 10, and for $e=0.99$ it is 100; these are bounds, not the sensitivity at every orbital phase. Iteration counts also depend on $M$, the initial guess, and the stopping criterion. Normalize $M$ to $[0,2\pi]$ for the initial-guess choices above. Because the residual is increasing and changes sign on $[0,2\pi]$ for an interior $M$, bisection provides a useful independent check on Newton’s answer.

Practical Algorithm for Finding Initial Brackets¶

To find initial brackets for bisection:

def find_brackets(f, x_min, x_max, n_scan=100):
    """Return sign-change intervals and (x, x) pairs for sampled zeros."""
    if isinstance(n_scan, (bool, np.bool_)) or not isinstance(n_scan, (int, np.integer)) or n_scan <= 0:
        raise ValueError("n_scan must be a positive integer")
    if not (np.isfinite(x_min) and np.isfinite(x_max) and x_min < x_max):
        raise ValueError("Need finite increasing scan limits")
    nodes = np.linspace(x_min, x_max, n_scan + 1)
    values = np.array([float(f(x)) for x in nodes])
    if not np.all(np.isfinite(values)):
        raise ValueError("Non-finite value in scan")
    brackets = []
    for i, (x, fx) in enumerate(zip(nodes, values)):
        if fx == 0.0:
            brackets.append((x, x))
        elif i > 0 and values[i-1] != 0.0 and (values[i-1] < 0) != (fx < 0):
            brackets.append((nodes[i-1], x))
    return brackets

print(find_brackets(lambda x: x - 0.5, 0.0, 1.0))

[(np.float64(0.5), np.float64(0.5))]

A scan can miss an even-multiplicity root between samples or several crossings within one cell. A sampled zero is reported once, as a zero-width pair; the reading’s bisection accepts it. Project 3’s bisect requires a < b, so retain these sampled roots separately and send only positive-width brackets to that solver. Neither a sign scan nor bisection proves continuity or verifies that a computed zero did not arise from underflow.

Debugging Root Finding¶

Common root-finding bugs and their diagnosis:

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

You can now compare three kinds of evidence: a shrinking bracket, a local iteration model, and a residual interpreted through sensitivity. A method choice should identify which of these supports the reported root.

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 synthesis returns to a problem in which an integral defines the root residual. The two algorithms have separate jobs, but their errors contribute to one final answer.

Next: Part 2 - Quadrature