brent.f90

此Fortran程序用布伦特方法求解一维的非线性方程,Func为非线性方程的对应函数表达式.

		   !布伦特方法求解非线性方程
           DOUBLE PRECISION FUNCTION ZBRENT(FUNC,X1,X2,TOL,Z,UZ)
           DOUBLE  PRECISION FUNC,X1,X2,TOL,Z,UZ,A,B,C,D,E,FA,FB,FC,XM,TOL1,S,P,Q,R
		   PARAMETER(ITMAX=2000,EPS=3E-8)
           A=X1
           B=X2
		   FA=FUNC(A,Z,UZ)
		   FB=FUNC(B,Z,UZ)
		   IF(FA*FB.GT.0.D0) PAUSE '必须给出正确的有根区间'
		   FC=FB
		   DO  11  ITER=1,ITMAX
		     IF(FB*FC.GT.0)THEN
			   C=A
			   FC=FA
			   D=B-A
			   E=D
			 ENDIF
			 IF(DABS(FC).LT.DABS(FB))THEN
			   A=B
			   B=C
			   C=A
			   FA=FB
			   FB=FC
			   FC=FA
			 ENDIF
			 TOL1=2.D0*EPS*DABS(B)+0.5D0*TOL
			 XM=0.5D0*(C-B)
			 IF((DABS(XM).LE.TOL1).OR.(FB.EQ.0.D0))THEN
			   ZBRENT=B
			   RETURN
			 ENDIF
			 IF((DABS(E).GE.TOL1).AND.(DABS(FA).GT.DABS(FB)))THEN
			   S=FB/FA
			   IF(A.EQ.C)THEN
			     P=2.D0*XM*S
				 Q=1.D0-S
			   ELSE
			     Q=FA/FC
				 R=FB/FC
				 P=S*(2.D0*XM*Q*(Q-R)-(B-A)*(R-1.D0))
				 Q=(Q-1.D0)*(R-1.D0)*(S-1.D0)
			   ENDIF
			   IF(P.GT.0) Q=-Q
			   P=DABS(P)
			   IF((2.D0*P).LT.MIN(3.D0*XM*Q-DABS(TOL1*Q),DABS(E*Q)))THEN
			     E=D
				 D=P/Q
			   ELSE
			     D=XM
				 E=D
			   ENDIF
			 ELSE
			   D=XM
			   E=D
			 ENDIF
			 A=B
			 FA=FB
			 IF(DABS(D).GT.TOL1)THEN
			   B=B+D
			 ELSE
			   B=B+DSIGN(TOL1,XM)
			 ENDIF
			 FB=FUNC(B,Z,UZ)
11         CONTINUE
           PAUSE '迭代次数过多'
		   ZBRENT=B
		   RETURN
		   END