There are two things that are not entirely correct with your Newton-Raphson technique ... but are definitely correctable! After we fix this, there is no need for the error VPA

you are talking about.

---

Mistake # 1 - Updating an iteration

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

_{(source: mit.edu )}

For the next iteration, you use the previous iteration value. What you are doing is using a loop counter and plugging that into yours f(x)

, which is wrong. This should be the value of the previous iteration .

Mistake # 2 - Mixing character values ​​with numeric values

If you look at how you code your function, you are defining it symbolically, but you are trying to substitute numeric values ​​into your function. This, unfortunately, does not fly with MATLAB. If you really want to substitute values ​​in, you must use subs

. This will replace the actual value for you as a function, x

or whatever independent variable your function uses. Once you do this, your value will still be of type sym

. You have to use this as a double to be able to use it numerically.

---

Also, for efficiency, there is no need to make it an y

array. Just make it the only value that gets updated every iteration. With all of the above, your code is updated to look like this. Keep in mind, I took the derivative of the function before the loop to reduce the amount of computation you need to do. I also split the terms numerator and denominator for Newton-Raphson iterations to clarify the situation and make it more palatable for subs

. Without many words:

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

      

        
        
        
      

    

---

Note that I replaced i

with idx

in the loop. The reason is that it is actually not recommended to use i

or j

as loop indices, since these letters are reserved for representing complex numbers. If you look at this Shai post , you will see that it is actually slower to use these variables as loop indices: Using i and j as variables in Matlab

Anyway, to test it out, let's say 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 is consistent with our knowledge sin(x)

, since the intersection points sin(x)

are in integer multiples pi

. x0 = 2

located next to pi

so it works as we expected.

---

A little bonus for you

In your original code, the values ​​of the root were stored on each iteration in y

. If you really want to do this, you will have to modify your code to look something like this. Remember I pre-highlighted 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 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 following values ​​are the assumptions at each iteration of the Newton-Raphson root. Note that the last element of the array is our last iteration, which is roughly equal pi

.