%
%   bisection method to find a root
%   problem: swinging mass on a string
%   inputs:  T (tension force), L (length), ms (mass of string)
%            m (hanging mass)
%            x1, x2 (start values for bracketing)
%   output:  x (frequency)

%    Note:  this is a very demanding problem for the bisection method
%    because of the cot() singularities.  The  method 'thinks' that there
%    are additional roots where the cot() jumps from -inf to +inf.  That
%    makes the plotting difficult, as well as finding an initial bracket.
%    Newton-Raphson will always converge, but maybe not the root you want
%    to find.

clear;                 % clean slate
x1=0.714;  x2=0.715;   %  second root, initial guess
x1=0; x2=0.5;          %  first root

m = 200;                %  hanging mass in kg
T = m*9.81;             %  200 kg times g --> tension in Newtons
L = 25;                 %  rope length in m
ms = 40;                %  mass of rope
c = sqrt( T / (ms/L));  %  calculate wave speed on the string 
prec=1.0e-7;            %  desired precision

%  define an inline function for which we will find the root

F=inline('cot(2*pi*f*L/c)-2*pi*f*m*L/(c*ms)','f','L','c','m','ms');

%...... plot function (not good for this case because of singularities)
%       misses zeroes unless the step is very small
xx=0:0.001:2;  yy=F(xx,L,c,m,ms);  %  make arrays for plotting
clip=1; yy(yy>clip)=clip;  yy(yy<-clip)=-clip; %.. clip values because cot() goes to +-inf 
plot(xx,yy); 
%..... do a fine grained search
dx=1.0e-4; x=0; f1=F(x,L,c,m,ms); while x<3; x=x+dx; f2=F(x,L,c,m,ms);
    if f1*f2 < 0; fprintf('root between %f and %f\n',x-dx,x); end;
    f1=f2; end

f1=F(x1,L,c,m,ms);  f2=F(x2,L,c,m,ms);  count=0;  %  initial values
if f1*f2 > 0; disp('root not bracketet'); return; end;  %  not a good initial guess
while 1
    xm=0.5*(x1+x2);  count=count+1; %   middle of interval
    fm=F(xm,L,c,m,ms);              %   function value at midpoint
    if fm < 0;  x2=xm; f2=fm;       %   new right point
    else  x1=xm; f1=fm;             %   new left point
    end
    if x2-x1 < prec; break; end       %   converged
end

xm=0.5*(x1+x2);                     %   mddle is better
fprintf('converged:  %d iterations  frequency:  %f Hz, perid is %f seconds\n',count,xm,1/xm);


