Fortran Coder

查看: 91|回复: 7

[求助] 寻找极点

[复制链接]

22

帖子

10

主题

0

精华

熟手

F 币
121 元
贡献
71 点
发表于 2019-2-27 17:25:38 | 显示全部楼层 |阅读模式
QQ123截图20190227093804.png 调用子程序求极点,现在极点找不到,yourf2v是这个函数。我现在的工作就是把这个函数换成我给的函数,但是结果不对,不知道问题出在哪里了?
回复

使用道具 举报

1359

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂伙计

F 币
407 元
贡献
365 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2019-2-27 18:18:55 | 显示全部楼层
请复制粘贴代码,给全(包含 yourf2v 和 minsimplex2)不要截图。(因为截图没法复制粘贴)

22

帖子

10

主题

0

精华

熟手

F 币
121 元
贡献
71 点
 楼主| 发表于 2019-2-28 09:05:09 | 显示全部楼层
[Fortran] 纯文本查看 复制代码
program polesearching
        implicit none
!!! two variables
        integer Npar,maxiter,inter
        parameter(Npar=2)
        real*8 par(Npar),del(Npar),error,fvalue,yourf2v
!! one variable to search zero
        real*8 yourf1v,x,lower,upper,x0,rr,repsdzero
        integer mxfdzero
        external yourf1v
!! two variable for subroutine: nelmin
        real*8 xmin(Npar),step(Npar),reqmin
        integer konvge,kcount,icount,numres,ifault
!! two variable for subroutine: MINIMIZE_2D
        real*8 par0(Npar)
        integer info
        external subyourf2v
      maxiter=1500
        inter=0
        error= 1.d-8
!!
        par(1)= 1011.0d0
        par(2)= 0.0d0
        del(1)= 0.001d0
        del(2)= 0.001d0
        print*, yourf2v(par)
        call minsimplex2(Npar,par,del,error,maxiter,fvalue,yourf2v,inter)
        write(*,'(A10,2F20.10,E12.2)') 'solution', par(1),par(2),fvalue
        end program polesearching
        function yourf2v(par)
        implicit none
        real*8 yourf2v
        integer Npar
        parameter(Npar=2)
        real*8 par(Npar),x,y
        complex*16 ui,s,E
        real*8 mp,mq,mpi,mk,meta,metap,Fpi,mu
      real*8 ct,st
        real*8 asL10n1, asL10n2, asL10n3
      integer nchan
      parameter (nchan=3)
        complex*16 gij10(nchan,nchan),gij102(nchan,nchan),
     -tijgij10(nchan,nchan),ro10(nchan,nchan),
     -i33tijgij10(nchan,nchan),ss(nchan,nchan),
     -pwTij10TollerM33MLo(nchan,nchan),s01(nchan,nchan),s11(nchan,nchan)
     -, EtauNijMlo10(nchan,nchan),DeltaLoguNijMlo10(nchan,nchan),
     -EtauNijMlo1d2(nchan,nchan),DeltaLoguNijMlo1d2(nchan,nchan)
        integer info    !qiu inverse  matrix
        integer ipvt(nchan)  !qiu inverse  matrix
      complex*16  work(nchan),detin(2),deterin  !qiu inverse  matrix
        complex*16 pwTij10p2M33MAnl(nchan,nchan)
        complex*16 i33umat(nchan,nchan)
        real*8 pi
        common /pi/ pi
      real*8 asl10
        common /lec10/ asl10
        common /physicalmass/mpi,mk,meta,metap,fpi,mu
      complex*16 GfdrI
        external GfdrI
        pi= dacos(-1.0d0)
      mpi=137.3d0
        mk=495.6d0
        meta=547.9d0
        metap=957.7d0
        Fpi =92.4d0
        mu=770.0d0
        asL10 = -1.44d0
        asL10n1 = asl10
      asL10n2 = asL10
   asL10n3 = asL10
st=-0.279d     !theta is -16.2
        ct=0.96d0    !theta is -16.2
      ui=(0.d0,1.d0)
      E=par(1)+ui*par(2)
        s=E**2
        gij10(1,1)=GfdrI(asL10n1,s,mpi**2,meta**2)
      gij10(1,2)=(0.0d0,0.0d0)
        gij10(1,3)=(0.0d0,0.0d0)
        gij10(2,1)=(0.0d0,0.0d0)
        gij10(2,2)=GfdrI(asL10n2,s,mk**2,mk**2)
        gij10(2,3)=(0.0d0,0.0d0)
        gij10(3,1)=(0.0d0,0.0d0)
        gij10(3,2)=(0.0d0,0.0d0)
        gij10(3,3)=GfdrI(asL10n3,s,mpi**2,metaP**2)
      pwTij10p2M33MAnl(1,1)=(mpi**2*(ct**2 - 2.0*Sqrt(2.0)*ct*st
     - + 2.0*st**2))/(3.0*Fpi**2)
      pwTij10p2M33MAnl(1,2)=  (ct*(3.0*meta**2 + 8.0*mk**2 + mpi**2
     - - 9.0*s) + 2.0*Sqrt(2.0)*(2.0*mk**2 + mpi**2)*st)/
     -(6.0*Sqrt(6.0)*Fpi**2)
        pwTij10p2M33MAnl(1,3)=(mpi**2*(Sqrt(2.0)*ct**2 -
     -ct*st - Sqrt(2.0)*st**2))/(3.0*Fpi**2)
        pwTij10p2M33MAnl(2,1)=(ct*(3.0*meta**2 + 
     -8.0*mk**2 + mpi**2 - 9.0*s) + 
     - 2.0*Sqrt(2.0)*(2.0*mk**2 + mpi**2)*st)/(6.0*Sqrt(6.0)*Fpi**2)
        pwTij10p2M33MAnl(2,2)=s/(4.0*Fpi**2)
        pwTij10p2M33MAnl(2,3)=   -(2.0*Sqrt(2.0)*ct*(2.0*mk**2 + mpi**2) - 
     -     (3.0*metaP**2 + 8.0*mk**2 + mpi**2 - 9.0*s)*st)/
     -  (6.0*Sqrt(6.0)*Fpi**2)
        pwTij10p2M33MAnl(3,1)=(mpi**2*(Sqrt(2.0)*ct**2 -
     -ct*st - Sqrt(2.0)*st**2))/(3.0*Fpi**2)
        pwTij10p2M33MAnl(3,2)=  -(2.0*Sqrt(2.0)*ct*(2.0*mk**2 + mpi**2) - 
     -     (3.0*metaP**2 + 8.0*mk**2 + mpi**2 - 9.0*s)*st)/
     -  (6.0*Sqrt(6.0)*Fpi**2)
        pwTij10p2M33MAnl(3,3)=(mpi**2*(2.0*ct**2 + 2.0*Sqrt(2.0)*ct*st +
     -st**2))/(3.0*Fpi**2)
      i33umat(1,1)=(1.0d0,0.0d0)
        i33umat(1,2)=(0.0d0,0.0d0)
        i33umat(1,3)=(0.0d0,0.0d0)
        i33umat(2,1)=(0.0d0,0.0d0)
        i33umat(2,2)=(1.0d0,0.0d0)
        i33umat(2,3)=(0.0d0,0.0d0)
        i33umat(3,1)=(0.0d0,0.0d0)
        i33umat(3,2)=(0.0d0,0.0d0)
        i33umat(3,3)=(1.0d0,0.0d0)
      tijgij10=MATMUL(pwTij10p2M33MAnl,gij10)
      i33tijgij10=i33umat- tijgij10
      call zgefa( i33tijgij10,nchan,nchan,ipvt,info)    !! call both
      call zgedi(i33tijgij10,nchan,nchan,ipvt,detin,work,11) !! call  both 
        deterin=detin(1)*10.**detin(2) !! determinant
        print*, deterin
        pwTij10TollerM33MLo=MATMUL(i33tijgij10,pwTij10p2M33MAnl) !no problem
        yourf2v=deterin
        end
SUBROUTINE minsimplex2 
     1   (NVAR,PUNT,DELT,ERROR,MAXITER,VALOR,FUNCIO,INTER)
        IMPLICIT REAL*8 (a-h,o-z)
        PARAMETER(nmax=30)
        DIMENSION PUNT(nvar),DELT(nvar)
        DIMENSION Y(nmax+1), P(nmax,nmax+1)
        common /iierror/ierror !(mio)
        EXTERNAL FUNCIO
        DO j=1,nvar+1
        DO i=1, nvar
        p(i,j)=PUNT(i)
        ENDDO
        IF(j .GT. 1) THEN
        p(j-1,j)=punt(j-1)+delt(j-1)
        ENDIF
        ENDDO
10      CONTINUE
        CALL SIMPLEX2 
     1  (P,Y,nvar,nmax,error,ierror,maxiter,funcio,millor,ipitjor)
        IF (INTER .EQ. 1) THEN
        WRITE(*,*)
        WRITE(*,*)' ============================================='
        WRITE(*,'(1X,3A20)') 'Par. Inic.','Par. Obtingu.', 'ERROR'
        DO i=1,NVAR
        WRITE(*,'(1X, 3F20.8)') PUNT(i),P(i,millor), 
     1   ABS(P(i,millor)-P(i,ipitjor))/2
        ENDDO
        WRITE(*,'(1X,A13,1P,E20.8,5x,a13,E20.8)') 
     1   ' FUNC. MILLOR', y(millor), ' FUNC. PITJOR',y(ipitjor)
        IF (ierror .EQ. 0) THEN
        WRITE(*,*) ' SIMPLEX arrib?a converg妌cia'
        WRITE(*,*) ' Entreu 1 -> dividir error per 10, 0 -> finir'
        READ(*,*) icas
        IF (icas .EQ. 1) THEN
        error=error/10
        WRITE(*,*) ' Nou error=',error
        GOTO 10
        ENDIF
        ELSE
        WRITE(*,*) ' No s''ha aconseguit la converg妌cia'
        WRITE(*,*) ' Entreu 1 per ITERAR, 0 per  FINIR'
        READ(*,*) icas
        IF (icas .EQ. 1) GOTO 10
        ENDIF
        ENDIF
c       Actualitzar DELTA i PUNT
        valor=Y(millor)
        DO i=1, nvar
        PUNT(i)=P(i,millor)
        vmin=P(i,1)
c        print *,'kkk'
        vmax=P(i,1)
        DO j=2,nvar+1
        vmin=MIN(vmin,p(i,j))
        vmax=MAX(vmax,p(i,j))
        ENDDO
        DELT(I)=vmax-vmin
        ENDDO
        RETURN
        END
        SUBROUTINE simplex2 
     1  (P,Y,NVAR,NDIM,ERROR,IERROR,MAXITER,FMV,MILLOR,PITJOR)
        IMPLICIT REAL*8 (a-h, o-z)
        PARAMETER(nmax=30, unitat=1.d0)
        DIMENSION p(ndim,ndim+1), y(ndim+1), iordre(nmax+1)
        DIMENSION cdg(nmax+1), pn(nmax+1), pp(nmax+1)
        INTEGER millor, pitjor , quasipitjor
        common /espolo/iespolo !(mio)
        IF (nvar .GT. nmax) THEN
        WRITE(*,*) ' Mass variables in SIMPLEX'
        STOP
        ENDIF
        IERROR=0
        nvert=nvar+1        
        DO 1000 iter=1, maxiter
        CALL ordenar2(y,iordre,nvert)
        millor=iordre(1)
        pitjor=iordre(nvert)
        quasipitjor=iordre(nvert-1)
        den = ABS(y(pitjor)+y(millor))/2.d0
c(mio)        IF ((y(pitjor)-y(millor)) .LT. error*den) RETURN
        aakk=abs(y(pitjor)/y(millor))-1.d0 !(mio)
        errorkk=1.d-6 !(mio)
        IF (aakk.LT.errorkk) then
        iespolo=0
        return
        endif
c        print *,'a',iter,y(millor),y(pitjor),y(pitjor)/y(millor)
c     +          ,(y(pitjor)-y(millor))/den
        if (y(millor).LT.error) then
        iespolo=1
        return
        endif
        IF (iter .EQ. maxiter) GOTO 1000
c       Moviments ... Primer busque CDG
        DO i=1,nvar
        cdg(i)=0
        DO j=1, nvert
        IF (j .NE. pitjor) cdg(i)=cdg(i)+P(i,j)
        ENDDO
        cdg(i)=cdg(i)/nvar
        ENDDO
        CALL combinacio2(-1.d0, P(1,pitjor), 2.d0, CDG, pn, nvar)
        F1=fmv(pn)
        F2=fmv(pp)
        IF (f1 .GT. f2) THEN
        y(pitjor)=f2
        DO i=1,nvar
        p(i,pitjor)=pp(i)
        ENDDO
        GOTO 1000
        ENDIF
        GOTO 2000
40      if (F1 .lt .Y(QUASIPITJOR)) goto 2000
        IF (f1 .GT. y(pitjor)) GOTO 3000
4000    y(pitjor)=f1
        DO i=1,nvar
        p(i,pitjor)=pn(i)
        ENDDO
ccc (m) hay problemas si pitjor es igual a cdg:
             akk1=abs(P(1,pitjor)/cdg(1)-1.d0)
        akk2=abs(P(2,pitjor)/cdg(2)-1.d0)
        if ((akk1.LT.1.d-6).and.(akk2.LT.1.d-6)) then
             iespolo=0
        return
        endif
c       Contracci?del pitjor punt (interior o exterior)
3000    CALL combinacio2(0.5d0, P(1,pitjor), 0.5d0, CDG, pn, nvar)
        F1=fmv(pn)
c       Acceptable?
        IF (f1 .LT. y(quasipitjor)) GOTO 2000
c       No, per?si 俿 MILLOR que PITJOR l'accepte i contrac de nou
        IF (F1 .LT. y(pitjor)) GOTO 4000

c       CAS 4: Contracci?de tot el splex
        DO j=1,nvert
        IF (j .NE. millor) THEN
        CALL combinacio2(0.5d0,P(1,j),0.5d0,P(1,millor),P(1,j),nvar)
        Y(j)=fmv(pn)
        ENDIF
        ENDDO
        goto 1000
c       Acceptar pn i comen嘺r de nou
2000    CONTINUE
        y(pitjor)=f1
        DO i=1,nvar
        p(i,pitjor)=pn(i)
        ENDDO
1000    CONTINUE
c       S'arriba ac?si no convergeix en MAXITER
c!(mio)        WRITE(*,*) ' Massa iteracions en SIMPLEX'
        ierror=1
        RETURN
        END
!!!!!!!!!!!!!!!!!
        SUBROUTINE combinacio2(a, P, b, Q, R, n)
        IMPLICIT REAL*8 (a-h, o-z)
        DIMENSION p(n), q(n), r(n)
        DO 10 i=1, n
10      r(i)=a*p(i)+b*q(i)
        RETURN
        END
!!!!!!!!!!!!!!!
        SUBROUTINE ordenar2(f,iordre,n)
        IMPLICIT REAL*8 (a-h,o-z)
        DIMENSION F(n), IORDRE(n)
        iordre(1)=1
        DO 10 j=2,n
        x=f(j)
        DO 20 i=j-1,1,-1
        IF (x .GT. f(iordre(i)) ) GOTO 10
20      iordre(i+1)=iordre(i)
        i=0
10      iordre(i+1)=j
        RETURN
        END
        

1359

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂伙计

F 币
407 元
贡献
365 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2019-2-28 09:44:55 | 显示全部楼层
缺少 GFDRI 函数,ZGEFA 和 ZGEDI 函数。

请先确保自己的代码能通过编译。

22

帖子

10

主题

0

精华

熟手

F 币
121 元
贡献
71 点
 楼主| 发表于 2019-2-28 10:17:43 | 显示全部楼层
fcode 发表于 2019-2-28 09:44
缺少 GFDRI 函数,ZGEFA 和 ZGEDI 函数。

请先确保自己的代码能通过编译。

确实是,因为他限制字节。现在就是调用一个求极点的子程序,需要改变得是这个函数就行。

1359

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂伙计

F 币
407 元
贡献
365 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2019-2-28 10:34:33 | 显示全部楼层
你可以上传附件。

22

帖子

10

主题

0

精华

熟手

F 币
121 元
贡献
71 点
 楼主| 发表于 2019-2-28 16:01:43 | 显示全部楼层

寻找极点

polesearching.for (18.51 KB, 下载次数: 4)

22

帖子

10

主题

0

精华

熟手

F 币
121 元
贡献
71 点
 楼主| 发表于 2019-2-28 16:52:12 | 显示全部楼层
问题已解决谢谢大家的关注
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2019-3-25 07:08

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表