Skip to article frontmatterSkip to article content

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 base it on the multidimensional Taylor series,

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

where J\mathbf{J} is called the Jacobian matrix of f\mathbf{f} and is defined by

J(x)=[f1x1f1x2f1xnf2x1f2x2f2xnfnx1fnx2fnxn]=[fixj]i,j=1,,n.\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_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \cdots & \frac{\partial f_n}{\partial x_n} \end{bmatrix} = \left[ \frac{\partial f_i}{\partial x_j} \right]_{\,i,j=1,\ldots,n}.

Because of the Jacobian’s role in (4.5.3), we may write J(x)\mathbf{J}(\mathbf{x}) as f(x)\mathbf{f}{\,}'(\mathbf{x}). Like any derivative, it is a function of the independent variable x\mathbf{x}.

The terms f(x)+J(x)h\mathbf{f}(\mathbf{x})+\mathbf{J}(\mathbf{x})\mathbf{h} in (4.5.3) represent the linear part of f\mathbf{f} near x\mathbf{x}. If f\mathbf{f} is actually linear, i.e., f(x)=Axb\mathbf{f}(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b}, then the Jacobian matrix is the constant matrix A\mathbf{A} and the higher-order terms in (4.5.3) disappear.

4.5.2The multidimensional Newton iteration

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

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

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

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

which can be rearranged into

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

Note that J1f\mathbf{J}^{-1}\mathbf{f} now plays the role that f/ff/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.2 presented earlier.

4.5.4Exercises

  1. ✍ Suppose that

    f(x)=[x1x2+x221x1x23+x12x22+1].\mathbf{f}(\mathbf{x}) = \begin{bmatrix} x_1x_2+x_2^2-1 \\[1mm] x_1x_2^3 + x_1^2x_2^2 + 1 \end{bmatrix}.

    Let x1=[2,1]T\mathbf{x}_1=[-2,1]^T. Use Newton’s method to find x2\mathbf{x}_2.

  2. ✍ Suppose that f(x)=Axb\mathbf{f}(\mathbf{x}) = \mathbf{A}\mathbf{x} - \mathbf{b} for a constant n×nn\times n matrix A\mathbf{A} and constant n×1n\times 1 vector b\mathbf{b}. Show that Newton’s method converges to the exact root in one iteration.

  3. Two curves in the (u,v)(u,v) plane are defined implicitly by the equations ulogu+vlogv=0.3u\log u + v \log v = -0.3 and u4+v2=1u^4 + v^2 = 1.

    (a) ✍ Write the intersection of these curves in the form f(x)=0\mathbf{f}(\mathbf{x}) = \boldsymbol{0} for two-dimensional f\mathbf{f} and x\mathbf{x}.

    (b) ✍ Find the Jacobian matrix of f\mathbf{f}.

    (c) ⌨ Use Function 4.5.2 to find an intersection point starting from u=1u=1, v=0.1v=0.1.

    (d) ⌨ Use Function 4.5.2 to find an intersection point starting from u=0.1u=0.1, v=1v=1.

  4. Two elliptical orbits (x1(t),y1(t))(x_1(t),y_1(t)) and (x2(t),y2(t))(x_2(t),y_2(t)) are described by the equations

    [x1(t)y1(t)]=[5+10cos(t)6sin(t)],[x2(t)y2(t)]=[8cos(t)3+12sin(t)],\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 tt represents time.

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

    x1(t) = -5+10*cos(t);   y1(t) = 6*sin(t);
    plot(x1,y1,0,2pi,aspect_ratio=1,legend=false)
    x2(t) = 8*cos(t);   y2(t) = 3 + 12*sin(t);
    plot!(x2,y2,0,2pi)

    (b) ✍ Write out a 2×22\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.

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

    x5=λx25,y4=λy16,z3=λz9,1=125x2+116y2+19z2\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 λ.

    (a) Write out this system in the form f(u)=0\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)(5,4,3), and which is the farthest?

  6. ⌨ Any three noncollinear points in the plane determine a unique circle. Suppose the points are given as (xi,yi)(x_i,y_i) for i=1,2,3i=1,2,3. We can define the circle in terms of its center (a,b)(a,b) and radius rr. Then

    fi(a,b,r)=(axi)2+(byi)2r2f_i(a,b,r) = (a-x_i)^2 + (b-y_i)^2 - r^2

    should be made zero for all i=1,2,3i=1,2,3. This defines a nonlinear system f(v)=0\mathbf{f}(\mathbf{v})=\boldsymbol{0} for v=[a,b,r]\mathbf{v}=[a,b,r].

    Use Function 4.5.2 on this system to find the circle passing through (5,0)(-5,0), (1,3)(1,-3), and (4,2)(4,2). Make a plot that shows you found the correct circle.