MATLAB TUTORIAL for the Second Course. Part 2.5: Newton's Methods

We consider the problem \( {\bf x} = {\bf g} \left( {\bf x} \right) , \) where \( {\bf g}: \, \mathbb{R}^n \mapsto \mathbb{R}^n \) is the function from n-dimensional space into the same space. The fixed-point iteration is defined by selecting a vector \( {\bf x}_0 \in \mathbb{R}^n \) and defining a sequence \( {\bf x}_{k+1} = {\bf g} \left( {\bf x}_k \right) , \quad k=0,1,2,\ldots . \)

A point \( {\bf x} \in \mathbb{R}^n \) is a fixed point of \( {\bf g} \left( {\bf x} \right) \) if and only if \( {\bf x} = {\bf g} \left( {\bf x} \right) . \) A mapping \( {\bf f} \left( {\bf x} \right) \) is a contractive mapping if there is a positive number L less than 1 such that

Theorem: Let \( {\bf g} \left( {\bf x} \right) \) be a continuous function on \( \Omega = \left\{ {\bf x} = (x_1 , x_2 , \ldots x_n ) \, : \ a_i < x_i < b_i \right\} \) such that

• \( {\bf g} \left( \Omega \right) \subset \Omega ; \)
• \( {\bf g} \left( {\bf x} \right) \) is a contraction.

• g has a fixed point \( {\bf x} \in \Omega \) and the sequence \( {\bf x}_{k+1} = {\bf g} \left( {\bf x}_k \right) , \quad k=0,1,2,\ldots ; \quad {\bf x}_0 \in \Omega \) converges to x;
• \( \| {\bf x}_k - {\bf x} \| \le L^k \,\| {\bf x}_0 - {\bf x} \| . \) ■

Example. We apply the fixed point iteration to find the roots of the system of nonlinear equations

If we use the starting value (1.5, 0.5), then

We want to determine why our iterative equations were not suitable for finding the solution near both fixed points (0, 1) and (1.88241, 0.778642). To answer this question, we need to evaluate the Jacobian matrix:

Theorem: Assume that the functions f(x,y) and g(x,y) and their first partial derivatives are continuous in a region that contain the fixed point (p,q). If the starting point is chosen sufficiently close to the fixed point and if

Let \( {\bf g} \, : \, \mathbb{R}^n \mapsto \mathbb{R}^n \) be a continuous function. Suppose we are given a system of equations

In planar case, we have two equations

Example. We apply the Wegstein algorithm to find the roots of the system of nonlinear equations

For systems, it is known the global method called the descent method of which Newton’s iteration is a special case. Newton’s method will be applied once we get close to a root.

Let us consider the system of nonlinear algebraic equations \( {\bf f} \left( {\bf x} \right) =0 , \) where \( {\bf f} \,:\, \mathbb{R}^n \mapsto \mathbb{R}^n , \) and define the scalar multivariable function

• \( \phi ({\bf x}) \ge 0 \) for all \( {\bf x} \in \mathbb{R}^n . \)
• If x is a solution of \( {\bf f} \left( {\bf x} \right) =0 , \) then φ has a local minimum at x.
• At an arbitrary point x₀, the vector \( - \nabla \phi ({\bf x}_0 ) \) is the direction of the most rapid decrease of φ.
• φ has infinitely many descent directions.
• A direction u is descent direction for φ at x if and only if \( {\bf u}^T \nabla \phi ({\bf x}) < 0 . \)
• Special descent directions:
  (a) Steepest descent method: \( {\bf u}_k = -JF^T ({\bf x}_k )\, {\bf f} ({\bf x}_k ) , \)
  (b) Newton’s method: \( {\bf u}_k = -JF^{-1} ({\bf x}_k )\, {\bf f} ({\bf x}_k ) . \) ■

For scalar problems \( f (x) =0 , \) the descent direction is given by \( - \phi' (x) = -2\,f(x) \,f' (x) \) while Newton’s method yields \( -f(x) / f' (x) \) which has the same sign as the steepest descent method.

The question that arises for all descent methods is how far should one go in a given descent direction. To answer this question, there are several techniques and conditions to guarantee convergence to a minimum of φ. For a detailed discussion consult the book by

function [p,iter, err] = newdim(F,JF,p,delta,epsilon,max1)
% Input:    F is the nonlinear system saved in the m-file F.m
%           JF is the Jacobian of F saved as the m-file JF.m
%           p is the initial giuess at the solution
%           delta is the error bound
%           epsilon is the tolerance for F(p)
%           max1 is the number of iterations
% Output:   p is the seidel approximation to the solution
%           iter is the number of iterations required
%           err is the error estimate for p
Y=feval(F,p);
for k=1:max1
    J=feval(JF,p);
    Q=p-(J\Y')';
    Z=feval(F,Q);
    err=norm(Q-p);
    relerr=err/(norm(Q)+epsilon);
    p=Q;
    Y=Z;
    iter=k;
    if (err < delta)|(relerr < delta)|(abs(Y) < epsilon)
        break
    end
end

Complete