MatLab - Finding root of f(x) = x - tan(x) with bisection method

Question

I have written a code for the bisection algorithm in MatLab. I have based this on the pseudocode given in my textbook. The algorithm has worked just fine on all my problems so far, but when I'm asked to find a root of f(x) = x - tan(x) on the interval [1,2] I have some troubles. My code is as follows:

function x = bisection(a,b,M)
f = @(x) x - tan(x);
u = f(a);
v = f(b);
e = b-a;
x = [a, b, u, v]
if (u > 0 && v > 0) || (u < 0 && v < 0)
    return;
end;
for k = 1:M
    e = e/2;
    c = a + e;
    w = f(c);
    x = [k, c, w, e]
    if (abs(e) < 10^(-5) || abs(w) < eps)
        return;
    end
    if (w < 0 && u > 0) || (w > 0 && u < 0)
        b = c;
        v = w;
    else
        a = c;
        u = w;
    end
end

If I run this algorithm on the interval [1,2] with, say, 15 iterations, my final answer is:

x =

  1.0e+004 *

    0.0015    0.0002   -3.8367    0.0000

which is obviously way off as I wish to get f(c) = 0 (the third entry in the vector above).

If someone can give me any help/tips on how to improve my result, I would greatly appreciate it. I am very new to MatLab, so treat me as a novice :).

Answer 1

Let's have a look at the sequence of c which is generated by the bisection method:

c =  1.5000
c =  1.7500
c =  1.6250
c =  1.5625
c =  1.5938
c =  1.5781
c =  1.5703
c =  1.5742
c =  1.5723
c =  1.5713
c =  1.5708
c =  1.5706
c =  1.5707
c =  1.5707
c =  1.5708

You can see that it converges to pi/2. The tan() has a singularity at this point and so does x - tan(x). The bisection method converges to this singularity as is also stated here, for example. This is why the function value at f(c) is not close to zero. In fact it should go to (plus/minus) infinity.

Other suggestions:

I like your bisection method. To make it usable in a more general fashion, you could incorporate these changes:

• Make the variable f a function parameter. Then you can use your method with different functions without having to re-write it.
• The very first time you assign x = [a, b, u, v] but later on x(1) is the number of iterations. I would make this more consistent.
• You can easily test if f(a) and f(b) have different signs by looking at the sign of the product p = f(a)*f(b). For p > 0 the signs are equal and there is no root, for p < 0 the signs differ and for p == 0 either f(a) or f(b) is zero and you already have found the root.