Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Roots - Bisection method

Subject: Roots - Bisection method

From: Susan

Date: 4 Jan, 2011 10:23:04

Message: 1 of 6

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

Subject: Roots - Bisection method

From: Husam Aldahiyat

Date: 4 Jan, 2011 10:57:04

Message: 2 of 6

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

Subject: Roots - Bisection method

From: Husam Aldahiyat

Date: 4 Jan, 2011 11:05:21

Message: 3 of 6

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;

Subject: Roots - Bisection method

From: Susan

Date: 4 Jan, 2011 12:13:04

Message: 4 of 6

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?

Subject: Roots - Bisection method

From: Husam Aldahiyat

Date: 4 Jan, 2011 13:05:05

Message: 5 of 6

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)),'*')

Subject: Roots - Bisection method

From: Susan

Date: 4 Jan, 2011 13:44:05

Message: 6 of 6

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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us