clear all
format short
set(0, 'defaultaxesfontsize', 12)
set(0, 'defaultlinelinewidth', 1.5)
set(0, 'defaultFunctionLinelinewidth', 1.5)
set(0, 'defaultscattermarkerfacecolor', 'flat')
gcf;
set(gcf, 'Position', [0 0 600 350], 'Theme', 'light')
addpath ../FNC_matlab

Starting MATLAB ...

Executing ...

4.5.Newton for nonlinear systems¶

The rootfinding problem becomes much more difficult when multiple variables and equations are involved.

Particular problems are often posed using scalar variables and equations.

While the equations of Example 4.5.1 are easy to solve by hand, in practice even establishing the existence and uniqueness of solutions for any particular system is typically quite difficult.

4.5.1Linear model¶

To extend rootfinding methods to systems, we will keep to the basic philosophy of constructing easily managed models of the exact function. As usual, the starting point is a linear model. We first need to define what it means to take a derivative of a vector-valued function of a vector variable.

Definition 4.5.2 (Jacobian matrix)

The Jacobian matrix of $\mathbf{f}(\mathbf{x})$ , where $\mathbf{f}$ is $m$ -dimensional and $\mathbf{x}$ is $n$ -dimensional, is the $m\times n$ matrix

\mathbf{J}(\mathbf{x}) = \begin{bmatrix} \rule[2mm]{0pt}{1em}\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n}\\[2mm] \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n}\\[1mm] \vdots & \vdots & & \vdots\\[1mm] \rule[-3mm]{0pt}{1em} \frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n} \end{bmatrix} = \left[ \frac{\partial f_i}{\partial x_j} \right]_{\,i=1,\ldots,m,\, j=1,\ldots,n.}

(4.5.3)

Example 4.5.2

Let

\begin{split} f_1(x_1,x_2,x_3) &= -x_1\cos(x_2) - 1\\ f_2(x_1,x_2,x_3) &= x_1x_2 + x_3\\ f_3(x_1,x_2,x_3) &= e^{-x_3}\sin(x_1+x_2) + x_1^2 - x_2^2. \end{split}

(4.5.4)

Then

\mathbf{J}(x) = \begin{bmatrix} -\cos(x_2) & x_1 \sin(x_2) & 0\\ x_2 & x_1 & 1\\ e^{-x_3}\cos(x_1+x_2)+2x_1 & e^{-x_3}\cos(x_1+x_2)-2x_2 & -e^{-x_3}\sin(x_1+x_2) \end{bmatrix}.

(4.5.5)

If we were to start writing out the terms in (4.5.7), we would begin with

\begin{split} f_1(x_1+h_1,x_2+h_2,x_3+h_3) &= -x_1\cos(x_2)-1 -\cos(x_2)h_1 + x_1\sin(x_2)h_2 + O\bigl(\| \mathbf{h} \|^2\bigr) \\ f_2(x_1+h_1,x_2+h_2,x_3+h_3) &= x_1x_2 + x_3 + x_2h_1 +x_1h_2 + h_3 + O\bigl(\| \mathbf{h} \|^2\bigr), \end{split}

(4.5.6)

and so on.

A multidimensional Taylor series begins with the linear approximation

\mathbf{f}(\mathbf{x}+\mathbf{h}) = \mathbf{f}(\mathbf{x}) + \mathbf{J}(\mathbf{x})\mathbf{h} + O(\| \mathbf{h} \|^2),

(4.5.7)

where $\mathbf{J}$ is the Jacobian matrix of $\mathbf{f}$ at $\mathbf{x}$ , and $\mathbf{h}$ is a small perturbation from $\mathbf{x}$ .

The terms $\mathbf{f}(\mathbf{x})+\mathbf{J}(\mathbf{x})\mathbf{h}$ in (4.5.7) represent the linear part of $\mathbf{f}$ near $\mathbf{x}$ , while the $O(\| \mathbf{h} \|^2)$ term represents the omitted higher-order terms in the Taylor series. If $\mathbf{f}$ is actually linear, i.e., $\mathbf{f}(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b}$ , then the Jacobian matrix is the constant matrix $\mathbf{A}$ and the higher-order terms in (4.5.7) disappear.

4.5.2The multidimensional Newton iteration¶

With a method in hand for constructing a linear model for the vector system $\mathbf{f}(\mathbf{x})$ , we can generalize Newton’s method. Specifically, at a root estimate $\mathbf{x}_k$ , we set $\mathbf{h} = \mathbf{x}-\mathbf{x}_k$ in (4.5.7) and get

\mathbf{f}(\mathbf{x}) \approx \mathbf{q}(\mathbf{x}) = \mathbf{f}(\mathbf{x}_k) + \mathbf{J}(\mathbf{x}_k)(\mathbf{x}-\mathbf{x}_k).

(4.5.8)

We define the next iteration value $\mathbf{x}_{k+1}$ by requiring $\mathbf{q}(\mathbf{x}_{k+1})=\boldsymbol{0}$ ,

\begin{split} \boldsymbol{0} &= \mathbf{f}(\mathbf{x}_k) + \mathbf{J}(\mathbf{x}_k)(\mathbf{x}_{k+1}-\mathbf{x}_k),\\ \end{split}

(4.5.9)

which can be rearranged into

\mathbf{x}_{k+1} = \mathbf{x}_k - \bigl[\mathbf{J}(\mathbf{x}_k)\bigr]^{-1} \mathbf{f}(\mathbf{x}_k).

(4.5.10)

Note that $\mathbf{J}^{-1}\mathbf{f}$ now plays the role that $f/f'$ had in the scalar case; in fact, the two are the same in one dimension. In computational practice, however, we don’t compute matrix inverses.

An extension of our series analysis of the scalar Newton’s method shows that the vector version is also quadratically convergent in any vector norm, under suitable circumstances and when the iteration converges at all.

4.5.3Implementation¶

An implementation of Newton’s method for systems is given in Function 4.5.2. Other than computing the Newton step using backslash and taking vector magnitudes with norm, Function 4.5.2 is virtually identical to the scalar version Function 4.3.1 presented earlier.

Algorithm 4.5.2 (newtonsys)

newtonsys.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function x = newtonsys(f,x1)
% NEWTONSYS   Newton's method for a system of equations.
% Input:
%   f        function that computes residual and Jacobian matrix
%   x1       initial root approximation (n-vector)
% Output       
%   x        array of approximations (one per column, last is best)

    % Stopping parameters.
    funtol = 1000 * eps;    % stop for small || f(x) ||
    xtol = 1000 * eps;      % stop for small || x change ||
    maxiter = 40;           % stop after this many iterations

    x = x1(:);  
    [y, J] = f(x1);
    dx = Inf;
    k = 1;

    while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
        dx = -(J \ y);    % Newton step
        x(:, k+1) = x(:, k) + dx;

        k = k+1;
        [y, J] = f(x(:, k));
    end

    if k==maxiter
        warning('Maximum number of iterations reached.')
    end
end

Example 4.5.3 (Convergence of Newton’s method for systems)

A system of nonlinear equations is defined by its residual and Jacobian.

Tip

This function needs to be defined within a script file or in a file of its own with the .m extension.

f45_nlsystem.m

function [f, J] = nlsystem(x)
    f = zeros(3, 1);   % ensure a column vector output
    f(1) = exp(x(2) - x(1)) - 2;
    f(2) = x(1) * x(2) + x(3);
    f(3) = x(2) * x(3) + x(1)^2 - x(2);
    J(1, :) = [-exp(x(2) - x(1)), exp(x(2) - x(1)), 0];
    J(2, :) = [x(2), x(1), 1];
    J(3, :) = [2 * x(1), x(3)-1, x(2)];
end

Since our system function is defined in an external file here, we need to use @ in order to reference it as a function argument.

nlsystem = @f45_nlsystem;
x1 = [0; 0; 0];    % column vector!
x = newtonsys(nlsystem, x1);
num_iter = size(x, 2)

Let’s compute the residual of the last result in order to check the quality.

r = x(:, end)
back_err = norm(nlsystem(r))

We take the sequence errors in the first component of the solution, applying the log so that we can look at the exponents.

log10( abs(x(1, 1:end-1) - r(1)) )'

This sequence looks to be nearly doubling at each iteration, which is a good sign of quadratic convergence.

4.5.4Exercises¶

Exercise 4.5.4

Two elliptical orbits $(x_1(t),y_1(t))$ and $(x_2(t),y_2(t))$ are described by the equations

\begin{bmatrix} x_1(t) \\ y_1(t) \end{bmatrix} = \begin{bmatrix} -5+10\cos(t) \\ 6\sin(t) \end{bmatrix}, \qquad \begin{bmatrix} x_2(t)\\y_2(t) \end{bmatrix} = \begin{bmatrix} 8\cos(t) \\ 3+12\sin(t) \end{bmatrix},

where $t$ represents time.

(a) ⌨ Make a plot of the two orbits with the following code:

t = linspace(0, 2*pi, 500)';
x1 = -5 + 10*cos(t);   y1 = 6*sin(t);
clf
plot(x1, y1); hold on;  
x2 = 8*cos(t);   y2 = 3 + 12*sin(t);
plot(x2, y2);
axis equal; grid on;

(b) ✍ Write out a $2\times 2$ nonlinear system of equations that describes an intersection of these orbits. (Note: An intersection is not the same as a collision—they don’t have to occupy the same point at the same time.)

(c) ✍ Write out the Jacobian matrix of this nonlinear system.

(d) ⌨ Use Function 4.5.2 to find all of the unique intersections. Add them to the plot from part (a).

Exercise 4.5.5

⌨ Suppose one wants to find the points on the ellipsoid $x^2/25 + y^2/16 + z^2/9 = 1$ that are closest to and farthest from the point $(5,4,3)$ . The method of Lagrange multipliers implies that any such point satisfies

\begin{split} x-5 &= \frac{\lambda x}{25}, \\[1mm] y-4 &= \frac{\lambda y}{16}, \\[1mm] z-3 &= \frac{\lambda z}{9}, \\[1mm] 1 &= \frac{1}{25}x^2 + \frac{1}{16}y^2 + \frac{1}{9}z^2 \end{split}

for an unknown value of $\lambda$ .

(a) Write out this system in the form $\mathbf{f}(\mathbf{u}) = \boldsymbol{0}$ . (Note that the system has four variables to go with the four equations.)

(b) Write out the Jacobian matrix of this system.

(c) Use Function 4.5.2 with different initial guesses to find the two roots of this system. Which is the closest point to $(5,4,3)$ , and which is the farthest?