1. Problem: Find n-Dimensional Zero of a Vector-Valued Function of a Vector Argument f(x),

   f(z) = 0,

   where x = [x _i]_n×1, z = [z _i]_n×1 and f(x) = [f _i(x)]_n×1 (Vectors are indicated by boldfaced variables and functions).

2. Taylor Approximation: Letting the (k+1)th iteration be such that x^(k+1) = x^(k)+delta_x^(k) ~ z, so that

   0 = f(x^(k+1)) = f(x^(k)+delta_x^(k)) ~ f(x^(k)) + J^(k)*delta_x^(k),

   assuming the change |delta_x^(k)| is sufficiently small, where

   • x^(k) is the kth Newton approximation,
   • f^(k) = f(x^(k)) is the kth value of the n-dim. function f,
   • J^(k) = ([f _i,j]_n×n)^(k) is the kth Jacobian function,
   • and f _i,j is the jth partial derivative (with respect to x _j) of the ith component of f.
3. Conversion to Linear Algebra Problem:

   J^(k)*delta_x^(k) = - f^(k),

   can be converted to the Augmented Matrix:

   ["A" | "b" ] = [ J^(k) | - f^(k) ].

   which can be solved by Forward Gaussian Elimination to get an approximation for delta_x^(k) which can be used to get the new Newton n-dimensional iterate,

   x^(k+1) = x^(k) + delta_x^(k).

4. Numerical Example:

   * Intersection of a circle of radius 2 and center (0,0) with a exponential function function y = 1 - e^x. So let the component functions be f₁(x₁, x₂) and f₂(x₁, x₂) or alternately f(x,y) and g(x,y) without the subscripts.

   * Vector Function with n=2:

   [ f(x,y) ] = [ x² + y² - 4 ]
   [ g(x,y) ] = [ e^x + y  - 1 ]

   * 2×2 Jacobian J:

   [ J₁₁   J₁₁ ] = [ 2x   2y ]
   [ J₂₁   J₂₂ ] = [ e^x      1]

   * Starting Guess (MATLAB Notation):

   [x; y]⁽⁰⁾ = [0; -2];

   * 0th Interate Interation: Initial Augmented Matrix using Forward Gaussian Elimination and Virtual Full Pivoting (pr = row pivot vector and pc = column pivot vector): The maximal full pivot is found in the element (1,2) of J⁽¹⁾:

   [J | - f || pr | pc ]⁽⁰⁾ =

   [0   -4 | 0 || 1 | 1] -> [ 0   -4.000      | 0 || 1 | 2 ]
   [1   +1 | 2 || 2 | 2] -> [ 1   (-0.25)_m | 2 || 2 | 1 ]

   The multiplier (1/(-4))_m = (-0.25)_m is calculated and is not really need at this step, except to do the elimination: J₂₂ -> J₂₂ - (-0.25)*J₁₂ = 0, since the two zeros in the second row due not cause a change in J₂₁ and "J₂₃" = - g⁽⁰⁾ = +2

   Backsubstitution is trivial due to zeros in the pivot (1st) row, yielding

   delta_x⁽⁰⁾ = [2; 0] => x⁽¹⁾ = x⁽⁰⁾ + delta_x⁽⁰⁾ = [0; -2] + [2; 0] = [2; -2];

   for the first iterate answer.

   * 1st Iterate Elimination: The first pivot is found in element (2,1) of J⁽¹⁾ after substiting x⁽¹⁾ = [x⁽¹⁾; y⁽¹⁾] into the augmented matrix function [J(x,y) | -[f(x,y); g(x,y)] ]:

   [J | - f || pr | pc ]⁽¹⁾ =(4r)

   [4.000   -4.000 | -4.000 || 1 | 1] -> [ (0.5413)_m   -4.541 | -1.624 || 1 | 2 ]
   [7.389   +1.000 | -4.389 || 2 | 2] -> [ 7.389          +1.000 | -4.389 || 2 | 1 ]

   The elimination step used was J_1,j = J_1,j - (0.5413)*J_2,j for j=1, 3, where J_1,3 is the element of the RHS (Right Hand Side) in the first row.

   Backsubstitution yields

   delta_x⁽⁰⁾ =(4r) [-0.6424; 0.3576]; =>

   x⁽²⁾ = x⁽¹⁾ + delta_x⁽¹⁾ =(4r) [2; -2] + [-0.6424; 0.3576] =(4r) [1.358; -1.642];

   for the second iterate answer.

   * 2nd Iterate Elimination: The first pivot is found in element (2,1), again, of J⁽²⁾, after substituting x⁽²⁾ = [x⁽²⁾; y⁽²⁾] into the augmented matrix function:

   [J | - f || pr | pc ]⁽²⁾ =(4r)

   [2.716   -3.284 | -0.5903|| 1 | 1] -> [ (0.6986)_m   -3.983 | +0.3302 || 2 | 1 ]
   [3.888   +1.000 | -1.246 || 2 | 2] -> [ 3.888          +1.000 | -1.246    || 1 | 2 ]

   The elimination step used was J_1,j = J_1,j - (0.6986)*J_2,j for j=1, 3, where J_1,3 is the element of the RHS (Right Hand Side) in the first row.

   Backsubstitution yields

   delta_x⁽¹⁾ =(4r) [-0.2992; -0.0.08290] =>

   x⁽²⁾ = x⁽¹⁾ + delta_x⁽¹⁾ =(4r) [1.358; -1.642] + [-0.2992; -0.0.08290] =(4r) [1.059; -1.725];

   Two norm measures of the closeness of the iterates that might be used for stopping criteria are the infinity (inf or maximal) norms:

   ||delta_x⁽¹⁾|| = ||x⁽²⁾-x⁽¹⁾|| =(4r) 0.2992;

   ||f⁽¹⁾|| = ||[f⁽¹⁾; g⁽¹⁾]|| =(4r) 1.246.

   For more terms, this procedure can be implemented in MATLAB.

Remarks:

• These calculation were done in 4 digit rounding exam precision as indicated by the modified equal sign "=(4r)".

• Also, this webpage was done in pure HTML, so the spacing and symbols will not be as good as typeset graphics webapage or a text.

• See also Heath's short and different Section 5.6.2 on N-Dimensional Newton's Method. In addition, Gerald and Wheatley in the Section 2.12 on Systems of Nonlinear Equations treats the same example somewhat differently.