How to compute Newton Raphson coding on Matlab?

I tried to compute Newton Raphson Algorithm on Matlab to calculate the root estimate for y but the final answer x pop is same as the initial guess which i think is wrong. I don't know where the mistake is in my script. Is there anything missing in my script?

if true
  % code
end
x=0.00462; % initial guess
emax=0.0000001; % maximum error (accuracy)
imax=10000; % maximum iteration
dmin=0.000001; % minimum slope
deltax=inf; % approximate error
i=1; % step number
while (i<imax) && (deltax>=emax);
  y = x.^5-16.5*x.^4+102.85*x.^3-299.475*x.^2+401.1634*x-194.0612;
  y1 = 5*x.^4-66*x.^3+308.55*x.^2-598.95*x+401.1634;
  x1=x-(y/y1);
  deltax = abs((x1-x)/x1);
  if abs(y1)<=dmin;
      disp('Too small slope');
      break;
  end
  i=i+1;
  x1=x;
end
disp(x)

[Jim Riggs]

The first thing missing from your script is comments. I suggest that you add comments which describe what each line is doing. Then see if that makes sense.

Answer 1

[Jim Riggs]

Here is how I would code it:

xguess = 0;     % initial guess
emax   = 1e-7;  % Maximum error
imax   = 1000;   % maximum number of iterations
maxstep = 0.5;   % maximum x-step size
minslope = .01;  % minimum value for the slope (prevents divide by zero)

func = @(x) x.^5-16.5*x.^4+102.85*x.^3-299.475*x.^2+401.1634*x-194.0612;  % User Function
dfunc= @(x) 5*x.^4-66*x.^3+308.55*x.^2-598.95*x+401.1634;                 % User Function derivative
lim  = @(x,xlim) max(-xlim,min(xlim,x));   % limiter function

i=1;
x=xguess;

while (i<imax) && (abs(deltax)>=emax    % needs abs(deltax) (can be + or -)
    yguess = fun(x);              % Get the function value for the x guess
    slope = dfunc(x);             % Get the function slope for x-guess
    if(abs(slope)<minslope)       % Don't allow slope to go to zero (trap divide by zero condition)
        slope = minslope*sign(slope);
    end
    yerr = yguess;                % Since the desired value is zero, the function value represents the error value 
    deltax = -yerr/slope;         % calculate the x-step
    xstep = lim(deltax,maxstep);  % limit the x-step (deltax) to +/- maxstep
    x = x + xstep;                % apply the (limited) x-step to the x-guess value and repeat the iteration
    i=i+1;                        % update the increment counter
end

Notice the use of "abs" on the error threshold test.

Also, notice that I used a limiter on the x-step. This is helpful (even necessary) where the slope is very small and the calculated x-step could be very large. For example, if you use an initial guess of 1.5, the function value is 5.0464 and the slope is -0.4616. This would give an x-step of -5.0464/-0.4616 = 10.9324, and the new x value would be 12.4324. In this case, it would not find the root at 2.1141. Using a limit of 0.5 results in a new x-value of 2.0, and the algorithm would successfully converge on the right answer.

Plotting the function shows you where the initial guess values should be. I used initial guess values of 0, 2, 3.5, 4.5, and 6.

I get roots at 1.1238 (9 iterations), 2.1141 (5 iterations), 3.4394 (5 iterations), 4.3009 (6 iterations), and 5.5219 (6 iterations)