Root Finding Without the Headache
Newton-Raphson is the first method most people learn, and it is also the first method that will quietly destroy your day if you hand it bad input. It needs a derivative, which is fine if you can write one analytically. It fails if your starting guess lands near a flat region. I once spent an entire afternoon chasing a root that simply did not exist because the function had been evaluated with single-precision floats and the plateau at the bottom looked flat enough to fool the iteration into converging on nonsense. Switching to a bisection bracket first, then handing a refined estimate to Newton-Raphson, solved that problem immediately. A good numerical analyst does not reach for the most elegant algorithm first. They reach for the one that will tell them honestly when it is broken. If you are implementing a solver yourself, always include a fallback strategy and track the residual at each step. Convergence is a necessary but insufficient condition for correctness. A sequence can converge beautifully to the wrong value. Brent's method deserves more attention than it gets. It combines bisection reliability with secant speed and does not require a derivative. Most scientific libraries already implement it, so there is rarely a reason to write your own unless you have a very specific constraint. The tradeoff is that Brent's method is slower than pure Newton-Raphson when the function is well-behaved and your derivative is available. That speed penalty is usually irrelevant compared to the risk of spending hours debugging a diverging Newton iteration.
Linear Systems: Stop Using Gaussian Elimination Blindly
Gaussian elimination works until your matrix is large, ill-conditioned, or sparse. At that point it either fails numerically or wastes enormous time and memory. I once ran a dense LU factorization on a matrix with approximately forty thousand rows and watched it eat nearly eighteen gigabytes of RAM while finishing in about three hours. The same system solved in under four minutes using an iterative Krylov-subspace method after I permuted the matrix to improve its sparsity pattern. Condition number matters far more than most people realize. A system with a condition number above one million is likely too sensitive for double-precision arithmetic to handle reliably without getting help from preconditioning or higher precision. You can check the condition number cheaply with an SVD or a condition estimator built into most linear algebra libraries. If the number is large, the result you get is an educated guess, not a precise answer.
When Iterative Methods Beat Direct Methods
Iterative solvers like GMRES, BiCGSTAB, or CG are designed for sparse systems. CG only works for symmetric positive definite matrices, which covers a meaningful subset of practical problems including many finite-element discretizations. GMRES handles general sparse systems but requires storing basis vectors, which means memory grows with the number of iterations unless you use restarting. Restarting helps with memory but can slow convergence or cause stagnation on difficult problems. Preconditioning is where most of the actual work happens. A poor preconditioner can make an iterative method slower than a direct solver. An effective one can reduce iteration counts from thousands to dozens. Diagonal preconditioning is trivial to implement and often good enough to get started. For more demanding problems, incomplete Cholesky or additive Schwarz preconditioners are reasonable choices depending on your matrix structure.
Get the Full Details
Interpolation: Polynomials Are Not Your Friend at Scale
Polynomial interpolation sounds appealing until you add more points and watch the oscillations explode near the edges. This is the Runge phenomenon and it is not theoretical. It shows up routinely when people try to fit a high-degree polynomial to equispaced data over a wide interval. The error grows without bound as degree increases. Splines solve this by using piecewise low-degree polynomials with continuity constraints at the knots. Cubic splines are the standard choice because they balance smoothness against computational cost. Clamped boundary conditions, where you specify the derivative at the endpoints, usually produce better results than natural splines when you actually know the endpoint behavior. If your data is noisy rather than smooth, interpolation is the wrong tool regardless of whether you use splines or polynomials. Smoothing splines or least-squares fitting with a lower-degree basis are more appropriate. Forcing a curve through every data point adds noise into your model instead of removing it. I have seen people interpolate noisy simulation outputs and then complain that their differentiated results were garbage. They were differentiating noise.
Ordinary Differential Equations: Stiffness Will Bite You
Explicit methods like the classic fourth-order Runge-Kutta are easy to understand and easy to implement. They are also completely unsuitable for stiff systems without impractically small step sizes. Stiffness occurs when your system has widely separated time scales, which is extremely common in chemical kinetics, control systems, and thermal problems. An explicit method might require microsecond steps to remain stable even though the solution is changing slowly on the macro scale. Implicit methods handle stiffness naturally because their stability regions extend much further into the left half plane. The downside is that each step requires solving a system of equations, which costs more per step but allows drastically larger steps. BDF methods, available in solvers like CVODE and scipy.integrate.odeint with the BDF option, are the standard choice for stiff ODEs. Radau methods are another option and tend to be more accurate at the cost of extra function evaluations. There is no universal solver selection rule. I routinely test two or three methods on a new problem before committing to one. The right choice depends on accuracy requirements, available Jacobian information, and how much time you will spend tuning tolerances rather than debugging physics.
Numerical Integration: Quadrature Rules Have Real Limitations
Gaussian quadrature is efficient for smooth integrands on finite intervals. It gives exact results for polynomials up to a certain degree with a minimal number of function evaluations. It fails miserably on integrals with singularities, infinite domains, or discontinuous integrands unless you transform the problem first. Adaptive quadrature routines like QUADPACK handle many of these cases by subdividing intervals and applying higher-order rules where the integrand varies rapidly. Monte Carlo integration sounds impractical until the dimension goes above five or six, at which point deterministic quadrature becomes impossible due to the curse of dimensionality. Monte Carlo error scales as one over the square root of sample count regardless of dimension, which is simultaneously its greatest strength and its worst weakness. You need orders of magnitude more samples than a good quadrature rule on a low-dimensional problem, but you also need far fewer samples than deterministic methods on a high-dimensional one. Variance reduction techniques matter enormously in practice. Control variates, stratified sampling, and importance sampling can reduce required sample counts by factors of ten or more. Antithetic variates are especially easy to implement and often provide noticeable improvement with no additional complexity beyond generating negatively correlated pairs.
Implementation Choices That Matter More Than You Expect
The choice of data structures and memory layout affects numerical performance as much as the choice of algorithm. Row-major versus column-major ordering changes cache behavior dramatically on modern hardware. A matrix-vector multiply that fits in L1 cache can be three to five times faster than one that does not, even when both perform the same number of FLOPs. Fixed-point arithmetic still has legitimate uses in embedded systems and high-frequency trading where deterministic rounding behavior matters more than raw precision. IEEE 754 floating point is the default for a reason, but it introduces rounding errors that accumulate differently depending on the order of operations. Summing a vector from smallest to largest magnitude reduces roundoff error compared to the naive left-to-right approach. Kahan summation improves accuracy further with minimal overhead. Testing numerical code is harder than testing regular software because every output depends on continuous parameters rather than discrete states. Boundary case coverage is insufficient. You need convergence tests, sensitivity analysis, and comparison against analytic solutions or established reference implementations whenever possible. A test suite that only checks exact outputs for fixed inputs will miss the cases that actually matter in production.
When Everything Else Fails
Sometimes the numerical problem is fundamentally unstable. No amount of algorithmic tweaking will fix an ill-posed problem. Reedermensing criteria, Tikhonov regularization, and Bayesian inference frameworks provide ways to stabilize ill-conditioned problems, but they introduce assumptions that you must justify on physical or statistical grounds rather than mathematical convenience. If you regularize without understanding what the regularization term implies about your solution, you are not solving your problem. You are solving a different problem that happens to be easier. The honest approach is to quantify uncertainty at every stage. Report condition numbers. Report residuals. Report convergence diagnostics alongside your final answer. A result without error estimates is not a result, it is a number that someone might mistake for truth.
