Newton-Raphson for Systems

You have one equation and one unknown: draw the tangent, follow it to zero, repeat. That is the original Newton-Raphson method — a simple geometric trick that converges blindingly fast when it works.

But most engineering and scientific problems do not arrive as a single equation. A power-flow network, a chemical reaction at equilibrium, or the kinematics of a robot arm may involve dozens of nonlinear equations in dozens of unknowns, all tangled together. No single tangent line exists here — yet Newton's insight can be rescued.

The key is the Jacobian matrix: the $n \times n$ array of all first-order partial derivatives of your $n$ functions with respect to your $n$ unknowns. At each iteration, Newton-Raphson for systems asks — if we pretend everything is linear right here, where would the system be zero? — and then steps there. The answer requires solving a linear system, not just dividing by a number, but the reward is the same dazzling quadratic convergence: the number of correct decimal digits roughly doubles every step.

Related reading: Newton-Raphson for a single equation gives the one-dimensional foundation this article builds on.

The demo below solves the two-equation system

$$f_1(x,y) = x^2 + y^2 - 4 = 0$$

$$f_2(x,y) = x \cdot y - 1 = 0$$

Click anywhere on the canvas to choose a starting point, then press Step to advance one Newton-Raphson iteration or Run to converge automatically. Each arrow shows the jump $-J^{-1}\mathbf{f}$ from the current point to the next.

<!-- {{c_intro}} -->
<div class="hint">{{hint_text}}</div>
<canvas id="cnv" width="340" height="280" title="{{canvas_title}}"></canvas>
<div class="controls">
  <button id="btn-step" type="button">{{btn_step}}</button>
  <button id="btn-run"  type="button">{{btn_run}}</button>
  <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
</div>
<div id="status" class="status"></div>
<div id="iter-log" class="iter-log"></div>

/* {{c_layout}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; background: #fff; }
.hint { font-size: .85rem; color: #444; margin-bottom: .5rem; line-height: 1.5; }
canvas { display: block; border: 1px solid #cdd9e3; border-radius: 8px;
         cursor: crosshair; margin-bottom: .5rem; max-width: 100%; }
.controls { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
.status { font-size: .9rem; font-weight: 600; min-height: 1.3em; }
.status.ok  { color: #0a7d33; }
.status.bad { color: #c92f3c; }
.status.info { color: #1d3557; }
.iter-log { font-size: .78rem; color: #555; margin-top: .3rem; max-height: 80px;
            overflow-y: auto; line-height: 1.6; }

// Code not found

Notice how the steps grow short very quickly — that is quadratic convergence at work. The system has four roots (the intersections of a circle of radius 2 and a hyperbola $xy = 1$); the root you reach depends entirely on where you start.

Why is the method so fast — and what can go wrong?

The quadratic convergence theorem. If $\mathbf{f} : \mathbb{R}^n \to \mathbb{R}^n$ is twice continuously differentiable and its Jacobian $J(\mathbf{x}^*)$ at the true root $\mathbf{x}^*$ is non-singular, then there exists a ball around $\mathbf{x}^*$ such that any starting point inside it satisfies

$$\|\mathbf{x}^{k+1} - \mathbf{x}^*\| \le C \|\mathbf{x}^k - \mathbf{x}^*\|^2$$

for some constant $C > 0$. Start with an error of $10^{-2}$ and two steps later the error is near $10^{-8}$ — seven extra digits for free.

Each step costs $O(n^3)$. The iteration $\mathbf{x}^{k+1} = \mathbf{x}^k - J(\mathbf{x}^k)^{-1}\mathbf{f}(\mathbf{x}^k)$ is never computed by literally inverting $J$; instead you solve the linear system $J(\mathbf{x}^k)\,\Delta\mathbf{x} = -\mathbf{f}(\mathbf{x}^k)$ via LU decomposition, which costs $O(n^3)$ floating-point operations. For large sparse systems the cost drops, but the $n^3$ wall is the reason variants like inexact Newton and quasi-Newton methods (which approximate $J$ cheaply) exist.

Failure modes.

• Singular or near-singular Jacobian: the linear system has no unique solution; the method diverges or oscillates.
• Poor starting point: far from any root, the linearization is a bad approximation and the iterates can wander.
• Multiple roots: the basin of attraction for each root has a fractal boundary — small changes in the start can land you at a completely different root, or even send the iterates to infinity.

These pitfalls echo the challenges of non-convex optimization: global guarantees are hard; local convergence is excellent.

Anywhere a model produces a set of nonlinear equations that must all be zero at once, Newton-Raphson for systems is the first tool engineers reach for:

• Power-flow analysis: an electric grid with $n$ buses yields $2n$ nonlinear equations in voltage magnitudes and angles. Newton-Raphson solves them in a handful of iterations — fast enough for real-time control.
• Chemical equilibrium: the concentrations of species in a reactor satisfy a system of nonlinear conservation and equilibrium equations. A Newton-Raphson solve finds the steady state without enumerating all possibilities.
• Robot kinematics: computing the joint angles that place a robot end-effector at a target position (inverse kinematics) is a classic nonlinear system. Newton-Raphson drives the residual to zero in real time on modern hardware.
• Computer graphics: implicit surface rendering, physics simulation, and cloth dynamics all reduce to nonlinear systems that Newton-Raphson handles at each time step.
• Financial models: option pricing under stochastic volatility (e.g., Heston model) calibration requires fitting several parameters simultaneously by solving a nonlinear system.

The speed advantage is decisive: where a simple grid search over an $n$-dimensional parameter space is exponential in $n$, Newton-Raphson delivers double-precision accuracy in $O(n^3 \log(1/\varepsilon))$ time.

Newton-Raphson for systems distills a profound engineering philosophy: linearize, solve, repeat. The Jacobian captures how each equation responds to each unknown locally; a single linear solve finds the best correction; and the quadratic convergence guarantee means you almost never need more than ten iterations to reach machine precision — regardless of how many equations you have.

The method is not magic. A bad starting point, a singular Jacobian, or a highly non-smooth function can all derail it. But paired with a good initial guess (from physics intuition, a coarser model, or a simple fixed-point iteration), it is the closest thing numerical analysis has to a universal solver for non-convex optimization problems that happen to have a zero you can write down.