24小时热门版块排行榜

返回列表

当前主题已经存档。

329543477

木虫 (正式写手)

应助: 0 (幼儿园)
金币: 2073.9
散金: 136
红花: 10
沙发: 1
帖子: 607
在线: 101.1小时
虫号: 668790
注册: 2008-12-05
性别: MM
专业: 化学反应工程

[交流] 【求助】帮忙修改用powell方法解题的fortran程序

利用所给的powell方法求解如下函数的极小值点：f(x1,x2)=(x1+x2)**2+(x1-1)**2
初始点为列矩阵（2，1），搜索方法为（1，0）（0，1）
powell方法的程序如下：
SUBROUTINE POWELL(P,XI,N,NP,FTOL,ITER,FRET)
   PARAMETER (NMAX=20,ITMAX=200)
   DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX)
   FRET=FUNC(P)
   DO 11 J=1,N
      PT(J)=P(J)
11 CONTINUE
   ITER=0
1    ITER=ITER+1
   FP=FRET
   IBIG=0
   DEL=0.
   DO 13 I=1,N
      DO 12 J=1,N
      XIT(J)=XI(J,I)
12    CONTINUE
      CALL LINMIN(P,XIT,N,FRET)
      IF(ABS(FP-FRET).GT.DEL)THEN
      DEL=ABS(FP-FRET)
      IBIG=I
      ENDIF
13 CONTINUE
   IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN
   IF(ITER.EQ.ITMAX) PAUSE 'Powell exceeding maximum iterations.'
   DO 14 J=1,N
      PTT(J)=2.*P(J)-PT(J)
      XIT(J)=P(J)-PT(J)
      PT(J)=P(J)
14 CONTINUE
   FPTT=FUNC(PTT)
   IF(FPTT.GE.FP)GO TO 1
   T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2
   IF(T.GE.0.)GO TO 1
   CALL LINMIN(P,XIT,N,FRET)
   DO 15 J=1,N
      XI(J,IBIG)=XIT(J)
15 CONTINUE
   GO TO 1
   END

   FUNCTION BRENT(AX,BX,CX,F,TOL,XMIN)
   PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
   A=MIN(AX,CX)
   B=MAX(AX,CX)
   V=BX
   W=V
   X=V
   E=0.
   FX=F(X)
   FV=FX
   FW=FX
   DO 11 ITER=1,ITMAX
      XM=0.5*(A+B)
      TOL1=TOL*ABS(X)+ZEPS
      TOL2=2.*TOL1
      IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3
      IF(ABS(E).GT.TOL1) THEN
      R=(X-W)*(FX-FV)
      Q=(X-V)*(FX-FW)
      P=(X-V)*Q-(X-W)*R
      Q=2.*(Q-R)
      IF(Q.GT.0.) P=-P
      Q=ABS(Q)
      ETEMP=E
      E=D
      IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR.
   *       P.GE.Q*(B-X)) GOTO 1
      D=P/Q
      U=X+D
      IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X)
      GOTO 2
      ENDIF
1    IF(X.GE.XM) THEN
      E=A-X
      ELSE
      E=B-X
      ENDIF
      D=CGOLD*E
2    IF(ABS(D).GE.TOL1) THEN
      U=X+D
      ELSE
      U=X+SIGN(TOL1,D)
      ENDIF
      FU=F(U)
      IF(FU.LE.FX) THEN
      IF(U.GE.X) THEN
         A=X
      ELSE
         B=X
      ENDIF
      V=W
      FV=FW
      W=X
      FW=FX
      X=U
      FX=FU
      ELSE
      IF(U.LT.X) THEN
         A=U
      ELSE
         B=U
      ENDIF
      IF(FU.LE.FW .OR. W.EQ.X) THEN
         V=W
         FV=FW
         W=U
         FW=FU
      ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
         V=U
         FV=FU
      ENDIF
      ENDIF
11 CONTINUE
   PAUSE 'Brent exceed maximum iterations.'
3    XMIN=X
   BRENT=FX
   RETURN
   END

   FUNCTION F1DIM(X)
   PARAMETER (NMAX=50)
   COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
   DIMENSION XT(NMAX)
   DO 11 J=1,NCOM
      XT(J)=PCOM(J)+X*XICOM(J)
11 CONTINUE
   F1DIM=FUNC(XT)
   RETURN
   END

   SUBROUTINE LINMIN(P,XI,N,FRET)
   PARAMETER (NMAX=50,TOL=1.E-4)
   EXTERNAL F1DIM
   DIMENSION P(N),XI(N)
   COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
   NCOM=N
   DO 11 J=1,N
      PCOM(J)=P(J)
      XICOM(J)=XI(J)
11 CONTINUE
   AX=0.
   XX=1.
   BX=2.
   CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM)
   FRET=BRENT(AX,XX,BX,F1DIM,TOL,XMIN)
   DO 12 J=1,N
      XI(J)=XMIN*XI(J)
      P(J)=P(J)+XI(J)
12 CONTINUE
   RETURN
   END

   SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNC)
   PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20)
   FA=FUNC(AX)
   FB=FUNC(BX)
   IF(FB.GT.FA)THEN
      DUM=AX
      AX=BX
      BX=DUM
      DUM=FB
      FB=FA
      FA=DUM
   ENDIF
   CX=BX+GOLD*(BX-AX)
   FC=FUNC(CX)
1    IF(FB.GE.FC)THEN
      R=(BX-AX)*(FB-FC)
      Q=(BX-CX)*(FB-FA)
      U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
      ULIM=BX+GLIMIT*(CX-BX)
      IF((BX-U)*(U-CX).GT.0.)THEN
      FU=FUNC(U)
      IF(FU.LT.FC)THEN
         AX=BX
         FA=FB
         BX=U
         FB=FU
         GO TO 1
      ELSE IF(FU.GT.FB)THEN
         CX=U
         FC=FU
         GO TO 1
      ENDIF
      U=CX+GOLD*(CX-BX)
      FU=FUNC(U)
      ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN
      FU=FUNC(U)
      IF(FU.LT.FC)THEN
         BX=CX
         CX=U
         U=CX+GOLD*(CX-BX)
         FB=FC
         FC=FU
         FU=FUNC(U)
      ENDIF
      ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN
      U=ULIM
      FU=FUNC(U)
      ELSE
      U=CX+GOLD*(CX-BX)
      FU=FUNC(U)
      ENDIF
      AX=BX
      BX=CX
      CX=U
      FA=FB
      FB=FC
      FC=FU
      GO TO 1
   ENDIF
   RETURN
   END

自己尝试编写的球目标函数的程序如下：
    parameter(ndim=2,ftol=1.e-6)
real xi(ndim,ndim),pp(ndim)
data xi/1,0,0,1/
data pp/2.,1./
np=ndim
call powell(pp,xi,ndim,np,ftol,iter,fret)
print *,pp
print *,'fret=',fret
end

   function f(x1,x2)
real x1,x2,f(x1,x2)
   f(x1,x2)=(x1+x2)**2+(x1-1)**2
   end

求助高手指点，程序改怎么修改？

[ Last edited by dongdong3881 on 2009-12-12 at 09:24 ]

回复此楼