Prev
Unit 2 · Nonlinear Equations
Systems of Nonlinear Equations
65m100 XP
Next
🔍Nonlinear Equations

Systems of Nonlinear Equations

The Jacobian, Newton's method in n dimensions, finite-difference Jacobians, damping, Broyden and power flow

The Northeast blackout, 9 November 1965
A cold November evening; the transmission lines into southern Ontario are heavily loaded. At 5:16 p.m. a small power swing reaches a protective relay on one of the lines out of the Sir Adam Beck hydroelectric station No. 2, near Niagara Falls. The relay's setting is wrong for that load, and it does what it was set to do: it trips the line. The power that line carried shifts instantly onto its neighbours, which overload and trip in turn, and within five minutes the failure has cascaded through the interconnected grid of Ontario, New York and New England. Within minutes about 30 million people have lost power. Some will be without it for up to 13 hours. Preventing a cascade starts with a question that sounds simple: given the generation and the loads, what voltage will each bus in the network settle at, and how much power will flow in each line? The answer is the solution of the power-flow equations — two nonlinear equations per bus, full of sines and cosines, every one coupled to its neighbours. A large grid has thousands of buses. Two years after the blackout, in November 1967, William Tinney and Clifford Hart of the Bonneville Power Administration published Power Flow Solution by Newton's Method. Combined with ordering techniques that exploit how few lines connect each bus, it made large networks practical to solve in a handful of iterations. Grid operators have used Newton–Raphson power flow ever since.
Write the system as F(x)=0\mathbf F(\mathbf x) = \mathbf 0 with x∈Rn\mathbf x \in \mathbb R^n and F:Rn→Rn\mathbf F: \mathbb R^n \to \mathbb R^n. Newton's 1669 idea carries over unchanged: near a guess, replace F\mathbf F by its linear approximation and solve that instead. The derivative of a vector function is a matrix — the Jacobian.
J(x)=[∂F1/∂x1⋯∂F1/∂xn⋮⋱⋮∂Fn/∂x1⋯∂Fn/∂xn],F(x+Δ)≈F(x)+J(x) ΔJ(\mathbf x) = \begin{bmatrix} \partial F_1/\partial x_1 & \cdots & \partial F_1/\partial x_n \\ \vdots & \ddots & \vdots \\ \partial F_n/\partial x_1 & \cdots & \partial F_n/\partial x_n \end{bmatrix}, \qquad \mathbf F(\mathbf x + \boldsymbol\Delta) \approx \mathbf F(\mathbf x) + J(\mathbf x)\,\boldsymbol\Delta

Row i of the Jacobian is the gradient of equation i. The linearisation is the multivariable Taylor expansion truncated after its linear term.

💡Newton's method for systems
Setting the linearisation to zero gives the Newton step:
  1. Evaluate F(xk)\mathbf F(\mathbf x_k) and J(xk)J(\mathbf x_k).
  2. Solve the linear system J(xk) Δk=−F(xk)J(\mathbf x_k)\,\boldsymbol\Delta_k = -\mathbf F(\mathbf x_k).
  3. Update xk+1=xk+Δk\mathbf x_{k+1} = \mathbf x_k + \boldsymbol\Delta_k.
Never form J−1J^{-1}; solve the system. If J(x∗)J(\mathbf x^*) is nonsingular and the start is close enough, convergence is quadratic, exactly as in one dimension: ∥ek+1∥≤C ∥ek∥2\|\mathbf e_{k+1}\| \le C\,\|\mathbf e_k\|^2.
📐A circle and a hyperbola
Solve x2+y2=4x^2 + y^2 = 4, xy=1xy = 1. With F=(x2+y2−4, xy−1)\mathbf F = (x^2 + y^2 - 4,\ xy - 1), J=[2x2yyx].J = \begin{bmatrix} 2x & 2y \\ y & x \end{bmatrix}. From (2,0.5)(2, 0.5): F=(0.25,0)\mathbf F = (0.25, 0), and solving JΔ=−FJ\boldsymbol\Delta = -\mathbf F gives Δ=(−0.0667,0.0167)\boldsymbol\Delta = (-0.0667, 0.0167). The residual norms then go 2.5×10−1, 4.9×10−3, 3.5×10−6, 2.5×10−12, 4.4×10−162.5\times10^{-1},\ 4.9\times10^{-3},\ 3.5\times10^{-6},\ 2.5\times10^{-12},\ 4.4\times10^{-16} — quadratic. The answer is x=1.931851652578136x = 1.931851652578136; the exact value is 2+3=1.931851652578137\sqrt{2+\sqrt3} = 1.931851652578137. Note that det⁡J=2(x2−y2)\det J = 2(x^2 - y^2) vanishes on the lines y=±xy = \pm x: a start there gives a singular Jacobian, and Newton cannot take a step.
Using a term like nonlinear science is like referring to the bulk of zoology as the study of non-elephant animals.
– Stanislaw Ulam, quoted in Nonlinear system
Where the cost goes. Each Newton step for nn unknowns needs
  • the Jacobian: n2n^2 partial derivatives. If you cannot differentiate F\mathbf F by hand, approximate column jj by a finite difference, (F(x+hej)−F(x))/h(\mathbf F(\mathbf x + h\mathbf e_j) - \mathbf F(\mathbf x))/h — nn extra evaluations of F\mathbf F. Lesson 1-3 tells you to take h≈u ∣xj∣≈10−8∣xj∣h \approx \sqrt{u}\,|x_j| \approx 10^{-8}|x_j|.
  • a linear solve: about 23n3\tfrac23 n^3 operations for a dense matrix by Gaussian elimination, far less if JJ is sparse.
For a 100-unknown problem with a finite-difference Jacobian, five Newton iterations cost 506 evaluations of F\mathbf F. When evaluating F\mathbf F means running a simulation, that is the bill that matters.
💡Broyden's method
Charles Broyden's idea (1965): compute the Jacobian once, then update it cheaply instead of recomputing it. After a step s=xk+1−xk\mathbf s = \mathbf x_{k+1} - \mathbf x_k with change y=Fk+1−Fk\mathbf y = \mathbf F_{k+1} - \mathbf F_k, require the new approximation to satisfy the secant equation Bk+1s=yB_{k+1}\mathbf s = \mathbf y — the multi-dimensional version of the secant slope — while changing BB as little as possible: Bk+1=Bk+(y−Bks) sTsTs.B_{k+1} = B_k + \dfrac{(\mathbf y - B_k\mathbf s)\,\mathbf s^{\mathsf T}}{\mathbf s^{\mathsf T}\mathbf s}. A rank-one update, one evaluation of F\mathbf F per iteration. Convergence becomes superlinear rather than quadratic, so it takes more iterations — but far fewer evaluations. On a 100-unknown benchmark: Newton 5 iterations / 506 evaluations, Broyden 13 iterations / 114 evaluations.
Derivation
What Broyden's update does, and does not, change
It satisfies the secant equation. Multiply the update by the step s\mathbf s; the scalar sTs\mathbf s^{\mathsf T}\mathbf s cancels:
Bk+1s=Bks+(y−Bks)(sTs)sTs=Bks+y−Bks=y.B_{k+1}\mathbf s = B_k\mathbf s + \frac{(\mathbf y - B_k\mathbf s)(\mathbf s^{\mathsf T}\mathbf s)}{\mathbf s^{\mathsf T}\mathbf s} = B_k\mathbf s + \mathbf y - B_k\mathbf s = \mathbf y.
It changes nothing else. For any direction v\mathbf v perpendicular to the step, sTv=0\mathbf s^{\mathsf T}\mathbf v = 0, so
Bk+1v=Bkv+(y−Bks)(sTv)sTs=Bkv.B_{k+1}\mathbf v = B_k\mathbf v + \frac{(\mathbf y - B_k\mathbf s)(\mathbf s^{\mathsf T}\mathbf v)}{\mathbf s^{\mathsf T}\mathbf s} = B_k\mathbf v.
The one step we just took tells us how F\mathbf F changes along s\mathbf s, and nothing about other directions. Broyden's update uses exactly that information and keeps every other direction as it was — the smallest change consistent with the new evidence. In one dimension there is only one direction, and the update reduces to the secant slope of lesson 2-2.
Making Newton global. The full Newton step can overshoot just as it did in one dimension. The standard remedy is a line search: compute the Newton direction Δ\boldsymbol\Delta, then take xk+1=xk+αΔ\mathbf x_{k+1} = \mathbf x_k + \alpha\boldsymbol\Delta with α=1,12,14,…\alpha = 1, \tfrac12, \tfrac14, \dots until the residual ∥F∥\|\mathbf F\| actually decreases. Near the solution α=1\alpha = 1 is always accepted and quadratic convergence returns. Trust-region methods instead limit the length of the step; MINPACK's hybrd (Powell's hybrid method, 1970) — the default in scipy.optimize.root — blends Newton and steepest-descent directions within a trust region and updates the Jacobian with Broyden's formula.
📐Two-bus power flow
A slack bus holds ∣V1∣=1|V_1| = 1 per unit at angle 0. A line of impedance 0.01+0.05j0.01 + 0.05j per unit feeds bus 2, which draws P=0.8P = 0.8, Q=0.3Q = 0.3 per unit. The unknowns are bus 2's voltage magnitude ∣V2∣|V_2| and angle θ2\theta_2; the two equations say that the complex power injected at bus 2, V2 (YV)2‾V_2\,\overline{(Y V)_2}, equals minus the load. From the flat start ∣V2∣=1|V_2| = 1, θ2=0\theta_2 = 0 — every bus at nominal voltage, the convention since Tinney and Hart — Newton with a finite-difference Jacobian drives the mismatch 8.5×10−1→2.9×10−2→4.7×10−5→1.2×10−10→3.6×10−168.5\times10^{-1} \to 2.9\times10^{-2} \to 4.7\times10^{-5} \to 1.2\times10^{-10} \to 3.6\times10^{-16}. Bus 2 settles at 0.9757 pu, −2.17°. A real grid has the same structure with thousands of buses and a Jacobian that is almost entirely zeros.
Derivation
Where the two power-flow equations come from
Represent each bus voltage by a complex number, V=∣V∣eiθV = |V|e^{i\theta} (a magnitude and a phase angle), and the line by its admittance y=1/z=1/(0.01+0.05i)y = 1/z = 1/(0.01 + 0.05i).
Currents from voltages. The current flowing from bus 2 into the line is y(V2−V1)y(V_2 - V_1) by Ohm's law. Collected for both buses, I=YV\mathbf I = Y\mathbf V with the admittance matrix
Y=[y−y−yy].Y = \begin{bmatrix} y & -y \\ -y & y \end{bmatrix}.
Power from current. The complex power a bus injects into the network is S=V I‾S = V\,\overline{I} (voltage times the complex conjugate of current); its real part is active power PP, its imaginary part reactive power QQ. Bus 2 injects nothing of its own — it absorbs the load — so
V2 (YV)2‾=−(0.8+0.3i).V_2\,\overline{(Y\mathbf V)_2} = -(0.8 + 0.3i).
One complex equation is two real ones: real parts equal and imaginary parts equal. Two unknowns, ∣V2∣|V_2| and θ2\theta_2; two equations — the mismatch function of the lab. In a grid with NN buses, YY is N×NN \times N with a nonzero entry for every line, which is why it is sparse.
💡Following a solution: load stepping and arc-length continuation
Engineering systems usually come as families F(x,λ)=0\mathbf F(\mathbf x, \lambda) = \mathbf 0 with a load parameter λ\lambda — the load on a structure, the demand on a grid. The natural strategy is load stepping: solve at λ\lambda, and use that solution as the starting guess for λ+Δλ\lambda + \Delta\lambda, so Newton always starts close. It breaks down at a limit point, where JJ becomes singular and the solution curve turns back on itself. A shallow two-bar truss loaded at its apex (in normalised units, P=w3−3w2+2wP = w^3 - 3w^2 + 2w) can carry at most P=0.385P = 0.385, at deflection w=0.423w = 0.423. Stepping the load by 0.05, Newton converges in 5–6 iterations up to P=0.35P = 0.35; at P=0.40P = 0.40 it fails to converge in 50 — there is no nearby solution — and at P=0.45P = 0.45 it lands on w=2.18w = 2.18: the truss has snapped through to an inverted shape. Arc-length methods (Riks, 1979; Crisfield, 1981) make λ\lambda an unknown too and add one equation fixing the length of each step along the curve. The enlarged Jacobian stays nonsingular at the limit point, and the method follows the whole path: up to 0.385, down to −0.385-0.385 at w=1.577w = 1.577 — where the load would have to pull the other way to hold the truss — and up again. The same idea traces the 'nose curve' of a power grid to its voltage-collapse point (continuation power flow) and finds buckling and collapse loads in nonlinear finite element analysis.
Derivation
Where the truss's limit points are
For load stepping, the equation at load PP is F(w)=w3−3w2+2w−P=0F(w) = w^3 - 3w^2 + 2w - P = 0, and Newton needs its derivative
F′(w)=3w2−6w+2.F'(w) = 3w^2 - 6w + 2.
The load–deflection curve turns back where dP/dw=0dP/dw = 0 — exactly where F′(w)=0F'(w) = 0 and Newton's step −F/F′-F/F' is undefined. Solve the quadratic:
w=6±36−246=1±13=0.4226  or  1.5774.w = \frac{6 \pm \sqrt{36 - 24}}{6} = 1 \pm \frac{1}{\sqrt3} = 0.4226 \ \text{ or } \ 1.5774.
The loads there are P(0.4226)=+0.3849=2/(33)P(0.4226) = +0.3849 = 2/(3\sqrt3), the maximum the truss can carry, and P(1.5774)=−0.3849P(1.5774) = -0.3849. Between the two the slope dP/dwdP/dw is negative: more deflection with less load, a state no load-controlled test can hold. Any load above 0.3849 has only the solution beyond the second limit point — which is why load stepping jumps to w=2.18w = 2.18, the inverted truss.
🕰Basins of attraction
In 1879 Arthur Cayley asked a deceptively simple question: for Newton's method on z3=1z^3 = 1 in the complex plane, which starting points lead to which of the three roots? He solved the quadratic case and admitted the cubic was harder. It is: the boundaries between the three basins are fractals, with all three colours arbitrarily close to every boundary point. The interactive draws them. They are a vivid reminder that 'close enough to the root' can be a very strange set.
💡Where you will use this
Power Systems: Newton–Raphson load flow, exactly as in the lab, on the whole grid. Rocket Propulsion: chemical-equilibrium codes such as NASA's CEA solve for the combustion products by Newton on a system of element balances. Structures and Geotechnics: nonlinear finite elements assemble a 'tangent stiffness matrix' — a Jacobian — and take Newton steps on it. Spacecraft Attitude Control and GNSS Navigation: the position fix from four pseudoranges is a nonlinear system solved by Newton-type iteration.
📚Further reading
  • J. E. Dennis Jr. and Robert Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations (Prentice-Hall, 1983; SIAM Classics 1996) — Newton's method for systems, line searches and Broyden's update, derived carefully.
  • C. T. Kelley, Solving Nonlinear Equations with Newton's Method (SIAM, 2003) — finite-difference Jacobians, globalisation and Broyden in practice.
  • William Tinney and Clifford Hart, 'Power flow solution by Newton's method', IEEE Transactions on Power Apparatus and Systems PAS-86(11), 1449–1460 (1967) — the paper behind every load-flow program.
↓But there is a problem
Look back at step 2 of Newton's method for systems: solve JΔ=−FJ\boldsymbol\Delta = -\mathbf F. We have treated that as one line of code, np.linalg.solve. For two unknowns it is. For a grid of 10,000 buses, a structural model with a million degrees of freedom or a CFD mesh with a hundred million cells, that one line is almost the entire cost of the simulation — and done naively, it can also be where all the accuracy is lost. How do you solve Ax=bA\mathbf x = \mathbf b — quickly, and correctly? That is the next unit.