How to plot convergence graph using bisection method code

I am having trouble getting the convergence graph to plot using this code for bisection method.

function xc=bijection(f,a,b,tol)
f=@(x) x-cos(x);
a=0;
b=1;
fa=f(a);
fb=f(b);
tol=(0.5).*10.^(-12);
if sign(fa)*sign(fb)>0
    error('f(a)f(b)<0 not satisfied')
end

while (b-a)/2>tol
    c=(a+b)/2;
    fc=f(c);
    if fc==0
        break;
    elseif sign(fc)*sign(fa)<0
        b=c;
        fb=fc;
    else
        a=c;
        fa=fc;
    end
end
xc=(a+b)/2;
disp(c)
format long
plot(???)

I do not exactly know what I am suppose to graph. I have yet to learn this in class, I'm just trying to do a little extra. Any help would be greatly appreciated. Thank you.

Answer 1

[Geoff Hayes]

Sarah - as the bisection method is a root finding algorithm, I suspect that what you want to plot (from each iteration of the while loop) is either c or f(c) as you would want to show convergence to either the the root, c, or zero.

A couple of things - you have named your function bijection rather than bisection, and you have initialized your inputs in the function body whereas they should be initialized prior to calling the function and passed in as the inputs:

 f=@(x) x-cos(x);
 a=0;
 b=1;
 fa=f(a);
 fb=f(b);
 tol=(0.5).*10.^(-12);

Since the above variables are declared in the base workspace, now call your function from the command line as

 bisection(f,a,b,tol)

Now, you will want to collect your c and f(c) at each iteration of the while loop. So you can do this by either creating an array (for both) and re-sizing it at each iteration of the loop, or you can pre-size the arrays given the maximum number of iterations for your while loop. You don't currently have one, but it is a good thing to have so that you don't get stuck in the loop. Try the following

 maxIters = 100;
 atIter   = 0;
 cVals    = zeros(1,maxIters); 
 while (b-a)/2>tol && atIter < maxIters 
      c=(a+b)/2;
      fc=f(c);

      atIter        = atIter + 1;
      cVals(atIter) = c;

      if fc==0
          break;
      elseif sign(fc)*sign(fa)<0
          b=c;
          fb=fc;
      else
          a=c;
          fa=fc;
      end  
 end
 xc=(a+b)/2;
 disp(c)
 format long
 cVals = cVals(1:atIter);

See in the above how we create the array of cVals, pre-sizing it for 100 elements (or the max 100 iterations of the while loop). At each iteration of the loop, we update the cVals array and increment the atIter counter. In the iteration condition, we also check to see if we have reached the maximum number of iterations - if so, then the loop is finished.

Outside of the loop, we do

 cVals = cVals(1:atIter);

so that we remove all unused elements of cVals (for the case where we exit the loop prior to iterating maxIters times.

You can now plot the convergence to the root as

 plot(cVals)

or the convergence to zero as

 plot(f(cVals))

or both.

One final thing to consider, since you are using doubles, your equality check for

 fc==0

is not guaranteed to succeed. The above condition is fine for integers, but when it comes to using floating precision values, you need to do a tolerance check instead (like you are doing in the condition of the while loop). Something like

 if abs(fc)<eps

is probably sufficient.

[Sarah Harvey]

Thank you so much for catching the bijection/bisection mistake, and thank you for all the help. Now, you have me set a max number of iterations. What if I don't want a set number of iterations or is that not possible? The code I have now works and I'm pretty sure it gives me the correct answer, I'm just not sure how to plot. So my other question is: Is my code correct the way it is (minus the bijection/bisection mistake)? I am not exactly sure if the way you're showing me is an extension so that I can plot or if you're showing me another way. I hope that doesn't come off the wrong way, it is just that this is the first semester I've ever done this and used MATLAB, so this is pretty much new to me. Also, it is a presentation I have to present to the class; therefore, it has to be what I've learned so far. I'm just trying to extend it to plotting it just because I wasn't to learn it. However, again thank you for all the help, I am definitely trying what you said, just to learn that way also. Thank you!

[Geoff Hayes]

Sarah - when using while loops, it is good practice to set a maximum number of iterations for the loop so that you don't become stuck in it (which is a common problem). It may be that the inputs you provide will be such that this issue never arises but you may want to err on the side of caution. You can always change to this value from a 100 to 1000 or more.

As for plotting the convergence, using plot is typically the function to use when plotting a 2D line (which will show your convergence). If you haven't learned this function in class, then you may want to consult with your professor to find out an alternative way to plot data and to confirm which data you should be plotting.

[Sarah Harvey]

Thank you. I will try one with what you said above and see what the teacher says about presenting that. Also, he taught us about using plot for vectors, I'm just not sure how to use it for iterations, what do I put in plot(???)? I want to plot the line converging to the root and the original function in the same plot. It would be:

Plot(something,function)

However, I an a little confused on what that something is, and when I try plugging f or f(x) in for function it says error. Is there something I put in plot() to do what I'm trying to do?

[Geoff Hayes]

Sarah - I'm not sure if plotting the convergence to the root and the convergence to zero on the same graph/figure is appropriate given that one converges to zero and the other converges to the root which could be any number (you would not necessarily observe the convergence in a "pleasant" manner). I think that it would be better if you have two figures or subplots for each convergence.

As for the inputs to plot you would want to show the iterations along the x-axis and the root (or function evaluation) along the y-axis. So you could do as I have shown you, or if you want to be explicit, you could do

 plot(1:length(cVals),cVals);

and

 plot(1:length(cVals),f(cVals));

[Sarah Harvey]

Ok, I think I get it now. Thank you again for all the help. I'll keep mine without the plot, for the presentation, and I will try the other one, with plot, using what you have said, for my own practice. Thank you!

[MAJID KHAN]

sarah harvey i need help in this