Newton Method Implementation in C++ – Numerical Computing in C++ Tutorial

In this C++ numerical computing tutorial, we explain how to implement from scratch a simple Newton method (Newton-Raphson method) for numerically solving nonlinear algebraic equations. The YouTube tutorial is given below.

The Newton method can be used to solve the nonlinear algebraic equations having the following form

(1)  

Motivation for studying the Newton method

• In its more complex form, the Newton method serves as the basis of a number of optimization algorithms.
• For example, to solve a nonlinear model predictive control problem for quadrotors, robots, or UAVs, you would most likely use an algorithm whose idea is similar to the Newton method. The same thing applies to other engineering fields as well as to finance.
• By first learning how to implement a “simple” Newton method, you will have a solid base for understanding how to implement optimization methods.

Basic Idea of the Newton Method

The first step is to define the function . For example, consider the nonlinear equation

(2)  

The left-hand side of this equation is the function . That is,

(3)  

To solve the nonlinear equation

(4)  

We introduce the Newton iteration:

(5)  

which is initialized for selected by the user. In (5) is the first-derivative of .

Consider the example

(6)  

Compute the first derivative of :

(7)  

The Newton iteration is

(8)  

However, we need to answer the following question: When to stop the iteration? There are two criteria for stopping the iterations

1. When the max number of iterations is reached.
2. When the residual tolerance is reached.

The residual can simply be defined as an absolute value of the function . That is, when is relatively small we can stop the Newton iteration. Combining 1. and 2. we can write a while loop for iterating through the Newton iteration.

C++ code for implementing the Newton method

The code given below implements the Newton iteration.

#include<iostream>
// we need cmath for fabs() function - absolute value of a floating point number
#include<cmath>
using namespace std;

// function that solves the nonlinear equation f(x)=0 using the Newton-Rapson method
// x is the initial guess that is improved 
// tolerance is the convergence tolerance for the residual |f(x)| 
// maxIterations is the max number of iterations for stopping the method
// verbose is for printing the solution while solving
void NewtonMethod(double &x, double tol=10e-5, int maxIter=1e2, bool verb=true)
{
    int iter=0;
    // fabs is the floating point absolute value
    double f=x*x+5*x+6; // this is the residual and at the same time the value of f
    while ((fabs(f)>tol) && (iter<maxIter))
    {
        if (verb)
        {
            cout<<"Iteration: "<<iter<<" x: "<<x<<" residual: "<<fabs(f)<<endl;
        }
        x=x-f/(2*x+5);
        f=x*x+5*x+6;
        iter++;
        
    }
    cout<<endl;
}

int main()
{
    double x= -10.0;
    double tolerance=1e-4;
    int iterations=500;
    bool verbose=true;
    double f;

    NewtonMethod(x,tolerance,iterations, verbose);
    f=x*x+5*x+6;
    cout<<"The solution is "<< x << " with the residual value of " <<fabs(f)<<endl;
    return 0;
}