Newton-Raphson Method in Matlab

Question

I am new to matlab and I need to create a function that does n iterations of the Newton-Raphson method with starting approximation x = a. This starting approximation does not count as an interation and another requirement is that a for loop is required. I have looked at other similar questions posted but in my case I do not want to use a while loop.

This is what my inputs are supposed to be:

mynewton(f,a,n) which takes three inputs: 
f: A function handle for a function of x.
a: A real number.
n: A positive integer.

And here is my code so far.

function r=mynewton(f,a,n)
syms x;
z=f(x);
y=a;
for i=1:n    
    y(i+1)=y(i)-(z(i)/diff(z(i)));
end
r=y
end

When I try to call the function I get an error saying:

Error in MuPAD command: DOUBLE cannot convert the input expression into a double    array.
If the input expression contains a symbolic variable, use the VPA function instead.
Error in mynewton (line 6)
y(i+1)=y(i)-(z(i)/diff(z(i)));

The question is how do I use this VPA function? Granted the rest of my code probably isn't 100% correct either but any help that addresses the vpa problem or fixes the other parts of my code would be greatly appreciated.

Thanks!

Answer 1

There are two things that aren't quite correct with your Newton-Raphson technique... but certainly fixable! After we fix this, there is no need for the VPA error that you're speaking of.

---

Error #1 - The iteration update

The first one is the iteration itself. Recall the definition of the Newton-Raphson technique:

_{(source: mit.edu)}

For the next iteration, you use the previous iteration's value. What you're doing is using the loop counter and substituting this into your f(x), which is not correct. It must be the previous iteration's value.

Error #2 - Mixing symbolic values with numeric values

If you take a look at how you're coding your function, you define your function symbolically, yet you are trying to substitute numeric values into your function. This unfortunately doesn't fly with MATLAB. If you actually want to substitute values in, you have to use subs. This will substitute an actual value for you as a function of x or whichever independent variable your function is using. Once you do this, your value is still a sym type. You need to cast this as double in order to be able to use this numerically.

---

Also for efficiency, there is no need to make y an array. Just make this a single value that updates itself at each iteration. With all of this being said, your code is updated to look like this. Mind you, I took the derivative of the function before the loop to decrease the amount of computation you need to take. I also split up the numerator and denominator terms for the Newton-Raphson iterations to make things clear, and to make this more palatable for subs. Without further ado:

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z); %// Edit - Include derivative
y = a; %// Initial root

for idx = 1 : n    
    numZ = subs(z,x,y); %// Numerator - Substitute f(x) for f(y)
    denZ = subs(diffZ,x,y); %// Denominator - Substitute for f'(x) for f'(y)
    y = y - double(numZ)/double(denZ); %// Update - Cast to double to get the numerical value
end
r = y; %// Send to output
end

---

Take note that I replaced i with idx in the loop. The reason why is because it is actually not recommended to use i or j as loop indices as these letters are reserved to represent complex numbers. If you take a look at this post by Shai, you'll see that it's actually slower to use these variables as loop indices: Using i and j as variables in Matlab

In any case, to test this out, suppose our function was y = sin(x), and my initial root was x0 = 2, with 5 iterations, we do:

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

3.1416

This agrees with our knowledge of sin(x), as the intercepts of sin(x) are located at integer multiples of pi. x0 = 2 is located near pi so this does work as we expected.

---

Little bonus for you

Your original code had the values of the root stored at each iteration in y. If you really want to do that, you'll have to modify your code so that it looks something like this. Bear in mind that I pre-allocated y to make things more efficient:

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z);
y = zeros(1,n+1); %// Pre-allocate output array 
y(1) = a; %// First entry is the initial root

for idx = 1 : n    
    numZ = subs(z,x,y(idx)); %// Remember to use PREVIOUS guess for next guess
    denZ = subs(diffZ,x,y(idx));
    y(idx+1) = y(idx) - double(numZ)/double(denZ); %// Place next guess in right spot  
end
r = y; %// Send to output
end

By running this code using the exact same parameters as above, we get:

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

    2.0000    4.1850    2.4679    3.2662    3.1409    3.1416

Each value in r tells you the guess of the root at that particular iteration. The first element of the array is the initial guess (of course). The next values are the guesses at each iteration of the Newton-Raphson root. Take note that the last element of the array is our final iteration, which is roughly equal to pi.

Answer 2

You need to use eval

So z(x)=f(x), then eval(z(1)) to evaluate the function at x = 1

Full code:

function r=mynewton(fun,a,n)
% e.g. for mynewton(@sin, 1, 2)
syms x;
z(x)=fun(x);  % sin(x)
y=a;
for i=1:n    
    y(i+1)=y(i)-(eval(z(i))/eval(diff(z(i))));
end
r=y
end

Not sure about the VPA part of the question, I usually ignore MUPAD errors :P