Roots and Iterative Methods

Root-finding algorithms: bisection, secant, and Newton-Raphson methods for physical transcendent equations.

What it is

The numerical resolution of nonlinear algebraic and transcendental equations $f(x) = 0$ (or multivariable vector systems $\mathbf{F}(\mathbf{x}) = \mathbf{0}$) is the foundational algorithmic process that computes physical coordinates where a continuous system reaches balance, zero flux, or resonance. Because most equations in physics lack closed-form analytic solutions, iterative schemes generate successive approximations $x_0, x_1, x_2, \dots, x_k$ converging to the true root $x^*$.

Root-finding methods fall into two primary classes:

1. Bracketing methods (guaranteed convergence):
   • Bisection Method: Grounded in Bolzano's intermediate value theorem. If $f(a) \cdot f(b) < 0$, it iteratively halves the enclosing interval. It exhibits guaranteed linear convergence, halving error bounds at each iteration ($e_{k+1} = \frac{1}{2} e_k$) and requiring $\approx 3.32$ steps per decimal digit of precision.
   • Regula Falsi (False Position): Linearly interpolates between bracket bounds, weighting the split point by function magnitudes.
2. Open methods (rapid local convergence):
   • Newton-Raphson Method: Approximates the nonlinear function locally via its tangent line: $x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$ It delivers local quadratic convergence ($|e_{k+1}| \approx C |e_k|^2$): significant correct digits roughly double each iteration, provided initial seed $x_0$ lies within the basin of attraction and $f'(x^*) \neq 0$.
   • Secant Method: Replaces analytic derivatives with finite-difference slopes between the two most recent iterates, achieving superlinear convergence of golden ratio order $\phi = \frac{1+\sqrt{5}}{2} \approx 1.618$.
   • Fixed-Point Iteration: Rewrites $f(x) = 0$ as $x = g(x)$, converging if the contraction mapping criterion $|g'(x^*)| < 1$ holds.

Why it appears in physics

Equilibrium, resonance, and eigenvalue constraints across physics generate transcendental equations:

• Astrophysics and celestial mechanics: Kepler's equation $M = E - e \sin E$ relates mean anomaly $M$ to eccentric anomaly $E$ in elliptical orbits; it has no closed-form elementary solution.
• Thermodynamics: The van der Waals cubic equation of state $(P + a/V_m^2)(V_m - b) = R T$ defines molar volume $V_m$, where real roots determine coexistence densities between vapor and liquid phases.
• Quantum optics and solid state: Wien's displacement constant arises by finding the peak of Planck's blackbody spectrum: $x e^x / (e^x - 1) - 5 = 0$; finite quantum wells impose transcendental dispersion constraints $\sqrt{z_0^2 - z^2} = z \tan(z)$ governing bound energy states.

Visual or geometric idea

• Bisection: Geometrically tracks a shrinking bracket bounding the root from both sides, guaranteeing the zero cannot escape.
• Newton-Raphson: From current coordinate $(x_k, f(x_k))$, trace the tangent line down to the horizontal axis ($y=0$); its intercept defines next iterate $x_{k+1}$. On smooth convex curves, iterations sprint toward the root.
• If $f'(x_k) \approx 0$ (near-horizontal slope), the projected tangent shoots far outside the target region, inducing explosive divergence or trapped limit cycles.

Minimal symbolic language

1. Bolzano bracket condition: $f(a) \cdot f(b) < 0 \implies \exists x^* \in (a, b) \quad \text{such that} \quad f(x^*) = 0$
2. Newton-Raphson recurrence: $x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$
3. Secant recurrence: $x_{k+1} = x_k - f(x_k) \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}$
4. Dual stopping criterion (residual and step): $|f(x_k)| < \epsilon_{\text{res}} \quad \text{and} \quad |x_{k+1} - x_k| < \epsilon_{\text{step}} (1 + |x_k|)$

What the result means

The computed root represents the physical equilibrium or state parameter (e.g. equilibrium position of a suspended catenary, plasma ignition temperature, cutoff frequency in a fiber). Iteration counts quantify algorithmic performance. When Newton-Raphson reaches target accuracy in 4 steps while Bisection takes 45, computational throughput increases tenfold.

When not to use it

1. Do not use Newton-Raphson near zero derivatives: If $f'(x) \approx 0$ near current estimates (local extrema or inflection points), division by tiny slopes causes catastrophic overflow or wild divergence.
2. Do not use Bisection on even-multiplicity roots: For tangential roots like $f(x) = (x - 2)^2 = 0$, $f(x) \geq 0$ everywhere, so no bracket satisfying $f(a) f(b) < 0$ can ever be formed.
3. Do not rely on a single naive stopping metric: Checking only $|f(x_k)| < 10^{-6}$ fails for shallow slopes (where small vertical residuals mask large coordinate errors); checking only $|x_{k+1} - x_k| < 10^{-6}$ fails on steep slopes.

Operational summary

• Plot or tabulate $f(x)$ to visually isolate roots and identify valid bracket intervals $[a, b]$ or initial seeds $x_0$.
• For absolute fail-safe convergence, use Bisection or Brent's hybrid method (combining Bisection, Secant, and inverse quadratic interpolation).
• When analytic derivatives are available and a close seed is known, deploy Newton-Raphson for quadratic speed.
• Always include an iteration safety cap ($k_{\text{max}} \sim 100$) to prevent infinite loops.
• Apply dual stopping checks verifying both $|f(x_k)|$ and step size $|\Delta x_k|$.

Procedure for Numerical Root-Finding

Before calculating

Before running root-finding routines:

• Verify physical variable domains (absolute temperature $T > 0$, molar volume $V_m > b$, mass $m > 0$).
• Plot or tabulate the function to ensure it crosses zero and locate any asymptotes, discontinuities, or flat regions in the search window.

Step 1: identify the mathematical object

Cast the physical equation into canonical form: $f(x) = 0$ Compute the analytical derivative $f'(x)$ if feasible:

• If $f'(x)$ is continuous and easily evaluated: prime candidate for Newton-Raphson.
• If $f(x)$ is a black-box simulation output or differentiating is expensive: use the Secant method or Bisection.

Step 2: choose the operation

Select the algorithm matching confidence levels:

• Guaranteed bracket $[a, b]$ with $f(a) f(b) < 0$: Bisection or Brent's method.
• Close initial estimate $x_0$ available: Newton-Raphson.

Step 3: preserve units and signs

• Iterate $x_k$ carries dimensions of the unknown observable ($[x]$).
• Function residual $f(x_k)$ carries equation dimensions ($[f]$).
• Derivative $f'(x_k)$ carries units of $[f]/[x]$.
• Step correction $\Delta x = \frac{f(x_k)}{f'(x_k)}$ strictly preserves dimensions of $[x]$, ensuring dimensional consistency in $x_{k+1} = x_k - \Delta x$.

Step 4: calculate (Newton-Raphson iterative cycle)

1. Initialize seed $x_0$, residual tolerance $\epsilon_f = 10^{-12}$, step tolerance $\epsilon_x = 10^{-12}$, and iteration cap $k_{\text{max}} = 50$.
2. For $k = 0, 1, 2, \dots$:
   • Evaluate $y_k = f(x_k)$ and $d_k = f'(x_k)$.
   • If $|d_k| < 10^{-15}$, abort: near-zero slope hazards division by zero.
   • Compute step: $\Delta x_k = y_k / d_k$.
   • Update estimate: $x_{k+1} = x_k - \Delta x_k$.
   • Test convergence: $\text{If } |f(x_{k+1})| < \epsilon_f \quad \text{and} \quad |\Delta x_k| < \epsilon_x (1 + |x_{k+1}|), \quad \text{terminate and return } x_{k+1}.$
3. If $k = k_{\text{max}}$ is reached without meeting criteria, flag a convergence failure.

Step 5: interpret the result

• Substitute root $x^*$ back into the original physics equation to confirm physical residual insignificance.
• Confirm physical viability (e.g. for a van der Waals cubic with three real roots, the smallest is the compressed liquid volume, the middle root is mechanically unstable with negative compressibility, and the largest is the gas phase).
• Report final figures consistent with the physical parameter precision.

Self-check checklist

• Is initial seed $x_0$ within the valid physical domain?
• Did I confirm $f'(x_0) \neq 0$?
• Does the routine enforce a maximum iteration limit?
• Did I check both residual $|f(x)|$ and coordinate step $|\Delta x|$?
• If multiple roots exist, did the solver converge to the physically meaningful branch?

Guided Example: Solving Kepler's Orbital Equation in Celestial Mechanics

Physical Situation

A communications satellite travels along an eccentric Earth orbit with orbital eccentricity $e = 0.600$. One-third of an orbit past perigee, its mean anomaly is $M = \frac{2\pi}{3} \approx 2.0943951 \text{ rad}$. We calculate the eccentric anomaly $E$ to fix the satellite's exact orbital position using Kepler's transcendent equation: $M = E - e \sin E$ This equation cannot be solved algebraically for $E$.

Mathematical Translation

Formulating the root-finding function: $f(E) = E - 0.600 \sin E - 2.0943951 = 0$ Analytical derivative: $f'(E) = 1 - 0.600 \cos E$ Since $e = 0.600 < 1$, $f'(E) \geq 1 - 0.600 = 0.400 > 0$ for all $E$, ensuring monotonic growth and unconditional convergence for Newton-Raphson. Initial estimate: $E_0 = M = 2.0943951 \text{ rad}$.

Step-by-Step Calculation

1. Iteration 0: $E_0 = 2.0943951$ $\sin(E_0) = \frac{\sqrt{3}}{2} \approx 0.8660254, \quad \cos(E_0) = -0.5000000$ $f(E_0) = 2.0943951 - 0.600(0.8660254) - 2.0943951 = -0.5196152$ $f'(E_0) = 1 - 0.600(-0.5) = 1.3000000$ $E_1 = E_0 - \frac{-0.5196152}{1.3000000} = 2.0943951 + 0.3997040 = 2.4940991$
2. Iteration 1: $E_1 = 2.4940991$ $\sin(E_1) \approx 0.6033702, \quad \cos(E_1) \approx -0.7974609$ $f(E_1) = 2.4940991 - 0.600(0.6033702) - 2.0943951 = +0.0376819$ $f'(E_1) = 1 - 0.600(-0.7974609) = 1.4784765$ $E_2 = 2.4940991 - \frac{0.0376819}{1.4784765} = 2.4940991 - 0.0254870 = 2.4686121$
3. Iteration 2: $E_2 = 2.4686121$ $f(E_2) \approx +0.000185, \quad f'(E_2) \approx 1.468763$ $E_3 = 2.4686121 - \frac{0.000185}{1.468763} = 2.4684860$ Residual is now $|f(E_3)| < 10^{-8}$.

Reading the Result

The eccentric anomaly converges quadratically to $E = 2.468486 \text{ rad}$ ($141.43^\circ$). Orbital radius is $r = a(1 - e \cos E) = a(1 - 0.6(-0.7818)) = 1.469 a$.

Common Error Avoided

Using pure fixed-point iteration $E_{n+1} = M + e \sin E_n$ without relaxation; while contractive for small $e$, it converges slowly linearly compared to the quadratic rate of Newton-Raphson.

Inverse Example: Recovering Intersite Energy Barrier from Arrhenius Inversion

Physical Situation

In materials physics, ionic hopping across solid electrolyte grain boundaries has a reaction rate $k$ governed by the Arrhenius relation: $k = A e^{-E_a / k_B T}$ An experiment measures reaction rate $k = 1.25 \times 10^3 \text{ s}^{-1}$ at $T = 300 \text{ K}$, with pre-exponential factor $A = 10^{13} \text{ s}^{-1}$.

Given Result

Measured rate: $k = 1.25 \times 10^3 \text{ s}^{-1}$. Attempt frequency: $A = 1.0 \times 10^{13} \text{ s}^{-1}$. Temperature: $T = 300 \text{ K}$.

Analytical Reconstruction

Formulating $f(E_a) = A e^{-E_a / k_B T} - k = 0$. Taking natural logarithms: $\ln k = \ln A - \frac{E_a}{k_B T} \implies E_a = k_B T \ln(A / k)$ Evaluating with $k_B T = 0.02585 \text{ eV}$: $E_a = 0.02585 \text{ eV} \times \ln(10^{13} / 1250) = 0.02585 \times \ln(8.0 \times 10^9) = 0.02585 \times 22.8027 \approx 0.5895 \text{ eV}$

What Physical Question It Can Answer

Does this barrier height indicate fast superionic conduction? An activation energy of $E_a \approx 0.59 \text{ eV}$ falls in the characteristic range for lithium battery solid-state electrolytes.

Boundary Verification

As $T \to 0$, $k \to 0$; as $E_a \to 0$, $k \to A$, matching transition state theory limits.

Common Mistakes in Roots and Iterative Methods

1. Running Newton-Raphson near local extrema or vanishing derivatives

It seems reasonable because Newton-Raphson is celebrated as the fastest general-purpose solver. It fails because if $f'(x_k) \approx 0$, division by near-zero slopes inflates step $\Delta x = f(x_k)/f'(x_k)$ by orders of magnitude, throwing the next estimate thousands of units away or generating inf or nan overflows. It is detected when iterate logs exhibit sudden wild swings in sign and magnitude instead of converging. It is corrected by verifying $|f'(x_k)| > \epsilon_{\text{min}}$ before dividing, or employing damped Newton schemes (Levenberg-Marquardt). Mini-example: solving $f(x) = x^3 - x = 0$ with seed $x_0 = 1/\sqrt{3}$ fails immediately because $f'(1/\sqrt{3}) = 3(1/3) - 1 = 0$, dividing by zero.

2. Using single-criterion termination based solely on residual $|f(x)|$

It seems reasonable because $f(x) \approx 0$ intuitively defines a root. It fails because for shallow or flat functions, $|f(x)|$ can easily drop below $10^{-6}$ while coordinate $x$ is still far from the true physical root. It is detected when tightening tolerances causes $x^*$ to shift noticeably despite prior small residuals. It is corrected by requiring simultaneous dual satisfaction: $|f(x_k)| < \epsilon_{\text{res}}$ and $|x_{k+1} - x_k| < \epsilon_{\text{step}} (1 + |x_k|)$. Mini-example: for $f(x) = (x - 100)^{10}$, at $x = 99.5$ the residual is $(0.5)^{10} \approx 0.00097$; a naive residual test accepts $x = 99.5$, missing the root at $100.0$ by half a unit.

3. Blindly accepting the first root found in multi-root physical systems

It seems reasonable because the solver exited cleanly and returned a valid real number. It fails because nonlinear physical laws (van der Waals equations, optical waveguide modes, mechanical equilibria) frequently feature multiple real roots, but only one corresponds to the stable physical state. It is detected when thermodynamic simulations yield negative heat capacities or molar volumes smaller than molecular core volumes. It is corrected by bounding search domains physically and testing thermodynamic stability (concavity of free energy). Mini-example: the van der Waals equation for steam at $100^\circ\text{C}$ and 1 atm has roots for liquid and vapor volumes; initializing near zero finds liquid when steam was intended.

4. Applying Bisection to tangential roots (even multiplicity)

It seems reasonable because Bisection is considered completely fail-safe. It fails because Bisection strictly demands Bolzano bracket satisfaction $f(a) \cdot f(b) < 0$. If a root touches zero without crossing (such as $f(x) = (x - c)^2$), the function maintains the same sign on both sides and no valid bracket can be constructed. It is detected when bracketing routines fail to find opposite signs across the search range. It is corrected by recasting root finding as an optimization problem minimizing $f(x)$ via Golden Section search or solving $f'(x) = 0$. Mini-example: finding the tangential grazing contact of a beam on a parabolic curve $y = x^2$ cannot be handled by Bisection.

5. Omitting iteration safety caps on open methods

It seems reasonable to assume open methods must converge if the model is continuous. It fails because inflection points can trap Newton-Raphson in periodic 2-cycles or chaotic attractors, cycling endlessly between coordinates without converging. It is detected when code hangs in an infinite while loop pegging a CPU core. It is corrected by enforcing a strict iteration cap ($k_{\text{max}} \approx 50 - 100$) triggering a convergence exception upon exhaustion. Mini-example: for $f(x) = x^3 - 5x$, seeding at $x_0 = 1.0$ yields $x_1 = -1.0$, which then returns $x_2 = +1.0$, cycling indefinitely between $+1$ and $-1$.

Quick detection rule

• Is $|f'(x)|$ approaching zero? Halt Newton-Raphson: switch to Bisection or Brent.
• Before returning a root: verify that $f(x^*) \approx 0$ and $x^*$ satisfies physical bounds ($T > 0, V > b$).

Physical Applications of Roots and Iterative Methods

1. Orbital astrophysics: solving Kepler's equation for exoplanet orbits

Appears in observational astronomy when modeling stellar radial velocity Doppler wobbles induced by massive eccentric exoplanets (monitored by spectrographs like ESPRESSO or HARPS). It is calculated by solving Kepler's transcendental equation at each observation epoch $t_i$: $M(t_i) = E(t_i) - e \sin E(t_i)$, where $M(t)$ is known mean anomaly and $E$ is unknown eccentric anomaly. The result yields true orbital anomaly $\nu(t) = 2 \arctan\left( \sqrt{\frac{1+e}{1-e}} \tan\frac{E}{2} \right)$, determining stellar radial velocity $v_r(t) = K [\cos(\nu + \omega) + e \cos\omega]$. The critical danger is applying naive seeds for highly eccentric exoplanets ($e > 0.95$, such as cometary or eccentric Jupiters); for extreme orbits, $E_0 = M$ lies outside the basin of attraction, causing Newton-Raphson to diverge or cycle chaotically. As a final sanity check, Danby's stabilized seed $E_0 = M + 0.85 e \operatorname{sgn}(\sin M)$ is adopted and derivative bounds $1 - e \cos E \geq 1 - e > 0$ are verified.

2. Chemical thermodynamics: vapor-liquid equilibria and phase molar volumes

Appears in chemical process engineering when sizing high-pressure flash drums and distillation columns using cubic equations of state (van der Waals, Peng-Robinson). It is calculated by computing the three real roots of the cubic molar volume polynomial $V_m$ at subcritical temperature and pressure ($T < T_c, P < P_c$). The result identifies the largest root $V_{\text{vapor}}$ (saturated vapor volume), the smallest root $V_{\text{liquid}}$ (condensed liquid volume), and discards the intermediate root as unphysical with negative compressibility ($(\partial P / \partial V)_T > 0$). The critical danger is seeding Newton-Raphson with ideal gas estimates when seeking liquid phase densities; the solver inevitably slides into the vapor root, miscalculating density by a factor of 100 and crashing pump head designs. As a final sanity check, fugacity coefficients $\phi$ are equated across roots to verify chemical potential equality $\mu_{\text{liq}} = \mu_{\text{vap}}$ at phase saturation.

3. Optics and semiconductors: bound energy states in quantum wells

Appears in semiconductor laser engineering (GaAs/AlGaAs heterostructures) for telecommunication optical transmitters. It is calculated by solving transcendental boundary-matching dispersion relations for finite 1D quantum wells of width $2a$ and barrier height $V_0$: $\sqrt{z_0^2 - z^2} = z \tan(z) \quad (\text{even states}), \qquad -\sqrt{z_0^2 - z^2} = z \cot(z) \quad (\text{odd states})$ where $z = \frac{a}{\hbar} \sqrt{2m E}$ and $z_0 = \frac{a}{\hbar} \sqrt{2m V_0}$. The result establishes quantized bound electronic energy levels $E_n$, tuning laser emission wavelength ($\lambda = h c / \Delta E$). The critical danger is evaluating Newton-Raphson steps near tangent poles at $z = \pi/2, 3\pi/2, \dots$, where vertical asymptotes cause explosive divergence. As a final sanity check, bracketing windows $[n\pi, n\pi + \pi/2]$ between poles are isolated and processed via Bisection or Brent's method.

4. Structural mechanics: catenary sag in high-voltage transmission lines

Appears in civil and electrical engineering when calculating mechanical sag and conductor tension across transmission towers separated by span $L$. It is calculated by solving the transcendental catenary equation: $\frac{2 a}{L} \sinh\left( \frac{L}{2 a} \right) = \frac{S}{L}$, where $S$ is wire arc length and $a = H / w$ is the catenary scaling parameter ($H$ is horizontal tension, $w$ is linear weight). The result determines parameter $a$ and maximum support tension $T_{\text{max}} = H + w y_{\text{sag}}$, verifying that wire mechanical stress remains safely below yield limits under winter ice and wind loads. The critical danger is relying on parabolic approximations ($y \approx \frac{x^2}{2a}$) on long spans; parabolas underestimate tower mechanical tension by over $15\%$, risking structural failure in storms. As a final sanity check, residual derivative monotonicity is checked and support reaction loads are tested against structural safety codes.

5. Fluid mechanics: Colebrook-White friction factor in turbulent pipes

Appears in hydraulic engineering when evaluating turbulent frictional pressure loss along pipelines and chemical cooling loops via the implicit Colebrook-White equation: $\frac{1}{\sqrt{f}} = -2.0 \log_{10}\left( \frac{\epsilon / D}{3.7} + \frac{2.51}{\operatorname{Re} \sqrt{f}} \right)$ where $f$ is Darcy friction factor, $\operatorname{Re}$ is Reynolds number, and $\epsilon/D$ is relative roughness. It is calculated via fixed-point iteration or Newton-Raphson in terms of $x = 1/\sqrt{f}$. The result determines the friction factor used to compute pump electrical ratings. The critical danger is initializing with values outside physical bounds ($0.008 \leq f \leq 0.10$), which can destabilize fixed-point iterations if logarithmic arguments approach zero. As a final sanity check, Haaland's explicit formula ($1/\sqrt{f} \approx -1.8 \log_{10}[(\epsilon/3.7D)^{1.11} + 6.9/\operatorname{Re}]$) provides seed $x_0$, enabling Newton-Raphson to reach machine precision in just one or two steps.

Guided Practice: Roots and Iterative Methods

Level 1: recognition

Select whether Bisection, Newton-Raphson, or the Secant method is most suitable for each situation:

1. Finding the root of smooth $f(x) = x^3 - 2x - 5$ on $[2, 3]$ with accessible exact derivative $f'(x) = 3x^2 - 2$.
2. Finding the root of $f(x)$ evaluated via a 5-minute hydrodynamic CFD code lacking derivatives.
3. Guaranteeing absolute unconditional convergence in nuclear reactor safety shutdown triggers.
4. Root finding for $f(x) = x^{1/3}$ starting at $x_0 = 1.0$.

Solution:

1. Newton-Raphson (analytic derivative available, rapid quadratic convergence).
2. Secant method or Brent's inverse quadratic interpolation (no derivative evaluations required, superlinear convergence $\phi \approx 1.618$).
3. Bisection (Bolzano's theorem provides unconditional deterministic convergence).
4. Bisection only: for $f(x) = x^{1/3}$, Newton stepping yields $x_{k+1} = -2 x_k$, causing unbounded geometric divergence ($\{1, -2, 4, -8, \dots\}$).

Level 2: direct application

Apply two full Newton-Raphson iterations to evaluate $\sqrt{7}$ by solving $f(x) = x^2 - 7 = 0$ starting from $x_0 = 3.0$.

Solution:

1. Model and derivative: $f(x) = x^2 - 7, \qquad f'(x) = 2x$ Iteration recurrence (Babylonian square-root formula): $x_{k+1} = \frac{1}{2}\left( x_k + \frac{7}{x_k} \right)$
2. First iteration ($k=0$): $x_1 = \frac{1}{2}\left( 3.0 + \frac{7}{3.0} \right) = \frac{1}{2}(3.0 + 2.3333333) \approx 2.6666667$
3. Second iteration ($k=1$): $x_2 = \frac{1}{2}\left( 2.6666667 + \frac{7}{2.6666667} \right) = \frac{1}{2}(2.6666667 + 2.6250000) \approx 2.6458333$
4. Comparison to exact value: $\sqrt{7} \approx 2.6457513$ Absolute error after two iterations is $|2.6458333 - 2.6457513| = 0.000082$ (relative error under $0.003\%$).

Level 3: tool selection

Solve satellite radiative equilibrium balance: $f(T) = \epsilon \sigma A T^4 + h A (T - T_{\text{space}}) - Q_{\text{solar}} = 0$ Physical temperatures are constrained between $T_{\text{min}} = 100 \text{ K}$ and $T_{\text{max}} = 400 \text{ K}$. What hybrid strategy balances speed with robustness?

Solution: Use Brent's hybrid algorithm:

1. Confirm bracket bounds by verifying $f(100) < 0$ and $f(400) > 0$.
2. Run fast Secant and inverse quadratic interpolation steps.
3. If an interpolation step projects outside current bracket bounds or fails to halve the interval every two iterations, revert dynamically to Bisection.
4. This couples open-method superlinear performance with bracketing stability, preventing runaway steps caused by steep $T^4$ curvature.

Level 4: physical interpretation

Planck's spectral blackbody density reaches its maximum at wavelength $\lambda_{\text{max}}$ satisfying the transcendental relation: $\frac{x e^x}{e^x - 1} - 5 = 0, \qquad x = \frac{h c}{\lambda k_B T}$ Show that this root generates Wien's displacement law $\lambda_{\text{max}} T = b$ and compute $b$ from numerical root $x^* \approx 4.96511423$.

Solution: From the dimensionless coordinate: $x^* = \frac{h c}{\lambda_{\text{max}} k_B T} \implies \lambda_{\text{max}} T = \frac{h c}{k_B x^*}$ Substituting fundamental SI physical constants:

• $h = 6.62607015 \times 10^{-34} \text{ J}\cdot\text{s}$
• $c = 2.99792458 \times 10^8 \text{ m/s}$
• $k_B = 1.380649 \times 10^{-23} \text{ J/K}$
• $x^* = 4.96511423$ $\lambda_{\text{max}} T = \frac{(6.62607015 \times 10^{-34}) \times (2.99792458 \times 10^8)}{(1.380649 \times 10^{-23}) \times 4.96511423} \approx 2.89777 \times 10^{-3} \text{ m}\cdot\text{K}$ Physical interpretation: Wien's empirical displacement constant $b \approx 2.898 \times 10^{-3} \text{ m}\cdot\text{K}$ is a direct mathematical consequence of Planck's spectral peak geometry.

Level 5: intentional mistake

A student solves $f(x) = \arctan(x) = 0$ using Newton-Raphson. The true root is $x^* = 0$. Testing the code with initial seed $x_0 = 1.5$, the iterates blow up: $x_1 \approx -1.69, x_2 \approx 2.32, x_3 \approx -5.11, x_4 \approx 37.9$. Explain this breakdown and calculate critical threshold $x_{\text{crit}}$.

Conceptual blunder: Differentiating gives $f'(x) = \frac{1}{1 + x^2}$. Newton-Raphson updates follow: $x_{k+1} = x_k - (1 + x_k^2) \arctan(x_k)$ For large $|x|$, the derivative decays toward zero ($f' \to 0$). Tangent projections overshoot the origin, landing further away than the starting coordinate ($|x_{k+1}| > |x_k|$).

Critical stability radius: Neutral 2-cycle oscillation occurs when $x_1 = -x_0$: $-x_0 = x_0 - (1 + x_0^2) \arctan(x_0) \implies \frac{2 x_0}{1 + x_0^2} = \arctan(x_0)$ Solving numerically yields critical boundary: $x_{\text{crit}} \approx 1.391745$

• For $|x_0| < 1.391745$: Newton-Raphson converges to zero.
• For $|x_0| = 1.391745$: iterates oscillate endlessly between $\pm 1.3917$.
• For $|x_0| > 1.391745$: iterates diverge to infinity. Starting at $x_0 = 1.5 > x_{\text{crit}}$ lands outside the basin of attraction.