Alex Pozo
Alex Pozo

Reputation: 168

Octave, The secant method

I am trying to implement the secant method using Octave, and in particular I made a function called secant("function","x0","x1","tolerance"). I used the while-loop for calculate the root of a function and I think the problem is in the loop.

But the result that I get, it is not what I expect. The correct answer is x=0,49438.

My code is the following:

tol = 10.^(-5); # tolerance
x0 = 0.4; # initial point
x1 = 0.6; # initial point
syms x;
fun = @(x) exp(sin(x)) - 2/(1+x.^2); # f(x) = exp(sin(x)) - 2/(1+x^2)
fprintf("Function f(x)\n");
fun
fprintf("\n");
fprintf("initial points: %f\n",x0,x1);
function[raiz,errel,niter] = secant(fun,x0,x1,tol)      
  niter = 0; # steps number
  errel = [];  # the vector of relative errors
  y0 = fun(x0); # f(x0)
  y1 = fun(x1); # f(x1)
  ra = 0.0; # ra is the variable of the function's root
  while abs(ra-x1)>= tol
    niter += 1;
    ra = x1 - ((x1-x0)./(y1-y0)).*y0; # formula of the secant method
    if abs((ra-x1))<tol
      raiz = ra;
      return;
    endif
    x0 = x1; y0 = y1; x1 = ra;
    y1 = fun(ra);
    errel(niter) = abs(ra-x1)/abs(ra); # Calcule of relative error
  endwhile
  if abs((ra-x1)/ra)<tol
    fprintf ('The method is over\n');
    fprintf ('\n');
  endif
  raiz = ra;
endfunction
[r,er,tot] = secant(fun,x0,x1,tol)

I appreciate the help you can give me

Upvotes: 0

Views: 654

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

It makes little sense to use the secant root in the loop condition. At first it is not defined, and inside the loop it is shifted into the support points for the next secant. Note that at the end ra and x1 contain the same value, making the loop condition trivial, in a wrong way.

Next, the secant has the equation

y(x) = y1 + (y1-y0)/(x1-x0)*(x-x_1)

check that with this formula y(x0)=y0 and y(x1)=y1. Thus the secant root is to be found at

x = x1 - (x1-x0)/(y1-y0)*y1

Finally, at no point are symbolic expressions used, defining x as symbolic variable is superfluous.

The break-out test inside the loop is also redundant and prevents a coherent state after the loop. Just remove it and remember that x1 contains the last approximation.

With all this I get an execution log as follows:

Function f(x)
fun =

@(x) exp (sin (x)) - 2 / (1 + x .^ 2)

initial points: 0.400000
initial points: 0.600000
The method is over

r = 0.494379048216965
er =

   2.182723270633349e-01   3.747373180587413e-03   5.220701832080676e-05   1.899377363987941e-08

tot = 4

Upvotes: 1

Related Questions