Introduction

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×nn \times n array of all first-order partial derivatives of your nn functions with respect to your nn 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.

Try It

The demo below solves the two-equation system

f1(x,y)=x2+y2−4=0f_1(x,y) = x^2 + y^2 - 4 = 0

f2(x,y)=x⋅y−1=0f_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−1f-J^{-1}\mathbf{f} from the current point to the next.

<!-- Newton-Raphson for a 2D system — solve f1(x,y)=0 and f2(x,y)=0 simultaneously -->
<div class="hint">Click the canvas to set a starting point, then step through Newton-Raphson iterations. Blue = start, orange = intermediate, green = converged root.</div>
<canvas id="cnv" width="340" height="280" title="Newton-Raphson iteration canvas"></canvas>
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run"  type="button">Run</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
</div>
<div id="status" class="status"></div>
<div id="iter-log" class="iter-log"></div>

/* layout and colour scheme */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* System of two nonlinear equations */
/* f1(x,y) = x^2 + y^2 - 4  (circle r=2) */
/* f2(x,y) = x*y - 1         (hyperbola) */
/* Jacobian of the system */
/* J = [[2x, 2y], [y, x]] */

const W = 340, H = 280;
const SCALE = 55; // pixels per unit
const OX = W / 2, OY = H / 2;

// canvas / math coordinate conversion
const toCanvas = (x, y) => [OX + x * SCALE, OY - y * SCALE];
const fromCanvas = (cx, cy) => [(cx - OX) / SCALE, (OY - cy) / SCALE];

function f1(x, y) { return x * x + y * y - 4; }
function f2(x, y) { return x * y - 1; }

// evaluate the 2x2 Jacobian at (x, y)
function jacobian(x, y) {
  return [[2 * x, 2 * y], [y, x]];
}

// solve the 2x2 linear system J * delta = b via Cramer's rule
function solve2x2(J, b) {
  const det = J[0][0] * J[1][1] - J[0][1] * J[1][0];
  if (Math.abs(det) < 1e-14) return null; // singular or near-singular Jacobian — step cannot be computed
  return [
    (b[0] * J[1][1] - b[1] * J[0][1]) / det,
    (J[0][0] * b[1] - J[1][0] * b[0]) / det
  ];
}

// iteration state
let pts = [];        // list of [x,y] iterates
let running = false;
let intervalId = null;

const cnv = document.getElementById('cnv');
const ctx = cnv.getContext('2d');

// draw axes, grid, and the two curves
function drawBackground() {
  ctx.clearRect(0, 0, W, H);
  // light grid lines
  ctx.strokeStyle = '#e8eef3';
  ctx.lineWidth = 1;
  for (let gx = -3; gx <= 3; gx++) {
    const [cx] = toCanvas(gx, 0);
    ctx.beginPath(); ctx.moveTo(cx, 0); ctx.lineTo(cx, H); ctx.stroke();
  }
  for (let gy = -3; gy <= 3; gy++) {
    const [, cy] = toCanvas(0, gy);
    ctx.beginPath(); ctx.moveTo(0, cy); ctx.lineTo(W, cy); ctx.stroke();
  }
  // coordinate axes
  ctx.strokeStyle = '#adb1b8';
  ctx.lineWidth = 1.2;
  ctx.beginPath(); ctx.moveTo(0, OY); ctx.lineTo(W, OY); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(OX, 0); ctx.lineTo(OX, H); ctx.stroke();

  // blue circle x^2+y^2=4
  ctx.strokeStyle = '#3a86ff55';
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  ctx.arc(OX, OY, 2 * SCALE, 0, 2 * Math.PI);
  ctx.stroke();

  // red hyperbola xy=1
  ctx.strokeStyle = '#ff595e55';
  ctx.lineWidth = 1.5;
  for (let branch = 0; branch < 2; branch++) {
    const sign = branch === 0 ? 1 : -1;
    ctx.beginPath();
    let first = true;
    for (let ix = -3; ix <= 3; ix += 0.02) {
      const iy = sign / ix;
      if (!isFinite(iy) || Math.abs(iy) > 3.5) { first = true; continue; }
      const [cx, cy] = toCanvas(ix, iy);
      if (first) { ctx.moveTo(cx, cy); first = false; }
      else ctx.lineTo(cx, cy);
    }
    ctx.stroke();
  }
}

// draw iteration path with arrows
function drawPath() {
  if (pts.length === 0) return;
  // draw arrows between consecutive iterates
  for (let i = 0; i + 1 < pts.length; i++) {
    const [ax, ay] = toCanvas(pts[i][0], pts[i][1]);
    const [bx, by] = toCanvas(pts[i + 1][0], pts[i + 1][1]);
    ctx.strokeStyle = '#1d3557cc';
    ctx.lineWidth = 1.8;
    ctx.beginPath(); ctx.moveTo(ax, ay); ctx.lineTo(bx, by); ctx.stroke();
    // arrowhead tip
    const angle = Math.atan2(by - ay, bx - ax);
    const hs = 7;
    ctx.fillStyle = '#1d3557cc';
    ctx.beginPath();
    ctx.moveTo(bx, by);
    ctx.lineTo(bx - hs * Math.cos(angle - 0.4), by - hs * Math.sin(angle - 0.4));
    ctx.lineTo(bx - hs * Math.cos(angle + 0.4), by - hs * Math.sin(angle + 0.4));
    ctx.closePath(); ctx.fill();
  }
  // colored dots for each iterate
  for (let i = 0; i < pts.length; i++) {
    const [cx, cy] = toCanvas(pts[i][0], pts[i][1]);
    ctx.fillStyle = i === 0 ? '#3a86ff' : i === pts.length - 1 ? '#0a7d33' : '#ff8800';
    ctx.beginPath(); ctx.arc(cx, cy, 5, 0, 2 * Math.PI); ctx.fill();
  }
}

function redraw() { drawBackground(); drawPath(); }

// perform one Newton-Raphson step; returns true if more steps are possible
function newtonStep() {
  if (pts.length === 0) {
    setStatus('Click the canvas to pick a starting point.', 'info');
    return false;
  }
  const [x, y] = pts[pts.length - 1];
  const F = [f1(x, y), f2(x, y)];
  const res = Math.sqrt(F[0] * F[0] + F[1] * F[1]);
  if (res < 1e-9) {
    setStatus(i18nFmt('Converged! Root ≈ (%x%, %y%)', { x: x.toFixed(6), y: y.toFixed(6) }), 'ok');
    return false;
  }
  const J = jacobian(x, y);
  const dx = solve2x2(J, [-F[0], -F[1]]);
  if (!dx) {
    setStatus('Jacobian is singular here — try a different starting point.', 'bad');
    return false;
  }
  const nx = x + dx[0], ny = y + dx[1];
  if (Math.abs(nx) > 6 || Math.abs(ny) > 6) {
    setStatus('Iterates left the visible region — try a different starting point.', 'bad');
    return false;
  }
  pts.push([nx, ny]);
  const nres = Math.sqrt(f1(nx, ny) ** 2 + f2(nx, ny) ** 2);
  appendLog(i18nFmt('#%n%: (%x%, %y%) res=%r%', { n: pts.length - 1, x: nx.toFixed(5), y: ny.toFixed(5), r: nres.toExponential(2) }));
  setStatus(i18nFmt('Step %n% · residual %r%', { n: pts.length - 1, r: nres.toExponential(2) }), 'info');
  redraw();
  if (nres < 1e-9) {
    setStatus(i18nFmt('Converged! Root ≈ (%x%, %y%)', { x: nx.toFixed(6), y: ny.toFixed(6) }), 'ok');
    return false;
  }
  return true;
}

function setStatus(t, cls) {
  const s = document.getElementById('status');
  s.textContent = t; s.className = 'status ' + cls;
}
function appendLog(line) {
  const el = document.getElementById('iter-log');
  el.innerHTML += line + '<br>';
  el.scrollTop = el.scrollHeight;
}

function doStep() { if (!newtonStep()) stopRun(); }
function startRun() {
  if (running) return;
  running = true;
  document.getElementById('btn-run').textContent = 'Stop';
  intervalId = setInterval(() => { if (!newtonStep()) stopRun(); }, 420);
}
function stopRun() {
  running = false;
  clearInterval(intervalId);
  document.getElementById('btn-run').textContent = 'Run';
}
function doReset() {
  stopRun();
  pts = [];
  document.getElementById('iter-log').innerHTML = '';
  setStatus('Click the canvas to pick a starting point.', 'info');
  redraw();
}

cnv.addEventListener('click', (e) => {
  stopRun();
  const r = cnv.getBoundingClientRect();
  const [x, y] = fromCanvas(e.clientX - r.left, e.clientY - r.top);
  pts = [[x, y]];
  document.getElementById('iter-log').innerHTML = '';
  setStatus(i18nFmt('Starting at (%x%, %y%) — press Step or Run.', { x: x.toFixed(3), y: y.toFixed(3) }), 'info');
  redraw();
});

document.getElementById('btn-step').addEventListener('click', doStep);
document.getElementById('btn-run').addEventListener('click', () => { running ? stopRun() : startRun(); });
document.getElementById('btn-reset').addEventListener('click', doReset);

setStatus('Click the canvas to pick a starting point.', 'info');
redraw();
Result
Console

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=1xy = 1); the root you reach depends entirely on where you start.

The Real Complexity

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

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

∥xk+1−x∗∥≤C∥xk−x∗∥2\|\mathbf{x}^{k+1} - \mathbf{x}^*\| \le C \|\mathbf{x}^k - \mathbf{x}^*\|^2

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

Each step costs O(n3)O(n^3). The iteration xk+1=xk−J(xk)−1f(xk)\mathbf{x}^{k+1} = \mathbf{x}^k - J(\mathbf{x}^k)^{-1}\mathbf{f}(\mathbf{x}^k) is never computed by literally inverting JJ; instead you solve the linear system J(xk) Δx=−f(xk)J(\mathbf{x}^k)\,\Delta\mathbf{x} = -\mathbf{f}(\mathbf{x}^k) via LU decomposition, which costs O(n3)O(n^3) floating-point operations. For large sparse systems the cost drops, but the n3n^3 wall is the reason variants like inexact Newton and quasi-Newton methods (which approximate JJ 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.

Where It Matters

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 nn buses yields 2n2n 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 nn-dimensional parameter space is exponential in nn, Newton-Raphson delivers double-precision accuracy in O(n3log⁡(1/ε))O(n^3 \log(1/\varepsilon)) time.

Conclusion

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.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/newton-raphson-systems/Content licensed under CC BY-NC 4.0.