odeint1.f90

FORTRANvisualfortran常用数值算法集及源码

SUBROUTINE ODEINT1(YSTART,NVAR,X1,X2,EPS,H1,HMIN,NOK
     ,NBAD,DERIVS,BSSTEP)
PARAMETER (MAXSTP=10000,NMAX=10,TWO=2.0,ZERO=0.0,
           TINY=1.E-30)
COMMON /PATH/ KMAX,KOUNT,DXSAV,XP(200),YP(10,200)
DIMENSION YSTART(NVAR),YSCAL(NMAX),Y(NMAX),DYDX(NMAX)
X=X1
H=SIGN(H1,X2-X1)
NOK=0
NBAD=0
KOUNT=0
DO 11 I=1,NVAR
    Y(I)=YSTART(I)
CONTINUE
XSAV=X-DXSAV*TWO
DO 16 NSTP=1,MAXSTP
    CALL DERIVS(X,Y,DYDX)
    DO 12 I=1,NVAR
        YSCAL(I)=ABS(Y(I))+ABS(H*DYDX(I))+TINY
    CONTINUE
    IF(KMAX.GT.0)THEN
        IF(ABS(X-XSAV).GT.ABS(DXSAV))  THEN
            IF(KOUNT.LT.KMAX-1)  THEN
                KOUNT=KOUNT+1
                XP(KOUNT)=X
                DO 13 I=1,NVAR
                    YP(I,KOUNT)=Y(I)
                CONTINUE
                XSAV=X
            ENDIF
        ENDIF
    ENDIF
    IF((X+H-X2)*(X+H-X1).GT.ZERO) H=X2-X
    CALL BSSTEP(Y,DYDX,NVAR,X,H,EPS,YSCAL,HDID,
                HNEXT,DERIVS)
    IF(HDID.EQ.H)  THEN
        NOK=NOK+1
    ELSE
        NBAD=NBAD+1
    ENDIF
    IF((X-X2)*(X2-X1).GE.ZERO) THEN
        DO 14 I=1,NVAR
            YSTART(I)=Y(I)
        CONTINUE
        IF(KMAX.NE.0) THEN
            KOUNT=KOUNT+1
            XP(KOUNT)=X
            DO 15 I=1,NVAR
                YP(I,KOUNT)=Y(I)
            CONTINUE
        ENDIF
        RETURN
    ENDIF
    IF(ABS(HNEXT).LT.HMIN)
             PAUSE 'Stepsize smaller than minimum.'
    H=HNEXT
CONTINUE
PAUSE 'Too many steps.'
RETURN
END