数学算法

				Function ZBRENT(X1, X2, TOL:real):real;
				var
				    ITMAX,ITER:integer;  ZZ,S,AAA,EPS,FA,FB,FC,C,TOL1,XM,P,Q,R,D,E,A,B:real;
				begin
				    ITMAX:=100;
				    EPS:=0.00000003;
				    A:=X1;
				    B:=X2;
				    FA:=FUN(A);
				    FB:=FUN(B);
				    If FB * FA > 0  Then
				        ShowMessage('Root must be bracketed for ZBRENT');
				    FC:=FB;
				    For ITER:=1 To ITMAX do
				    begin
				        If FB * FC > 0 Then
				        begin
				            C:=A;
				            FC:=FA;
				            D:=B - A;
				            E:=D;
				        end;
				        If Abs(FC) < Abs(FB) Then
				        begin
				            A:=B;
				            B:=C;
				            C:=A;
				            FA:=FB;
				            FB:=FC;
				            FC:=FA;
				        end;
				        TOL1:=2 * EPS * Abs(B) + 0.5 * TOL;
				        XM:=0.5 * (C - B);
				        If (Abs(XM) 				        begin
				            ZBRENT:=B;
				            Exit;
				        end;
				        If (Abs(E) >= TOL1) And (Abs(FA) > Abs(FB)) Then
				        begin
				            S:=FB / FA;
				            If A = C Then
				            begin
				                P:=2 * XM * S;
				                Q:=1 - S;
				            end
				            Else
				            begin
				                Q:=FA / FC;
				                R:=FB / FC;
				                P:=S * (2 * XM * Q * (Q - R) - (B - A) * (R - 1 ));
				                Q:=(Q - 1 ) * (R - 1 ) * (S - 1 );
				            end;
				            If P > 0  Then Q:=-Q;
				            P:=Abs(P);
				            If 3 * XM * Q - Abs(TOL1*Q) < Abs(E*Q) Then
				                AAA:=3  * XM * Q - Abs(TOL1 * Q)
				            Else
				                AAA:=Abs(E * Q);
				            If 2 * P < AAA Then
				            begin
				                E:=D;
				                D:=P / Q;
				            end
				            Else
				            begin
				                D:=XM;
				                E:=D;
				            end;
				        end
				        Else
				        begin
				            D:=XM;
				            E:=D;
				        end;
				        A:=B;
				        FA:=FB;
				        If Abs(D) > TOL1 Then
				            B:=B + D
				        Else
				        Begin
				            If XM>=0 then
				                ZZ:=1
				            Else
				                ZZ:=-1;
				            B:=B + Abs(TOL1) * ZZ;
				        end;
				        FB:=FUN(B);
				    end; 
				    ShowMessage('ZBRENT exceeding maximum iterations.');
				    ZBRENT:=B;
				end;