The idea: If $f(a)$ and $f(b)$ have opposite signs, then by the Intermediate Value Theorem, there must be a root between $a$ and $b$. Cut the interval in half. Check the sign of $f$ at the midpoint. The root is in whichever half has a sign change. Repeat. Each iteration halves the interval, so after $n$ iterations the error is at most $\frac{b - a}{2^n}$.

Strengths: Guaranteed convergence. No derivatives needed. Works on any continuous function. Weakness: Slow. To get 10 decimal digits from an interval of width 1, you need about 33 iterations.

Note how the number of iterations is computed in advance from the interval width and the desired tolerance. There is no convergence check inside the loop—bisection is so predictable that you know exactly how many steps you need.

The idea: Instead of always taking the midpoint, draw a straight line between $(a, f(a))$ and $(b, f(b))$ and find where it crosses the x-axis. This “false position” is usually a better estimate of the root than the midpoint, especially when the function is nearly linear near the root.

False position retains the bracketing guarantee of bisection (one endpoint always stays fixed), but converges faster for well-behaved functions. The downside: one endpoint can get “stuck,” leading to slow convergence in certain cases. Variants like the Illinois method fix this.

The idea: At your current guess $x_n$, approximate the function by its tangent line. Follow the tangent line to where it crosses the x-axis. That crossing is your next guess. Formally:

Newton’s method is astonishingly fast when it works. From a good starting guess, 5–6 iterations typically gives 15+ digits of accuracy. But it can fail spectacularly: if $f'(x_n)$ is near zero, the tangent line is nearly horizontal and the next guess flies off to infinity. It can also cycle, diverge, or converge to the wrong root.

The implementation includes safeguards: a minimum derivative threshold (min_dfx) to catch near-zero derivatives, and a convergence check based on the step size.

The idea: Standard Newton’s method can overshoot. If the function value at the proposed new point is larger in magnitude than at the current point, the full step was too aggressive. This refined version checks: did the step actually improve things? If not, it halves the step. This simple modification significantly improves robustness without sacrificing much speed.

This is a simple form of line search—a technique that becomes essential in optimization and multidimensional root-finding. More sophisticated versions (Armijo rule, Wolfe conditions) are used in gradient descent and quasi-Newton methods.

The idea: Use row operations to eliminate entries below the diagonal, transforming $A$ into an upper triangular matrix $U$. Then solve $Ux = b'$ by back substitution (starting from the last equation and working up). This is the method taught in every linear algebra course.

Why “naive”? Because it doesn’t pivot. If a diagonal entry (the “pivot”) is zero or very small, the algorithm divides by it and either fails or produces wildly inaccurate results. Real implementations always use pivoting.

The problem with naive elimination: Small pivots amplify rounding errors catastrophically. If $A[k,k]$ is $10^{-15}$, dividing by it magnifies errors by a factor of $10^{15}$.

Scaled partial pivoting: Before eliminating column k, find the row (among remaining rows) that has the largest relative entry in column k. “Relative” means divided by the maximum entry in that row—this prevents large rows from dominating. Swap that row to the pivot position. Use an index array to track swaps (cheaper than physically swapping rows).

The factored storage trick: Instead of processing $b$ during elimination, store the multipliers in the zeroed-out entries below the diagonal. This lets you factor $A$ once and then solve for any number of right-hand sides $b$ without repeating the $O(n^3)$ elimination.

This design is important: factoring $A$ costs $O(n^3)$, but each subsequent solve costs only $O(n^2)$. If you have 100 different right-hand sides for the same matrix, you save 99 full eliminations.

The idea: The $j$-th column of $A^{-1}$ is the solution to $Ax = e_j$, where $e_j$ is the $j$-th column of the identity matrix. Factor $A$ once with gauss(), then call solve() $n$ times with different identity columns.

A note on practice: In almost all applications, you should solve $Ax = b$ directly rather than computing $A^{-1}$ and then multiplying $A^{-1}b$. Computing the inverse is slower, less accurate, and rarely necessary. But there are situations (e.g., computing covariance matrices, certain control theory problems) where the inverse itself is needed.

The idea: Many problems (finite differences, spline interpolation, heat equations) produce matrices where only the diagonal and its immediate neighbors are nonzero. These tridiagonal systems have a special structure that Gaussian elimination can exploit: each row operation only affects one row below, so the entire elimination reduces to a simple forward sweep followed by back substitution.

The Thomas algorithm is one of the most elegant results in numerical linear algebra. A problem that would take millions of operations for a large dense matrix takes only a few thousand for a tridiagonal one.

Same idea as the factored Gaussian elimination above, adapted for the tridiagonal case. The forward elimination (which modifies d and produces multipliers) is done once. Then for each new right-hand side b, only the O(n) substitution phase runs.

Standard pivoting in a tridiagonal system is tricky: swapping rows can create fill-in entries that destroy the tridiagonal structure. This implementation is conservative—it only pivots when the swap won’t introduce new nonzeros (i.e., when the entry above the super-diagonal in the pivot row is already zero, or when we’re at the last row). This gives better numerical stability than no pivoting, while preserving the $O(n)$ complexity.

This handles a specific structure that arises in certain finite-difference discretizations and boundary value problems. The matrix has entries below the diagonal in the top half, entries above the diagonal in the bottom half, and a full three-entry row in the middle. The key insight is that the top half can be solved forward (top to bottom), the bottom half can be solved backward (bottom to top), and then the middle row uses both results.

The idea: Rearrange the $i$-th equation to solve for $x_i$:

The key difference from Jacobi: When computing $x_i^{(\text{new})}$, use the already-updated values $x_1^{(\text{new})}, \ldots, x_{i-1}^{(\text{new})}$ for the earlier components. This “immediate use” of new information typically accelerates convergence by a factor of about 2.

The tradeoff: Gauss-Seidel is inherently sequential (each update depends on the previous one), so it cannot be parallelized as easily as Jacobi.

Look carefully at the difference: Jacobi uses x0[j] (old values), Gauss-Seidel uses x[j] (current values, which includes already-updated entries). One character difference, significant performance impact.

The idea: Gauss-Seidel computes a “natural” update. SOR takes a weighted average between the old value and the Gauss-Seidel update, controlled by a relaxation parameter $\omega$:

Finding the optimal $\omega$ is problem-dependent. For some structured matrices (e.g., from discretizing Laplace’s equation), the optimal $\omega$ can be computed analytically and reduces the iteration count from $O(n)$ to $O(\sqrt{n})$. For a $1000 \times 1000$ system, that’s the difference between 1000 and 32 iterations.

The idea: At each step, approximate $x(t + h)$ by following the tangent line at $x(t)$:

Euler’s method is rarely used in practice because it’s inaccurate. But it’s the conceptual foundation for all higher-order methods. Understanding Euler is understanding the basic idea of numerical integration of ODEs.

The idea: Euler’s method uses the first term of the Taylor expansion of $x(t + h)$. The Taylor method uses more terms:

The Taylor method is powerful but impractical for most real-world use because computing higher-order derivatives is tedious and problem-specific. The Runge-Kutta methods achieve the same accuracy without needing any derivatives beyond the first.

The idea: Achieve 4th-order accuracy (matching the Taylor method of order 4) without computing any higher derivatives. Instead, evaluate $f$ at four carefully chosen points within each step and take a weighted average:

Why these weights? The coefficients $(1, 2, 2, 1)/6$ come from Simpson’s rule for numerical integration. $k_1$ evaluates the slope at the left endpoint, $k_2$ and $k_3$ evaluate at the midpoint (using progressively better estimates), and $k_4$ evaluates at the right endpoint. The weighting gives the midpoint estimates double importance, which is characteristic of Simpson’s rule.

In practice: RK4 is the default choice for non-stiff ODEs. It’s simple, accurate, and reliable. For stiff equations (where solutions change on vastly different timescales), implicit methods like backward Euler or BDF methods are needed instead.