zero.m

如何求解函数的零点的两个子函数

%function  ZERO (X1,X2,X0,ERROR,DIFF,IMAX,I,F,FLAG)

%     ***************************************************************
%
%             L. C. ANDREANI AND A. PASQUARELLO, PISA, 1986
%
%
%     Finds the zero of a function F(x) in a specified interval
%     by straight bisection. To be used after subroutine "SETACC".
%
%     Arguments:
%
%     X1,X2 (ENTRY) = Extrema of interval (unmodified)
%     X0    (EXIT)  = Zero of function
%     ERROR (ENTRY) = Error parameter (unmodified)
%     DIFF  (EXIT)  = actual error; abs(x0-x0true)<diff
%     IMAX  (ENTRY) = Maximum number of iterations (unmodified)
%     I     (EXIT)  = Actual number of iterations
%     F     (ENTRY) = Name of external function F(x) (unmodified)
%     FLAG  (EXIT)  = set equal to 1 if the function has equal sign
%                     at extrema, otherwise set to 0...
%
%     Other variables:
%
%     F0,F1,F2 = Function values
%     S0,S1,S2 = Function signs
%
%     ****************************************************************
%
   %   IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   clear all
   clc
   
   X1=6;
   X2=8;
   ERROR=1e-10
  
   IMAX=100000;
      I=0 ;
      FLAG=0 ;
      F1=X1^3-6*X1^2-1/6
      F2=X2^3-6*X2^2-1/6
      S1=sign(F1)
      S2=sign(F2)
      if (F1==0)    %  x1即区间左端点上函数值为零
         X0=X1 ;
         DIFF=0 ;
          VALUE=X0^3-6*X0^2-1/6
         return
      elseif (F2==0.) 
         X0=X2 ;     %  x1即区间右端点上函数值为零
         DIFF=0 ;
          VALUE=X0^3-6*X0^2-1/6
         return
      elseif (S1==S2) 
         FLAG=1 ;     % 该区间没有零点，FLAG=1就是指在这个区间里没有根。
         return
      end
for jjj=1:IMAX+1
    X0=(X2+X1)/2 ;
      F0=X0^3-6*X0^2-1/6
      S0=sign(F0)
      I=I+1 ;
      if (F0==0.) 
         DIFF=0 ;
         VALUE=X0^3-6*X0^2-1/6
         return 
         elseif  S0==S1 
         X1=X0 ;
         F1=F0 ;
         else
         X2=X0 ;
         F2=F0 ;
     end
      if ((abs(X2-X1)<ERROR)|(I==IMAX)) 
         X0=(X2+X1)/2 ;
         DIFF=abs((X2-X1)/2)
         if (DIFF>=ERROR) 
           fprintf( 'FROM MODULE ZERO: MAXIMUM STEPS TAKEN')
         end
          VALUE=X0^3-6*X0^2-1/6
         return
      end
end