Newton Raphson method - Symbolic function

Hello, I made a calculation according to Newton Raphson method, but using syms feature for taking derivative. I believe this is where my problem begins, I only need 4 iterations at total when solved, and parameter k is supposed to be 0,567 as a final value. But my result in terms of iteration is above 10000+ and result is -9 something.

Could you please spot my error/s and correct them?

%%

% finding root of f = (e^-x) - x with Newton Raphson method
clc,clear
  syms x

  f = (exp(-x))-x;  % f = (e^-x) - x
  df = diff(f);
  k=0;              %initial value of x0 (k in my notation) is defined as 0 in the problem
  ep = 100; 
  
  
  t = 10^-8; % result of epsilon value is supposed to be near this value for iterations to stop. 
           % epsilon =  absolute val of 100 * [(xi+1 - xi) / xi+1]
i=0;


      
    while ep > t
  a = vpa(subs(df,x,i));
  b = vpa(subs(f,x,i));
   
    
    k_New = k-a/b;
    ep = 100 * abs(k-k_New)/abs(k_New);
    k = k_New;
    i = i+1;
   
    end

disp(k)

[J. Alex Lee]

why are you bothering to use syms? that residual is so easily differentiable

[Onur Aytan]

I wanted to obtain something universial, applicable for different problems. I might try to make this script even more complicated with additional criterion in the future. Therefore I wanted to use differential concept and N-R method within the same script.

But matlab seems to resist on bothering ppl when we want to use "syms" feature apparently.

[J. Alex Lee]

That's a sensible response. For 1D problems though, you might consider testing your general solution against things like built-in fzero(), and then ask and answer if you gained any benefit with your custom solution.

For real world problems, I doubt that taking derivatives will be your true bottle neck.

I don't think there's anything wrong intrinsically with symbolic toolbox, I'm sure it's useful. I would just be weary of mixing it in with numerical methods, especially if the objective is efficient numerical solutions. Of course best choice comes down to what are your objectives and constraints.

Anyway, good for you for going through this exercise!

[Onur Aytan]

thanks, yeah the greatest benefit is the brain-storming and gaining max experience through such errors. Makes me even happier since I will never forget.

Answer 1

[Alan Stevens]

You have

  a = vpa(subs(df,x,i));
  b = vpa(subs(f,x,i));
   
    
    k_New = k-a/b;    %  but this is k - diff(f)/f  
                      %  you should have k - f/diff(f) or  k - b/a  for Newton-Raphson

[Onur Aytan]

Omg, so true, I must have missed that part while trying different options, cutting and pasting different expressions on my script.

Thanks a lot.

[Onur Aytan]

Though, about the iteration limit of 10^-8.... I still have doubts matlab would perceive that long digit. Right now my matlab is still processing strangely that in reality I only needed 4 iterations when I solve by hand.

[Alan Stevens]

I don't know why your code is so slow, though symbolic calculations are often much slower than numerical ones. You are better off here by using only numerical calculations. You can generalise the derivative as shown in the following

% finding root of f = (e^-x) - x with Newton Raphson method
  f = @(x)(exp(-x))-x;  % f = (e^-x) - x
  h = 10^-10;
  df = @(x) (f(x+h)-f(x))/h;
  k=0;              %initial value of x0 (k in my notation) is defined as 0 in the problem
  ep = 100; 
    
  t = 10^-8; % result of epsilon value is supposed to be near this value for iterations to stop. 
           % epsilon =  absolute val of 100 * [(xi+1 - xi) / xi+1]
i=0;
    
    while ep > t
  a = f(k);
  b = df(k);
   
    k_New = k-a/b;
    ep = 100 * abs(k-k_New)/abs(k_New);
    k = k_New;
    i = i+1;
   
    end

disp(k)

This produces 0.5671 with i = 5 at the end.

[Alan Stevens]

Come to think of it, I suspect your symbolics should have

a = vpa(subs(df,x,k));  %  i.e. k not i
  b = vpa(subs(f,x,k));

[Onur Aytan]

Thank you, your final replacement for derivatives, not using syms x helped. I will use your final replacement!

Though what was my fault about using syms x for taking derivative? I wanted to do something different and it was logically working in my mind.

About k and i,

"k" stands for "x" in original statement xi+1 = xi - f / df

I redifined x part as k to prevent confusion with sys x interception. It had nothing to do with i indice after all.