Jeremy Fisher
Jeremy Fisher

Reputation: 2782

Vector in Matlab not populating correct values

I am trying to test the efficiency of using two methods of finding the roots of a function in Matlab. In order to do this I am counting the iterations each method takes and then plotting it for various tolerance values (how close f(root) is to 0). For some reason, when I am calculating the iteration values and saving them to a vector - because there are different values for different tolerances - the vector seems to only be populated with 0. I am confused by this because the function I wrote should not be returning 0 as an iteration value. Here is all my code:

function cooldrive()
    a = 0.1; b = 3;
    tol = 0.01;
    L = 1:0.1:4;
    r_x = zeros(1, length(L));
    ctr = 1;
    for i=L
        [x, iter] = bis(a, b, tol, i);
        r_x(ctr) = x;
        ctr = ctr+1;
        b = x;
    end

    r_x = r_x.^2;
    figure
    plot(L, r_x, '.-b', 'markersize', 16)
    ylabel('Decay Rate, x^2'), xlabel('Length, L'), title('Cooling Rate')
    grid

    x = (a+b)/2;
    r_iter_newt = zeros(1, 8);
    r_iter_bis = zeros(1, 8);
    r_tol = zeros(1, 8);
    ctr = 1;
    for j=1:8
        tol = 10^(-1*j);
        [x_newt, iter_newt] = newt(x, tol, 1);
        [x_bis, iter_bis] = bis(a, b, tol, 1);
        b = x_bis;
        r_iter_newt(ctr) = iter_newt;
        r_iter_bis(ctr) = iter_bis;
        r_tol(ctr) = tol;
        ctr = ctr+1;
    end

    figure
    plot(r_tol, r_iter_newt, '.-r', 'markersize', 16)
    hold on
    plot(r_tol, r_iter_bis, '.-b', 'markersize', 16);
    hold off
return

%% Bisection Method 
function [x, iter] = bis(a, b, tol, L)
    iter = 0;
    z1 = coolfun(a, L);
    z2 = coolfun(b, L);
    if z1 * z2 > 0
        disp('Root may not exist')
        x = NaN;
        iter = NaN;
        return
    end

    x = (a+b)/2;

    while abs(coolfun(x, L)) > tol
        if (coolfun(a,L)*coolfun(x,L)) <= 0
            b = x;
        else
            a = x;
        end
        x = (a+b)/2;
        iter = iter + 1;
    end 
    %[x, iter] = bis(a, b, tol, L);
return 

%% Newton's Method
function [ox, iter] = newt(ix, tol, L)
    iter = 0;
    while abs(coolfun(ix, L)) > tol
        x_j = ix - (coolfun(ix, L)/coolfundx(ix, L));
        ix = x_j;
        iter = iter+1;
    end
    ox = ix;
return

%% Evaluate function
function val = coolfun(x, L)
    val = sin(x*L) + x*cos(x*L);
return

%% Evaluate Function's Derivative
function val = coolfundx(x, L)
    val = (L*cos(x*L)) - (L*x*sin(x*L)) + cos(x*L);
return

In the above code, the value of iter_bis is always 0, so the values in r_iter_bis are NaN, 0, 0, .... I can't see a reason why this is occurring.

Upvotes: 0

Views: 109

Answers (1)

David
David

Reputation: 8459

The problem is you test fora root to exist. You have

z1 = coolfun(a, L);
z2 = coolfun(b, L);
if z1 * z2 > 0
    disp('Root may not exist')

however since the function you are trying to solve is not monotonic, you can have both end-points positive, but intermediate points negative (or vice-versa). In this case multiple zeroes exist but your test says that they don't, and your code will stop there.

The problem is really with the method. The bisection method does not work for non-monotonic functions. You have to know, based on whether your function value is positive or negative, which way to move to get closer to the zero, but if your function is not monotonic then you don't know.

You either need to restrict the bisection search domain so that the function is monotonic, or use a different method.

Upvotes: 1

Related Questions