Systems of Equations: Newton with the Jacobian

From one equation to many: the shape of the problem

Everything so far in this rung lived on a line: one unknown x, one equation f(x) = 0, a curve crossing an axis. But the real world rarely hands you a single equation. A robot arm's two joint angles must satisfy two equations at once so the hand lands on a target. A circuit's node voltages must jointly satisfy Kirchhoff's laws. A chemical mixture's concentrations must balance several reactions simultaneously. In every case you have n unknowns stacked into a vector x = (x_1, ..., x_n), and n equations stacked into a vector function F(x) = (f_1(x), ..., f_n(x)), and you want the x where all of them vanish together: F(x) = 0. That is the system of nonlinear equations this guide is about.

Notice immediately what we lose. In one dimension, bisection rode on a sign change — f was below the axis here and above it there, so a root was trapped between. In n dimensions there is no axis to be on one side of, no clean notion of 'opposite signs across an interval', and no multidimensional bisection at all. The bracketing methods that gave us iron-clad guarantees simply do not generalize. What does generalize is the idea behind Newton's method: at your current guess, replace the curved function by its best straight-line approximation, solve that easy linear version exactly, and step there. In many dimensions the 'straight-line approximation' becomes a flat plane — a linear map — and solving it becomes solving a linear system.

The Jacobian: a derivative that is a matrix

In one dimension Newton's step came from the tangent line, whose slope was the single number f'(x). With n functions each depending on n variables, 'the slope' is no longer one number — it is a whole grid of partial derivatives. The Jacobian matrix J(x) collects them: its entry in row i, column j is the partial derivative of the i-th function f_i with respect to the j-th variable x_j. Row i is the gradient of f_i, telling you how that one equation responds as you nudge each variable; column j tells you how all the equations respond when you nudge variable x_j alone. For n = 2 it is a 2-by-2 grid of four partials; for n unknowns it is an n-by-n matrix.

The Jacobian plays exactly the role f'(x) played before. Near a point x_n, the function is well approximated by its linearization: F(x) is about F(x_n) + J(x_n) * (x - x_n). This is the multidimensional tangent — a flat plane hugging the curved surface at x_n. To find where this flat approximation hits zero, set the right-hand side to zero and solve for the step s = x - x_n. That gives the defining equation of the Newton step: J(x_n) * s = -F(x_n), a linear system A x = b in disguise, with the Jacobian as A, the negative function value as b, and the step s as the unknown.

one dimension:   x_{n+1} = x_n - f(x_n) / f'(x_n)

n dimensions:    solve   J(x_n) * s = -F(x_n)   for the step s
                 then    x_{n+1} = x_n + s

(division by f'(x_n) becomes: SOLVE a linear system with J(x_n))

Why each step is a linear solve — and why that matters

Here is the beautiful unification. The scalar Newton update divided by f'(x_n); the vector version cannot divide by a matrix, so 'dividing' becomes solving J(x_n) * s = -F(x_n). Suddenly everything you learned in the linear-systems rung comes flooding back into the root-finding rung. You do not compute the matrix inverse — you solve the system, almost always by LU factorization (Gaussian elimination), which factors J = L U once and then finds s with a fast forward-and-back substitution. If the Jacobian barely changes between iterations you can even reuse a stale factorization for several steps, trading a touch of convergence speed for far less work per step.

And the warnings from the linear-systems rung come back too — they are not optional fine print. The Newton step is only as trustworthy as the linear solve that produced it, and that solve inherits the condition number of the Jacobian. If J(x_n) is nearly singular — which happens precisely near a degenerate or multiple root, the multidimensional cousin of f'(x) being near zero — the step blows up or points in a wildly wrong direction, just as the one-dimensional Newton step exploded when f'(x_n) was tiny. Accuracy is conditioning times stability: a stable LU solve still loses about 8 of your roughly 16 double-precision digits when the Jacobian's condition number reaches 10^8, and no clever coding buys those digits back. Every step's answer is an approximation on the floating-point grid; the Jacobian's conditioning decides how good an approximation it can possibly be.

The algorithm, step by step

Stripped to its bones, Newton-for-systems is the scalar loop with 'divide' swapped for 'solve a linear system'. Start from an initial guess, build the function value and the Jacobian there, solve for the step, take it, and repeat until F(x) is small enough. The only genuinely new pieces of machinery are computing the n-by-n Jacobian at each iteration and doing a linear solve each time — both of which you already met in earlier rungs.

1. Choose an initial guess x_0. As in one dimension, Newton is only locally convergent — a good starting point matters enormously, and a bad one can diverge.
2. At the current x_n, evaluate the function vector F(x_n) and assemble the Jacobian matrix J(x_n) of partial derivatives.
3. Solve the linear system J(x_n) * s = -F(x_n) for the step s — by LU factorization (Gaussian elimination), never by forming the inverse.
4. Update x_{n+1} = x_n + s.
5. Check the stopping test: is ||F(x_{n+1})|| below tolerance, and is the step ||s|| small? If yes, stop; if not, go back to step 2. Cap the iterations so a divergent run cannot loop forever.

Getting the Jacobian right is the part people stumble on. Hand-deriving n^2 partial derivatives is error-prone, and a single wrong entry quietly destroys the quadratic convergence without breaking the program — it just limps along linearly, or stalls. So in practice the Jacobian is supplied by automatic differentiation when possible, or approximated by finite differences: nudge each variable by a small h and watch how F changes. But beware — the round-off floor from the floating-point guide applies here too. Too large an h gives truncation error; too small an h drowns the difference in round-off, since you are subtracting two nearly equal vectors. There is a sweet spot near the square root of the machine epsilon, and going smaller does not keep helping.

Speed, failure, and the cheaper cousins

The reward for all this machinery is the same spectacular speed Newton showed in one dimension: near a well-behaved root — a simple root where the Jacobian is nonsingular — convergence is quadratic, roughly doubling the number of correct digits per iteration. A few steps can take you from a rough guess to machine precision. That is why Newton-for-systems is the default engine inside almost every serious solver for nonlinear equations, from circuit simulators to weather models to the implicit time-steppers that stiff differential equations force on us.

But every failure mode from the scalar guide survives the jump to many dimensions, and a few new ones appear. From far away Newton can diverge, overshoot, or wander into a region where the Jacobian goes singular and the step explodes. The honest fixes are the multidimensional versions of the scalar ones: a damped or line-search Newton method, which takes the Newton direction but shortens the step (x_{n+1} = x_n + t*s with 0 < t <= 1) until F actually shrinks, trading a little of the quadratic speed for a much larger basin of starting points that converge. Globalized Newton methods like these are what make the quadratic local rate usable from a realistic, imperfect initial guess. Quadratic convergence is a promise about the neighbourhood of the root, never a promise from anywhere.

There is also a cost question Newton forces on you. Each step needs a fresh Jacobian (n^2 partial derivatives) and a linear solve (about O(n^3) work by Gaussian elimination for a dense system). When the Jacobian is expensive or unavailable, the multidimensional analogue of the secant method takes over: quasi-Newton methods such as Broyden's never compute J at all. They start from a rough Jacobian and update it cheaply using only the function values you are already collecting, much as the secant method replaced f'(x_n) by a slope through two recent points. You give up exact quadratic convergence for something fast and far cheaper per step — often the better trade when each Jacobian costs more than the iterations it saves.

Knowing when to stop, and the whole rung in one frame

Stopping is subtler than in one dimension, because there is no bracket whose width honestly bounds the error. The usual stopping criterion watches two things together: the residual ||F(x_n)|| (how close the equations are to balancing) and the step size ||x_{n+1} - x_n|| (how little the guess is still moving), each measured relative to scale, and stops only when both are small. Watching only the residual can be fooled by an ill-conditioned Jacobian, where a tiny residual hides a not-so-tiny error; watching only the step size can stop early when the method merely slows down. And the floor is the same as everywhere in this ladder: you cannot drive ||F|| below the level set by floating-point round-off and the problem's conditioning, so chasing extra digits past that point just burns iterations.

Step back and the whole rung clicks into one frame. Bisection gave a guarantee but only linear speed and only in one dimension. Fixed-point iteration revealed the contraction condition under which simple iteration converges at all. Newton's method bought quadratic speed by trusting the tangent, and the failure-and-secant guide showed how that trust can break and how to soften or sidestep it. This final guide carries the very same Newton idea into n dimensions, where the tangent slope becomes the Jacobian matrix and the division becomes a linear solve — which is why the conditioning and stability lessons of the linear-algebra rungs are not a separate subject but the beating heart of nonlinear root-finding. Solve A x = b well, and you can solve F(x) = 0 well; the two halves of numerical computing meet exactly here.