GuideNonlinear systemsintermediate

Newton for nonlinear systems

The vector version of Newton's method: the derivative becomes the Jacobian matrix, the quotient becomes one linear system per iteration, and the quadratic order is preserved.

From the derivative to the Jacobian

In scalar Newton one divides by f′(xk)f'(x_k), but there are no quotients between matrices: the role of f′f' is taken by the Jacobian matrix F′(X)F'(X), with entries [F′(X)]ij=∂fi∂xj(X)\bigl[F'(X)\bigr]_{ij}=\frac{\partial f_i}{\partial x_j}(X), and the quotient is replaced by the inverse:

x(k+1)=x(k)−[F′(x(k))]−1F(x(k)),k=0,1,2,…x^{(k+1)}=x^{(k)}-\bigl[F'(x^{(k)})\bigr]^{-1}F(x^{(k)}),\qquad k=0,1,2,\dots
Newton's method for systems.
DerivationDerivation: Newton for systems by linearizationView as its own page →

Linearize and set to zero

  1. The Jacobian matrix collects all first partial derivatives of the coordinate functions:

    F′(X)=[∂f1∂x1⋯∂f1∂xn⋮⋮∂fn∂x1⋯∂fn∂xn]F'(X)=\begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n}\\ \vdots & & \vdots\\ \frac{\partial f_n}{\partial x_1} & \cdots & \frac{\partial f_n}{\partial x_n} \end{bmatrix}

  2. The first-order multivariate Taylor expansion around the current iterate x(k)x^{(k)} is the vector version of the tangent line in the scalar derivation:

    F(X)≈F(x(k))+F′(x(k))(X−x(k))F(X)\approx F\bigl(x^{(k)}\bigr)+F'\bigl(x^{(k)}\bigr)\bigl(X-x^{(k)}\bigr)

  3. We seek the point that makes this linear approximation vanish (the analogue of the tangent crossing the axis). Setting it to zero and solving, with the Jacobian inverse instead of the quotient:

    x(k+1)=x(k)−[F′(x(k))]−1F(x(k))x^{(k+1)}=x^{(k)}-\bigl[F'(x^{(k)})\bigr]^{-1}F\bigl(x^{(k)}\bigr)

Example: building F and its Jacobian

ExampleA 2×2 system

Write the system exey+xcos⁡y=0e^xe^y+x\cos y=0, x+y=1x+y=1 in the form F(X)=0F(X)=0, and compute its Jacobian matrix.

  1. Move everything to the left-hand side:

    F(X)=[exey+xcos⁡yx+y−1]=[00]F(X)=\begin{bmatrix} e^xe^y+x\cos y\\ x+y-1 \end{bmatrix}=\begin{bmatrix}0\\0\end{bmatrix}

  2. Differentiating each component with respect to each variable:

    F′(X)=[exey+cos⁡yexey−xsin⁡y11]F'(X)=\begin{bmatrix} e^xe^y+\cos y & e^xe^y-x\sin y\\ 1 & 1 \end{bmatrix}

This system is solved with Newton in the solved exercise. Note that on the constraint x+y=1x+y=1 one has exey=ex+y=ee^xe^y=e^{x+y}=e: at the solution, xcos⁡(1−x)=−ex\cos(1-x)=-e.

Do not invert the Jacobian

In practice, [F′(x(k))]−1F(x(k))[F'(x^{(k)})]^{-1}F(x^{(k)}) is not computed by inverting the matrix: it is much cheaper to solve at each iteration the linear system

[F′(x(k))] u=F(x(k)),x(k+1)=x(k)−u\bigl[F'(x^{(k)})\bigr]\,u=F(x^{(k)}),\qquad x^{(k+1)}=x^{(k)}-u

Solving a nonlinear system therefore requires solving one linear system per iteration: the two topics are chained. The cost of that solve (n33+n2−n3\frac{n^3}{3}+n^2-\frac{n}{3} products/quotients by Gaussian elimination) dominates the method's total cost and stars in the efficiency analysis.

Quadrature adaptations

The quadrature-based scalar methods adapt directly with the same quotient → linear system substitution. With y(k)=x(k)−[F′(x(k))]−1F(x(k))y^{(k)}=x^{(k)}-[F'(x^{(k)})]^{-1}F(x^{(k)}) (one Newton step as predictor):

x(k+1)=x(k)−2[F′(y(k))+F′(x(k))]−1F(x(k))x^{(k+1)}=x^{(k)}-2\bigl[F'(y^{(k)})+F'(x^{(k)})\bigr]^{-1}F(x^{(k)})
Trapezoid method for systems.
x(k+1)=x(k)−[F′ ⁣(x(k)+y(k)2)]−1F(x(k))x^{(k+1)}=x^{(k)}-\left[F'\!\left(\frac{x^{(k)}+y^{(k)}}{2}\right)\right]^{-1}F(x^{(k)})
Midpoint method for systems.
x(k+1)=x(k)−6[F′(x(k))+4F′ ⁣(x(k)+y(k)2)+F′(y(k))]−1F(x(k))x^{(k+1)}=x^{(k)}-6\left[F'(x^{(k)})+4F'\!\left(\frac{x^{(k)}+y^{(k)}}{2}\right)+F'(y^{(k)})\right]^{-1}F(x^{(k)})
Simpson method for systems.