MATLAB: Odd results for a zero finding function after it converges… Please help!

newton/bisect methodzero approximation

This function is a hack and slash version of the Bisection method and Newton's method for finding a zero intersection of an inputted function, Fun, and your initial interval guess. I have no other choice in methods as my teacher explicitly stated to do so.
The function is pretty amazing in how it convergences on xsol, which is the x value at Fun=0…. But after it converges it produces oddly predictable patterns.
If x = 2 for y = 0, xsol finds x=2, but after it converges it produces xsol = -5. If x = 3, xsol = 3, and after it converged again xsol = -4.5 and so on so on x = 3, xsol = 3, but then goes to xsol = -4 after the initial convergence.
Note: my function has no safety features included as i didn't want to complicate it when asking for help. Because of this, if your if your interval is your solution it will mess up, if the derivative is 0, it will mess up, and if a value of an interval limit is zero, it will mess up. This is okay, as i can fix this later.
If none of these features are bothered, the function converges on the correct value before it goes off on a tantrum.
A quick explanation of my variables, a0 and b0 are the initial interval, tol is nothing as i deleted it for now, a is a row vector full of the new a0 interval values, b is also a row vector, but of the new b interval values and the current estimation for xsol.
Any help for explaining this odd pattern (and how to fix it) would be much appreciated!! I'm certainly not worried about function efficiency, only capabilities as I'm quite the matlab noob. If you read this far, you certainly deserve a Thank you!
%%%Script
      Fun = @(x) (x-2);
      a0 = 5;
      b0 = -12;
      tol = 10^-6;
      UntitledEight(Fun,a0,b0,tol);
    %%Function
    function UntitledEight (Fun,a0,b0,tol)
    fa = Fun(a0);
    fb = Fun(b0);
    if (a0 > b0)
        temp = a0;
        a0 = b0;
        b0 = temp;
    end
    a = [0;a0;zeros([101,1])];
    b = [a0;b0;zeros([101,1])];
    n = 100;
    for i = 2:n+2
        xsol = b(i,1);
        if (a(i,1) > b(i,1))
            temp2 = a(i,1);
            a(i,1) = b(i,1);
          b(i,1) = temp2;
      end
      m = (b(i,1) + a(i,1))/2;
      s = b(i,1)-(Fun(b(i,1))*(a(i,1)-b(i,1)))/(Fun(a(i,1))-Fun(b(i,1)));
       if (Fun(b(i,1)) == Fun(b(i-1,1)))
       b(i+1,1) = b(i+1,1) + m;       
              else if  (m < s && s < b(i,1))
              b(i+1,1) = b(i+1,1) + s;
                  else, b(i+1,1) = b(i+1,1) + m;
                  end
       end
       if (Fun(a(i,1))*Fun(b(i,1)) < 1)
           a(i+1,1) = a(i+1,1) + a(i,1);
       else a(i+1,1) = a(i+1,1) + b(i,1);
       end
end
end

Best Answer

incase anyone else finds themselves in this position, this is the final product. Works well
function xsol = NewSol_Corley(Fun,a0,b0,tol)
clear all;
clc;
Fun = @(x) (x-50)*(x+2)*(x+80);
a0 = 2;
b0 = 60;
tol = 10^-100;
% evaluating function at inital interval points
fa = Fun(a0);
fb = Fun(b0);
% error message if the interval doesn't contain a solution or contains an
% even number of solutions.
if (fa*fb > 0) 
    fprintf ('The brackets do not enclose a solution or enclose an even number of solutions');
    return;
end
% these if statements display xsol if one or both of the initial intervals is
% a solution
if (fa*fb == 0)
        % if both are solutions
        if (fa == 0 && fb == 0)
        fprintf('xsol = %d and xsol = %d (the inputs a0 & b0)\n\n',a0,b0);
        return;
            % if a0 is a solution
            elseif (fa == 0)
            fprintf('xsol = %d (the input a0)\n\n',a0');
            return;
                % if b0 is a solution
                else fprintf('xsol = %d (the input b0)\n\n',b0);
                return;
        end
end
% this ensures that a0 and b0 are initially in ascending order by switching
% them if they are not.
if (a0 > b0)
   temp = a0;
   a0 = b0;
   b0 = temp;
end
% a will be the matrix where all the left endpoints are stored.
a = [0;a0;zeros([101,1])];
% b will be the matrix where all the right endpoints are stored, as well as
% the current approximation of xsol.
b = [a0;b0;zeros([101,1])];
% n is the number of iterations.
n = 100;
% this function is the driving for loop, it begins at two for convience
% purposes only.
for i = 2:n+2
    %defining xsol's temporary value
    xsol = b(i);
      % m is the current midpoint for each iteration.
      m = (b(i) + a(i))/2;
      % s is the current result of the secant method for each iteration.
      s = b(i)-(Fun(b(i))*(a(i)-b(i)))/(Fun(a(i))-Fun(b(i)));
% the following statments match the requested critera for defining b(i+1).
% first critera, if the last approximation is the same as the current, then
% the next interval is altered with the m value (bisect method)
if (b(i) == b(i-1))   
        % this if / else statment insures that the correct value is altered


        % in such a way that the solution is still within the new interval.
        if (Fun(m)*Fun(a(i)) <= 0)
        b(i+1) = b(i+1) + m;
        a(i+1) = a(i);
            else
            a(i+1) = b(i+1) + m;
            b(i+1) = b(i);
        end
% second critera, if the function values of the current intervals multiplied are
% negative, then the endpoints are adjusted with the s value (secant method)
elseif (Fun(a(i))*Fun(b(i)) < 0)
      % this if / else statment insures that the correct value is altered
      % in such a way that the solution is still within the new interval.        
      if (Fun(s)*Fun(a(i)) <= 0)
      b(i+1) = b(i+1) + s;
      a(i+1) = a(i);
          else
          a(i+1) = b(i+1) + s;
          b(i+1) = b(i);
      end
% third critera inclues all other possible outcomes
else
      % this if / else statment insures that the correct value is altered
      % in such a way that the solution is still within the interval.        
      if (Fun(m)*Fun(a(i)) <= 0)
      b(i+1) = b(i+1) + m;
      a(i+1) = a(i);
          else
          a(i+1) = b(i+1) + m;
          b(i+1) = b(i);
      end
end
% this if statement checks if the current approximation is within tolerance
% and feeds back occordingly.
    if (abs(Fun(b(i+1))) < tol)
    fprintf('xsol was found in %d iterations\n\n', i);
    xsol = b(i+1)
    return
    end
    end
    % if the tolerance is not met within the interval, this outputs the current
    % best guess.
    fprintf('xsol was not found within the tolerance using %d iterations\n\n',n);
    xsolapprox = xsol
    end
Related Question