fbpx

Newton method: Approximation is key

0
Lesson 2: Bisection Method – Numerical Analysis and Methods

Newton Method: What you need to know

Newton method and its MATLAB implementation are detailed in this article. We build an intuitive approach to help you understand newton method and why it operates the way it does.

As in my previous post on the fixed point method and bisection method, the goal here is to find a root of

    \[f(x) = 0\]

.

The Newton method applies a local approximation of the function at any given iteration. Said informally, the function is seen as a bunch of infinitely linear piecewise functions. Yes, this is a strong application of taylor series.

How Newton Method Works

Newton does a very smart trick at every iteration of the method. At each iteration, a variable step size is to be adjusted

    \[x_{k+1} = x_k + \Delta_k\]

where

    \[f(x_k + \Delta_k)=0\]

Sounds pretty simple, all we need to do is find \Delta_k at every every iteration that nulls out f(x_k + \Delta_k), could be simple for some functions (like f(x) = x) but in general this is complicated. In other words, Newton works with a taylor series approximation of f(x) in the neighborhood of x_k around \Delta_k, i.e.

    \[f(x_k + \Delta_k) = f(x_k) + \Delta_k f'(x_k) + \mathcal{O}(\Delta_k^2)\]

Now setting f(x_k + \Delta_k)=0 and neglecting \mathcal{O}(\Delta_k^2) we have that

    \[f(x_k) + \Delta_k f'(x_k)=0\]

or in other words

    \[\Delta_k= - \frac{f(x_k)}{f'(x_k)}\]

. Therefore, we now have a closed form express of newtons method as

    \[x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}\]

MATLAB implementation of Newton method

This is easily achieved on MATLAB. First of all, let’s write a newton function so that we can re-use it over and over again

function [ x, error ] = newton( f, df, x0, Niter )
%
% Input:
%   f - input funtion
%   df - f'(x) which is derivative of f(x)
%   x0 - Initial guess
%   tol - tolerance
%   Niter - Number of iterations
%
% Output:
%   x     - is a vector giving root at each iteration
%   error - is a vector giving  error at each iteration

Now let’s write the newton main function.

function [ x, error ] = newton( f, df, x0, Niter )
%
% Input:
%   f - input funtion
%   df - f'(x) which is derivative of f(x)
%   x0 - Initial guess
%   tol - tolerance
%   Niter - Number of iterations
%
% Output:
%   x     - is a vector giving root at each iteration
%   error - is a vector giving  error at each iteration



x(1) = x0 - (f(x0)/df(x0));
error(1) = abs(x(1)-x0);
k = 2;
while (k <= Niter)
    x(k) = x(k-1) - (f(x(k-1))/df(x(k-1)));
    error(k) = abs(x(k)-x(k-1));
    k = k+1;
end

Now, let’s proceed to implement the main function for testing. As we did with the Bisection method, we will re-use the same previous main function

Niter = 20;

Note that in the above, we have assumed the stopping criterion of max number of iterations instead of min tolerance. Next, we can proceed to defining some functions to test, for example:

Progress .. You’re almost done with this article 🙂
%given functions
f1 = @(x)x^3 - 2*x -5; 
f2 = @(x)exp(-x) - x;
f3 = @(x)x*sin(x) -1;
f4 = @(x)x^3 - 3*x^2 +3*x - 1;

Since the Newton method is gradient-based (i.e. it requires full knowledge of the derivative), so we need to compute the derivatives by hand i.e.

    \[f_1'(x) = 3x^2 - 2\]

    \[f_2'(x) =-e^{-x} - 1\]

    \[f_3'(x) = \sin(x) + x\cos(x)\]

and

    \[f_4'(x) = 3x^2 - 6x +3\]

.
Up till now in the main function, we have defined

    \[f_1(x) = x^3 - 2x - 5\]

    \[f_2(x) = e^{-x} - x\]

    \[f_3(x) = xsin(x) - 1\]

and

    \[f_4(x) = x^3 - 3x^2 + 3x - 1\]

which in MATLAB translates to



We can now solve each function using

%derivative (used for newtons method)
df1 = @(x)3*x^2 - 2;
df2 = @(x)-exp(-x) - 1;
df3 = @(x)sin(x) + x*cos(x);
df4 = @(x)3*x^2 - 6*x +3;
%solve using newton
[x_1_newton,error_1_newton] = newton( f1, df1, 1, Niter );
[x_2_newton,error_2_newton] = newton( f2, df2, 1, Niter );
[x_3_newton,error_3_newton] = newton( f3, df3, 1, Niter );
[x_4_newton,error_4_newton] = newton( f4, df4, 1, Niter);

For comparison sake, we will also check whether we converged properly or not. We can either us the error function outputted by our bisection.m function or we can use MATLAB’s fsolve function as

%solving using a library routine
x_1 = fsolve(f1,1);
x_2 = fsolve(f2,1);
x_3 = fsolve(f3,1);
x_4 = fsolve(f4,1);

Right now, we are ready to plot

figure
subplot(2,2,1)
plot(x_1_newton,'m','Linewidth',1)
hold on
plot(x_1*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = x^3 - 2x - 5','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor

subplot(2,2,2)
plot(x_2_newton,'m','Linewidth',1)
hold on
plot(x_2*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = e^{-x} - x','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor


subplot(2,2,3)
plot(x_3_newton,'m','Linewidth',1)
hold on
plot(x_3*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = x\sin(x) - 1','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor

Our complete main function looks like this

Niter = 20;

%given functions
f1 = @(x)x^3 - 2*x -5; 
f2 = @(x)exp(-x) - x;
f3 = @(x)x*sin(x) -1;
f4 = @(x)x^3 - 3*x^2 +3*x - 1;

%their derivative (used for newtons method)
df1 = @(x)3*x^2 - 2;
df2 = @(x)-exp(-x) - 1;
df3 = @(x)sin(x) + x*cos(x);
df4 = @(x)3*x^2 - 6*x +3;
%solve using newton
[x_1_newton,error_1_newton] = newton( f1, df1, 1, Niter );
[x_2_newton,error_2_newton] = newton( f2, df2, 1, Niter );
[x_3_newton,error_3_newton] = newton( f3, df3, 1, Niter );
[x_4_newton,error_4_newton] = newton( f4, df4, 0.5, Niter);


%solving using a library routine
x_1 = fsolve(f1,1);
x_2 = fsolve(f2,1);
x_3 = fsolve(f3,1);
x_4 = fsolve(f4,1);

figure
subplot(2,2,1)
plot(x_1_newton,'m','Linewidth',1)
hold on
plot(x_1*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = x^3 - 2x - 5','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor

subplot(2,2,2)
plot(x_2_newton,'m','Linewidth',1)
hold on
plot(x_2*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = e^{-x} - x','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor


subplot(2,2,3)
plot(x_3_newton,'m','Linewidth',1)
hold on
plot(x_3*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = x\sin(x) - 1','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor

subplot(2,2,4)
plot(x_4_newton,'m','Linewidth',1)
hold on
plot(x_4*ones(Niter,1),'g--','Linewidth',2)
xlabel('Iteration number n','Interpreter','latex')
ylabel('x','Interpreter','latex')
title('Solving f(x) = x^{3} - 3x^{2} +3x - 1','Interpreter','latex')
legend('Newton','Library Routine','Interpreter','latex')
grid on
grid minor
newton method
newton method

Summary

In this article, we have detailed the newton method along with its conditions for proper operation. Moreover, we sliced down the MATLAB code so that we explain each and every block appearing in our youtube lecture.

Follow my other courses such as convex optimization.

Buy me a cup of coffee using the donate link below 😊☕️

PS: I’m on twitter. I retweet stuff around algorithms, python, MATLAB and mathematical optimization, mostly convex.

PSS: We are so close to 100K subscribers on YouTube. It would mean so much if you could share the channel and subscribe to help us sustain.