| 查看: 424 | 回复: 2 | |||
| 当前主题已经存档。 | |||
329543477木虫 (正式写手)
|
[交流]
【求助】帮忙修改用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 ] |
回复此楼» 猜你喜欢
【修正版】单环用磁约束低速电子实验方案（简化版）
已经有0人回复
-
桂林理工大学物理学专业招收调剂，还有三个名额！！！
已经有22人回复
-
物理学I论文润色/翻译怎么收费?
已经有179人回复
-
考博自荐
已经有4人回复
-
山东大学第二批博士研究生招生
已经有0人回复
-
中国科学院东莞材料科学与技术研究所-2026年博士招生-吴昊研究员-磁学与自旋电子学
已经有0人回复
-
《电磁学》教材推荐
已经有1人回复
-
【急招】合肥工大核聚变材料计算方向2026级工程博士生
已经有4人回复
-
大豆异黄酮分离
已经有0人回复
-
湖南大学材料学院急招2026年博士生，临时增加一名博士联培指标
已经有10人回复
-
天津理工大学晶体材料全国重点实验室刘红军教授课题组招收博士生1-2名
已经有0人回复
329543477
木虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 2073.9
- 散金: 136
- 红花: 10
- 沙发: 1
- 帖子: 607
- 在线: 101.1小时
- 虫号: 668790
- 注册: 2008-12-05
- 性别: MM
- 专业: 化学反应工程
2楼2009-09-01 16:59:39
329543477
木虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 2073.9
- 散金: 136
- 红花: 10
- 沙发: 1
- 帖子: 607
- 在线: 101.1小时
- 虫号: 668790
- 注册: 2008-12-05
- 性别: MM
- 专业: 化学反应工程
3楼2009-09-05 19:59:37
20