MODULE define_DDEs IMPLICIT NONE INTEGER, PARAMETER :: NEQN=1,NLAGS=1,NEF=2 ! Physical parameters: DOUBLE PRECISION, PARAMETER :: R=3.5D0,M=19D0 CONTAINS SUBROUTINE DDES(T,Y,Z,DY) DOUBLE PRECISION :: T DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Z INTENT(IN) :: T,Y,Z INTENT(OUT) :: DY DY(1) = R*Y(1)*(1D0 - Z(1,1)/M) RETURN END SUBROUTINE DDES SUBROUTINE HISTORY(T,Y) DOUBLE PRECISION :: T DOUBLE PRECISION, DIMENSION(NEQN) :: Y INTENT(IN) :: T INTENT(OUT) :: Y ! Use the JUMPS option to tell solver about ! the discontinuity at t = 0. IF (T == 0D0) THEN Y(1) = 19.001D0 ELSE Y(1) = 19D0 END IF RETURN END SUBROUTINE HISTORY SUBROUTINE EF(T,Y,DY,Z,G) DOUBLE PRECISION :: T DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Z DOUBLE PRECISION, DIMENSION(NEF) :: G INTENT(IN) :: T,Y,DY,Z INTENT(OUT) :: G ! Locate extrema as points where the derivative ! vanishes. The DIRECTION option is used to ! distinguish maxima and minima. G = (/ DY(1),DY(1) /) RETURN END SUBROUTINE EF END MODULE define_DDEs !****************************************************************** PROGRAM ex4p4p4 ! Example 4.4.4 in the SGT book. ! ! The DDE is defined in the module define_DDEs. The problem ! is solved here with DDE_SOLVER and its output written to a ! file. The auxilary function EX4P4P4.M imports the data into ! Matlab and plots it. USE define_DDEs USE DDE_SOLVER_M IMPLICIT NONE ! The quantities ! ! NEQN = number of equations ! NLAGS = number of delays ! NEF = number of event functions ! ! are defined in the module define_DDEs as PARAMETERs so ! they can be used for dimensioning arrays here. They are ! passed to the solver in the array NVAR. INTEGER, DIMENSION(3) :: NVAR = (/NEQN,NLAGS,NEF/) DOUBLE PRECISION, DIMENSION(1) :: DELAY = (/ 0.74D0 /) TYPE(DDE_SOL) :: SOL TYPE(DDE_INT) :: YINT TYPE(DDE_OPTS) :: OPTS DOUBLE PRECISION, PARAMETER :: T0=0D0,TFINAL=40D0 ! Prepare output points INTEGER, PARAMETER :: NOUT=1000 DOUBLE PRECISION, DIMENSION(NOUT) :: TINT ! Local variables: INTEGER :: I OPTS = DDE_SET(RE=1D-4,AE=1D-7,INTERPOLATION=.TRUE.,& JUMPS=(/ 0D0 /),DIRECTION=(/ +1,-1 /)) SOL = DDE_SOLVER(NVAR,DDES,DELAY,HISTORY, & TSPAN=(/ T0,TFINAL /),OPTIONS=OPTS,EVENT_FCN=EF) ! Was the solver successful? IF (SOL%FLAG == 0) THEN ! Form values for phase plane plot: TINT = (/ (T0+(I-1)*((TFINAL-T0)/(NOUT-1)), I=1,NOUT) /) YINT = DDE_VAL(TINT,SOL,DERIVATIVES=.TRUE.) ! Write the solution to a file for subsequent plotting ! in Matlab. OPEN(UNIT=6, FILE='ex4p4p4.dat') DO I = 1,NOUT WRITE(UNIT=6,FMT='(3D12.4)') TINT(I),YINT%YT(I,1),& YINT%DT(I,1) END DO ! Write the extrema to a file for plotting. OPEN(UNIT=7,FILE='ex4p4p4extr.dat') DO I = 1,SOL%NE WRITE(UNIT=7,FMT='(I10,2D12.4)') SOL%IE(I), & SOL%TE(I),SOL%YE(I,1) END DO PRINT *,' Normal return from DDE_SOLVER with results' PRINT *," written to the file 'ex4p4p4.dat' and the" PRINT *," extrema written to 'ex4p4p4extr.dat'." PRINT *,' ' PRINT *,' These results can be accessed in Matlab' PRINT *,' and plotted by' PRINT *,' ' PRINT *," >> [t,y,yp,te,ye,ie] = ex4p4p4;" PRINT *,' ' ELSE PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',& SOL%FLAG ENDIF STOP END PROGRAM ex4p4p4