Recall that is a root (or a zero) of a function , if . Newton's method (or the Newton–Raphson Method) is used to find an approximate root of a function by computing successive approximations using the formula

for

for some initial approximation to a root of .

The method is useful when the root of cannot be obtained analytically and is often used an alternative to the Bisection method. However, Newton's method requires that the function , and the initial guess , satisfy the following conditions:

- the derivative of can be calculated;
- the derivative is continuous around the root;
- the initial approximation is sufficiently close to a root of ; and
- the derivative at is nonzero.

In this activity, we will learn how to implement Newton's method in MATLAB.

Use the MATLAB Live Editor to edit and run this Live Script in your browser or your desktop.

- Read each section carefully.
- Run the code within a section before starting to read the next.
- To run the code from each section, click into the code section and then click on the Run Section button (from the toolstrip) or click on the blue stripe on left side of that section as shown below:

Remark: Run the code of each section from top to bottom, otherwise you may get an error.

- The end of a section is indicated with a thin line, like the next one -

Suppose we want to find the root of the function

.

To check out in which range the root is, we first plot in the range as follows:

clear

syms x

fplot(x^3-5,[-3 3])

grid on, xlabel('x'), ylabel('f(x)')

title('Plot of f(x) = x^3-5')

The graph indicates that intersects the x-axis in the range . This means that the root we are looking for is within this interval. So, we can choose as our initial approximation. However, we can choose any other value as long as it is sufficiently close to the root of .

Since , we have that . Newton's iteration scheme takes the form

for

To perform all the calculations, we need to define our own function for named fun:

function y = fun(x)

y = x^3-5;

end

and its derivative named funp:

function yprime = funp(x)

yprime = 3*x^2;

end

Now, we set as our initial approximation:

x1 = 1

Then the first 5 iterations are:

format long

x2 = x1 - fun(x1)/funp(x1)

x3 = x2 - fun(x2)/funp(x2)

x4 = x3 - fun(x3)/funp(x3)

x5 = x4 - fun(x4)/funp(x4)

x6 = x5 - fun(x5)/funp(x5)

Remark: By default, MATLAB shows the output with only 4 decimal places. In this context it is better to see the output with more decimals using the command format.

We can take the value of x6 as a good approximation of the root of . Run this section to see the value of :

fun(x6)

Notice that the output is very small, almost zero. Here MATLAB uses the scientific notation 3.6e-08. It means . If we want more precision, we need to calculate more iterations. However, the previous method is not very efficient, since we will need to rewrite line by line over and over again. This problem can be solved using a for or while loop.

With a for loop we can control the number of iterations we want to calculate. First, we must set our initial value:

x = 1

And now we just need to write the for loop:

for n = 1:5 % For every value of n from 1 to 5

x = x - fun(x)/funp(x) % Perform this calculation

end % Exit "for" loop

fun(x) % Calculate f(x) which is approx 0

That's it! Run this section and compare with the previous outputs.

Now modiy the code. Change the number of iterations to 10, 20 or 50 and then re-run this section. What do you notice? What do you wonder?

We could easily adapt the for loop as a while loop to obtain the same result:

x = 1;

n = 1;

while n <=5

x = x - fun(x)/funp(x)

n = n + 1;

end

fun(x)

However, we would like to add extra conditions to calculate a better approximation of the root. Here is where the while loop turns out to be a very useful tool.

For example, we are happy if the difference between successive iterates is smaller than , that is. We also do not want that the number of iterations gets too large, so we restrict n to . Thus, we write:

x=1; % Set starting value

nmax=25; % Set max number of iterations

eps=1; % Initialise error bound "eps"

n=0; % Initialise "n" (counts iterations)

while eps >= 1e-5 && n <= nmax % Set while-conditions

y = x - fun(x)/funp(x); % Compute next iterate using a different

% variable to calculate the error

eps = abs(y-x); % Compute error

x = y; % Update "x"

n = n + 1; % Update "n" by 1

end % Exit the while loop

n % Display the number of iterations

x % Display the last value of "x"

fun(x) % Calculate f(x) which is approx 0

Remark: Notice we use the scientific notation 1e-5. It means . In other words, 0.00001.

Run this section and compare the results with those obtained with the for loop. Is this a better approximation?

Let's practice what we just learned.

Use Newton's method to approximate the root of the function

in the interval . (Also make sure that f actually has a root in that interval). Define this function and its derivative in the appendix with an appropriate name, similarly to those in section 1.1. Use a for loop and then a while loop to compare the results. Which method is more accurate?

Write your code with a for loop here:

Write your code with a while loop here:

function y = fun(x)

y = x^3-5;

end

function yprime = funp(x)

yprime = 3*x^2;

end

% Write your functions after this comment %

% Write your functions before this comment %