bisectio.m

这是在网上下的一个东东

%%%%%%%%%% start of function bisection.m %%%%%%%%%%
%
% Function to calculate the values of the Lagrange 
% multipliers, beta, using bisection method
%
% Written by: Roseanna M. Neupauer
% Modification Date: April 24, 1999
%
% function beta=bisection(upper,expvalue,leftbegin,...
%     rightbegin,maxiter,tolbeta,nt,nearzero,large)
% 
% Inputs
%     upper        array containing upper limit of prior 
%                    distributions
%     expvalue     array containing prior expected value
%     leftbegin    lower limit of range for bisection method
%     rightbegin   upper limit of range for bisection method
%     maxiter      maximum number of iterations in beta-fitting
%     tolbeta      tolerance for beta-fitting
%     nt           number of elements in the time array 
%     nearzero     value below which the asymptotic 
%                    approximation to zero is used
%     large        value above which the asymptotic 
%                    approximation to infinity
%
% Outputs
%     beta         array containing the Lagrange multipliers
%
% Functions called
%     snfcn        evaluates the function whose zero will be 
%                    found using the bisection method

function beta=bisection(upper,expvalue,leftbegin,rightbegin,...
          maxiter,tolbeta,nt,nearzero,large)

beta=zeros(nt,1);

for i=1:nt
     u=upper(i);
     sn=expvalue(i);
     left=leftbegin;
     right=rightbegin;
     for j=1:maxiter
          bold=(left+right)/2;
          snfcnl=snfcn(left,u,sn,nearzero,large);
          snfcnr=snfcn(right,u,sn,nearzero,large);
          snfcnb=snfcn(bold,u,sn,nearzero,large);
          if (abs(2*sn - u)<eps) 
               beta(i)=0.;
               break
          elseif (snfcnl < 0 & snfcnr < 0)
               if (snfcnl > snfcnr)
                    beta(i)=leftbegin;
               else
                    beta(i)=rightbegin;
               end
               break
          elseif (snfcnl > 0 & snfcnr > 0)
               if (snfcnl > snfcnr)
                    beta(i)=rightbegin;
               else
                    beta(i)=leftbegin;
               end
               break
          else
               if (snfcnb*snfcnl) > 0 
                    left=bold;
               else
                    right=bold;
               end
               if (right-left) < tolbeta 
                    beta(i)=(left+right)/2;
                    break
               else
                    if (j == maxiter)
                      disp('Beta iterations did not converge')
                      beta(i)=(right+left)/2;
                    end
               end
          end
     end % for j
     if (beta(i) == leftbegin & u-sn ~= 0)
          beta(i)=-1/(u-sn);
     end
end % for i
%%%%%%%%%%  end of function bisection.m  %%%%%%%%%%