Thread: My Matlab function M-file for the Secant Method

Hi guys.

I've recently wrote this Function M-file to find the root of any function using the secant method. I have a problem though, in that whenever i use the file on a function, i get the error

???Error using ==> secant at 34
convergence not achieved

Can anyone see why my file is doing this?
Here is the script:

function [x, iter] = secant( f, x0, x1, nmax, tol )
%SECANT Secant method for root finding.
% [X, ITER] = SECANT( F, X0, X1, NMAX, TOL ) finds a zero, X,
% of the function F(X) using the secant method. X0 and X1 are the
% initial values required for the secant method, NMAX is the
% maximum number of iterations allowed, and TOL is the convergence
% criterion. ITER is the number of iterations required.

% Test to make sure nmax and tol are valid.
if nmax < 1
error('nmax must be greater than 0')
end
if tol < eps
error('tol must be greater than or equal to eps')
end
% nmax >=1 and tol >= eps at this point.

%Set initial valies for x, iter and err.
x = x1;
iter = 0;
err = abs(x1-x0);

% The invariant for this loop is "iter is the number of iterations
% completed at the start of the loop".

% At the start of the loop, x0 and x1 are the first two points which will be
% used to extrapolate a secant line to the next value of x. The current
% approximation to the solution (after 0 iterations) is x = x1 and err is an
% upper bound for the error in x.

while err > tol
% Test for maximum number of iterations.
if iter == nmax;
error('convergence not achieved')
end
% Carry out one secant step.
fx1 = f(x1);
fx0 = f(x0);
x = x1 - (fx1*((x1 - x0)/(fx1 - fx0))); % Secant forumula
x1 = x0; % x1 is new value for x0
x = x1; % x is new value for x1
% Increment iter so that the invariant is preserved.
iter = iter + 1;
end
% Here, the while condition is false so err <= tol.
% From the loop invariant, iter is the number of iterations taken.
end

Originally Posted by nugiboy

Hi guys.

I've recently wrote this Function M-file to find the root of any function using the secant method. I have a problem though, in that whenever i use the file on a function, i get the error

> In secant at 39
Warning: Divide by zero.
> In secant at 39
Warning: Divide by zero
> In secant at 39
Warning: Divide by zero.

Can anyone see why my file is doing this?
Here is the script:

function [x, iter] = secant( f, x0, x1, nmax, tol )
%SECANT Secant method for root finding.
% [X, ITER] = SECANT( F, X0, X1, NMAX, TOL ) finds a zero, X,
% of the function F(X) using the secant method. X0 and X1 are the
% initial values required for the secant method, NMAX is the
% maximum number of iterations allowed, and TOL is the convergence
% criterion. ITER is the number of iterations required.

% Test to make sure nmax and conv are valid.
if nmax < 1
error('nmax must be greater than 0')
end
if tol < eps
error('tol must be greater than or equal to eps')
end
% nmax >=1 and tol >= eps at this point.

%Set initial valies for x, iter and err.
x = x1;
iter = 0;
err = abs(x1-x0);

% The invariant for this loop is "iter is the number of iterations
% completed at the start of the loop".

% At the start of the loop, x0 and x1 are the first two points which will be
% used to extrapolate a secant line to the next value of x. The current
% approximation to the solution (after 0 iterations) is x = x1 and err is an
% upper bound for the error in x.

while err > tol
% Test for maximum number of iterations.
if iter == nmax;
error('convergence not achieved')
end
% Carry out one secant step.
fx1 = f(x1);
fx0 = f(x0);
x = x1 - (fx1*((x1 - x0)/(fx1 - fx0))); % Secant forumula
x1 = x0; % x1 is new value for x0
x = x1; % x is new value for x1
% Increment iter so that the invariant is preserved.
iter = iter + 1;
end
% Here, the while condition is false so err <= tol.
% From the loop invariant, iter is the number of iterations taken.
end

What do you thing are the value of x0 and x1 after this:

Code:

x = x1 - (fx1*((x1 - x0)/(fx1 - fx0))); % Secant forumula 
x1 = x0; % x1 is new value for x0
x = x1; % x is new value for x1

??

CB

Originally Posted by CaptainBlack

What do you thing are the value of x0 and x1 after this:

Code:

x = x1 - (fx1*((x1 - x0)/(fx1 - fx0))); % Secant forumula 
x1 = x0; % x1 is new value for x0
x = x1; % x is new value for x1

??

CB

I'm saying that the new secant line is extrapulated using the old x1 value and the new x value from the last iteration.
Is that not right?

ahh ive got them the wrong way around havn't i?
It should be

x = x1 - (fx1*((x1 - x0)/(fx1 - fx0))); % Secant forumula
x0 = x1; % x1 is new value for x0
x1 = x; % x is new value for x1

EDIT:
This still doesn't work. I'm getting:

??? Error using ==> secant at 34
convergence not achieved

Any ideas anyone? I need to know quite urgently!