Newton’s Method

Overview

Newton’s method is one of the most popular numerical methods, and is even referred by Burden and Faires as the most powerful method that is used to solve for the equation \(f(x) = 0\). This method originates from the Taylor’s series expansion of the function \(f(x)\) about the point \(x1\):

\begin{equation} f(x)=f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right)+\frac{1}{2 !}\left(x-x_{1}\right)^{2} f^{\prime \prime}\left(x_{1}\right)+\cdots \end{equation}

where \(f\), and its first and second order derivatives, \(f'\) and \(f''\) are calculated at \(x_1\). If we take the first two terms of the Taylor’s series expansion we have:

\begin{equation} f(x) \approx f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right). \end{equation}

We then set previous expression to zero (i.e \(f(x) = 0\) ) to find the root of the equation which gives us:

\begin{equation} f\left(x_{1}\right)+\left(x-x_{1}\right) f^{\prime}\left(x_{1}\right)=0. \end{equation}

Rearranging the previous expression we obtain the next approximation to the root, giving us:

\begin{equation} x=x_{2}=x_{1}-\frac{f\left(x_{1}\right)}{f^{\prime}\left(x_{1}\right)} \end{equation}

Thus generalizing previous expression we obtain Newton’s iterative method:

\begin{equation} x_{i}=x_{i-1}-\frac{f\left(x_{i-1}\right)}{f^{\prime}\left(x_{i-1}\right)}, i \in \mathbb{N} \end{equation}

where \(x_{i} \rightarrow \bar{x}\) (as \(i \rightarrow \infty\)), and \(x\) is the approximation to a root of the function \(f(x)\).

Note: As the iterations begin to have the same repeated values i.e. as \(x_i = x_{i+1} = \bar{x}\) this is an indication that \(f(x)\) converges to \(\bar{x}\). Thus \(x_i\) is the root of the function \(f(x)\).

The Method

Step 1:

Let \(\mathbf{x}^{(0)}=\left(x_{1}^{(0)}, x_{2}^{(0)}, \ldots, x_{n}^{(0)}\right)\) be a given initial vector.

Step 2:

Calculate \(J\left(\mathbf{x}^{(0)}\right)\) and \(\mathbf{F}\left(\mathbf{x}^{(0)}\right)\).

Step 3:

We now have to calculate the vector \(\mathbf{y}^{(0)}\), where

\begin{equation} \mathbf{y}= \left[ \begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{array} \right] \end{equation}

In order to find \(\mathbf{y}^{(0)}\), we solve the linear system \(J\left(\mathbf{x}^{(0)}\right) \mathbf{y}^{(0)}=-\mathbf{F}\left(\mathbf{x}^{(0)}\right)\), using Gaussian Elimination.

Note: Rearranging the system in Step 3, we get that \(\mathbf{y}^{(0)}=-J\left(\mathbf{x}^{(0)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(0)}\right)\). The significance of this is that, since \(\mathbf{y}^{(0)}=-J\left(\mathbf{x}^{(0)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(0)}\right)\), we can replace \(-J\left(\mathbf{x}^{(0)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(0)}\right)\) in our iterative formula with \(\mathbf{y}^{(0)}\). This result will yield that

\begin{equation} \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}-J\left(\mathbf{x}^{(k-1)}\right)^{-1} \mathbf{F}\left(\mathbf{x}^{(k-1)}\right)=\mathbf{x}^{(k-1)}-\mathbf{y}^{(k-1)} \end{equation}

Step 4:

Once \(\mathbf{y}^{(0)}\) is found, we can now proceed to finish the first iteration by solving for \(\mathbf{x}^{(1)}\). Thus using the result from Step 3, we have that

\begin{equation} \mathbf{x}^{(1)}=\mathbf{x}^{(0)}+\mathbf{y}^{(0)}= \left[\begin{array}{c} x_{1}^{(0)} \\ x_{2}^{(0)} \\ \vdots \\ x_{n}^{(0)} \end{array}\right] + \left[\begin{array}{c} y_{1}^{(0)} \\ y_{2}^{(0)} \\ \vdots \\ y_{n}^{(0)} \end{array}\right] \end{equation}

Step 5:

Once we have calculated \(\mathbf{x}^{(1)}\), we repeat the process again, until \(\mathbf{x}^{(k)}\) converges to \(\bar{x}\). This indicates we have reached the solution to \(\mathbf{F}(\mathbf{x})=\mathbf{0}\), where \(\bar{x}\) is the solution to the system.

Note: When a set of vectors converges, the norm \(\left\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\|=0\). This means that

\begin{equation} \left\|\mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\|=\sqrt{\left(x_{1}^{(k)}-x_{1}^{(k-1)}\right)^{2}+\cdots+\left(x_{n}^{(k)}-x_{n}^{(k-1)}\right)^{2}}=0 \end{equation}

Usage

imagine that we want to solve the following nonlinear system of equations:

\begin{equation}\left\{\begin{array}{l} f(x, y)=x^{2}+y^{2}-5 \\ g(x, y)=y-3 x+5 \end{array}\right.\end{equation}

then the code will look like this:

// example_newton.cpp

#include <iostream>
#include "../src/numerary.hpp" // Numerary library

using namespace std;
using namespace numerary;

/* System to solve */
void system(double *x, double *fv, int n) {
    fv[0] = x[0]*x[0] + x[1]*x[1] - 5;
    fv[1] = x[1] - 3*x[0] + 5;
}

/* The main function */
int main() {

    const int maxiter = 100; // Maximum interations for set method (100 by default)
    const double eps = 1.e-7; // eps value for method (1.e-7 by default)
    double *x = new double[2], *fv = new double[2];
    short int is_solved;

    // Initial point
    x[0] = 1; x[1] = 2;

    is_solved = Numerary::nonlinear_systems_of_equations(system, x, fv, 2, "newton", eps, maxiter);

    if (is_solved == 1) {
        cout << "Solved!" << endl;
        cout << "x = " << x[0] << endl;
        cout << "y = " << x[1] << endl;
    } else {
        cout << "Method not allowed!";
    }

    delete[] x;
    delete[] fv;

    return 0;
}