Hi all,

I have created a simple sine plot and I wish to implement the bisection method to yield all the intersection points along the x-axis (i.e. where y = 0);  however, my code is simply bisecting the whole range and splitting it instead of checking the range at intervals.  If run, the code will work, but to only converge at one point.

I would appreciate any help!!

Many thanks!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

clc;

fo = 10;
Fs = 1000;
Ts = 1/Fs;
i = 0:Ts:1-Ts;
n = length(i);
SI = 2*sin(2*pi*fo*i);

yy = 0;

plot(i,SI,i,yy);

I = @(x) 2*sin(2*pi*fo*x);

a = 0;
b = 0.999;
TOL = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NumIters = 0; % initialize iteration counter
MaxIters = 10;

for x = a:0.001:b

    fa = feval(I,a);
    fb = feval(I,b);

    sfa = sign(fa);
    sfb = sign(fb);

    if (sfa*sfb > 0)
     % this doesn't satisfy the sufficiency condition for existence of
     % a root
     error('The bisection method requires f(a)f(b) < 0 ');
    end

    % bisection loop
    % disp(sprintf('\t iterate \t\t interval \t\t\t\t\t p'));
    while (NumIters <= MaxIters)
    % use midpoint as the iterate  
    p = (a+b)/2;
    fp = feval(I,p);
    sfp = sign(fp);

    % display some data to the screen
    disp (sprintf ('\t %3d \t (%.10f,%.10f) \t %.10f', NumIters, a, b, p ) )
  
    if ( (b-a)<2*TOL | fp == 0)
        x = p;
    return % terminate this function early
    elseif (sfa*sfp < 0)
        b = p;
    else
        a = p;
    end

    NumIters = NumIters + 1;
    end

% if this while loop terminates before returning the root
disp(' max number of iterations reached ... ')
x = a;

end

You're going about this the wrong way. See what I've done here and edit accordingly.

clc;

fo = 10;
Fs = 1000;
Ts = 1/Fs;
i = 0:Ts:1-Ts;
n = length(i);
SI = 2*sin(2*pi*fo*i);

yy = 0;

plot(i,SI,i,yy);

I = @(x) 2*sin(2*pi*fo*x);

TOL = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MaxIters = 10;

for a = 0:.01:1
    
    % intervals tested are 0-.03, .01-.04, .02-.05 ...1-1.03
    b = a + .03;

    fa = feval(I,a);
    fb = feval(I,b);

    sfa = sign(fa);
    sfb = sign(fb);
    
    if (sfa*sfb > 0)
        continue
    end

     disp(['Interval: [', num2str(a),', ',num2str(b),']'])
    
    % bisection loop
    % disp(sprintf('\t iterate \t\t interval \t\t\t\t\t p'));
    an = a;
    bn = b;
    
    NumIters = 0; % initialize iteration counter
    
    while (NumIters <= MaxIters)
    % use midpoint as the iterate
    p = (an+bn)/2;
    fp = feval(I,p);
    sfp = sign(fp);

    % display some data to the screen
    disp (sprintf ('\t %3d \t (%.10f,%.10f) \t %.10f', NumIters, a, b, p ) )
  
     fa = feval(I,an);
    fb = feval(I,bn);

    sfa = sign(fa);
    sfb = sign(fb);
    
    if ( (b-a)<2*TOL | fp == 0)
        x = p;
    break % terminate this function early
    elseif (sfa*sfp < 0)
        bn = p;
    else
        an = p;
    end

    NumIters = NumIters + 1;
    end

% if this while loop terminates before returning the root
disp(' max number of iterations reached ... ')
x = p;

end 

BTW, my intervals overlap each other, so you will see a lot of duplicate roots. 
To change this, edit the for loop index into something else.

For example, to overcome overlapping, you can use the following lines:

    a = 0:0.03:1
    b = a + .03;

Husam,

Thank you very much for your repl(ies).

I have implemented what you suggested and indeed the code outputs periodic roots.  However, I am trying to plot these over my original signal SI by using the following 

if ( (bn-an)<2*TOL | fp == 0)
        x = p
        hold on;
        plot(i,SI,i,x)

etc..

However, I am not getting what I require which is single points at periodic intervals.  Have you any suggestions?

Yes.
clc;

fo = 10;
Fs = 1000;
Ts = 1/Fs;
i = 0:Ts:1-Ts;
n = length(i);
SI = 2*sin(2*pi*fo*i);

yy = 0;

plot(i,SI,i,yy);

I = @(x) 2*sin(2*pi*fo*x);

TOL = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MaxIters = 10;

X = [];

for a = 0:.03:.9
    
    % intervals tested are 0-.03, .01-.04, .02-.05 ...1-1.03
    b = a + .03;

    fa = feval(I,a);
    fb = feval(I,b);

    sfa = sign(fa);
    sfb = sign(fb);
    
    if (sfa*sfb > 0)
        continue
    end

     disp(['Interval: [', num2str(a),', ',num2str(b),']'])
    
    % bisection loop
    % disp(sprintf('\t iterate \t\t interval \t\t\t\t\t p'));
    an = a;
    bn = b;
    
    NumIters = 0; % initialize iteration counter
    
    while (NumIters <= MaxIters)
    % use midpoint as the iterate
    p = (an+bn)/2;
    fp = feval(I,p);
    sfp = sign(fp);

    % display some data to the screen
    disp (sprintf ('\t %3d \t (%.10f,%.10f) \t %.10f', NumIters, a, b, p ) )
  
     fa = feval(I,an);
    fb = feval(I,bn);

    sfa = sign(fa);
    sfb = sign(fb);
    
    if ( (b-a)<2*TOL | fp == 0)
        x = p;
    break % terminate this function early
    elseif (sfa*sfp < 0)
        bn = p;
    else
        an = p;
    end

    NumIters = NumIters + 1;
    end

% if this while loop terminates before returning the root
disp(' max number of iterations reached ... ')
x = p;
        X = [X;x];

end 
hold on
plot(X,zeros(size(X)),'*')

Many thanks for your help Hasam, objectives achieved. :)