MODULE define_DDEs IMPLICIT NONE INTEGER, PARAMETER :: NEQN=6,NLAGS=1,NEF=1 ! Physical parameters: DOUBLE PRECISION, PARAMETER :: D=5D0,c=0.5D0,G=1D0,n=100 ! Flag that indicates whether to evaluate the ODE model ! in subroutine DDES or the DDE model. LOGICAL :: useODEmodel CONTAINS SUBROUTINE DDES(T,Y,Z,DY) ! x = y(1), y = y(2), lambda = y(3), ! I_1 = y(4), I_2 = y(5), I_3 = y(6). DOUBLE PRECISION :: T DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY DOUBLE PRECISION, DIMENSION(NEQN,NLAGS) :: Z INTENT(IN) :: T,Y,Z INTENT(OUT) :: DY IF (useODEmodel) THEN ! Integrate the ODE model for T <= D: DY(1) = -Y(1)*Y(3) + G*Y(2) DY(2) = -DY(1) DY(4) = EXP(G*T)*Y(2) DY(5) = T*EXP(G*T)*Y(1)*Y(3) DY(6) = EXP(G*T)*Y(1)*Y(3) DY(3) = (c/n)*EXP(-G*T)*((DY(4)+DY(5))-G*(Y(4)+Y(5))) ELSE ! Integrate the DDE model for T >= D: DY(1) = -Y(1)*Y(3) + G*Y(2) DY(2) = -DY(1) DY(4) = EXP(G*T)*Y(2) - EXP(G*(T-D))*Z(2,1) DY(5) = D*EXP(G*T)*Y(1)*Y(3) - Y(6) DY(6) = EXP(G*T)*Y(1)*Y(3) - EXP(G*(T-D))*Z(1,1)*Z(3,1) DY(3) = (c/n)*EXP(-G*T)*((DY(4)+DY(5))-G*(Y(4)+Y(5))) ENDIF RETURN END SUBROUTINE DDES SUBROUTINE EF(T,Y,DY,Z,G) ! Event function to recognize when to change from the ! ODE model to the DDE model in subroutine DDES. 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 G(1) = T - D RETURN END SUBROUTINE EF SUBROUTINE CHNG(NEVENT,TEVENT,YEVENT,DYEVENT,HINIT, & DIRECTION,ISTERMINAL,QUIT) ! Function to change a flag so that the DDE model will ! be evaluated in subroutine DDES instead of the ODE model. INTEGER :: NEVENT INTEGER, DIMENSION(NEF) :: DIRECTION DOUBLE PRECISION :: TEVENT,HINIT DOUBLE PRECISION, DIMENSION(NEQN) :: YEVENT,DYEVENT LOGICAL :: QUIT LOGICAL, DIMENSION(NEF) :: ISTERMINAL INTENT(IN) :: NEVENT,TEVENT INTENT(INOUT) :: YEVENT,DYEVENT,HINIT,DIRECTION,& ISTERMINAL,QUIT IF (NEVENT == 1) THEN useODEmodel = .FALSE. END IF RETURN END SUBROUTINE CHNG END MODULE define_DDEs !****************************************************************** PROGRAM ex4p4p3 ! Example 4.4.3 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 EX4P4P3.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/) TYPE(DDE_SOL) :: SOL ! The fields of SOL are expressed in terms of the ! number of differential equations, NEQN, and the ! number of output points, NPTS: ! SOL%NPTS -- NPTS,number of output points. ! ! SOL%T(NPTS) -- values of independent variable, T. ! ! SOL%Y(NPTS,NEQN) -- values of dependent variable, Y, ! corresponding to values of SOL%T. ! When there is an event function, there are fields ! SOL%NE -- NE, number of events. ! SOL%TE(NE) -- locations of events ! SOL%YE(NE,NEQN) -- values of solution at events ! SOL%IE(NE) -- identifies which event occurred ! Physical parameters D,c,G,n are defined in the module ! define_DDEs as PARAMETERS. DOUBLE PRECISION, DIMENSION(NLAGS) :: LAGS=(/ D /) ! The simulation begins with an ODE model, so the ! constant HISTORY is just y(0). DOUBLE PRECISION, DIMENSION(NEQN) :: Y0= & (/ 0.8D0*n,0.2D0*n,0D0,0D0,0D0,0D0 /) DOUBLE PRECISION, DIMENSION(2) :: TSPAN= (/0D0, 4D0*D /) ! Local variables: INTEGER :: I,J ! Initialize flag that indicates whether to evaluate ! the ODE model in subroutine DDES or the DDE model. useODEmodel = .TRUE. SOL = DDE_SOLVER(NVAR,DDES,LAGS,Y0,TSPAN, & EVENT_FCN=EF,CHANGE_FCN=CHNG) ! Was the solver successful? IF (SOL%FLAG == 0) THEN ! Print integration statistics: CALL PRINT_STATS(SOL) ! Write the first three components of the solution ! to a file for subsequent plotting in Matlab. OPEN(UNIT=6, FILE='ex4p4p3.dat') DO I = 1,SOL%NPTS WRITE(UNIT=6,FMT='(4D12.4)') SOL%T(I),(SOL%Y(I,J),J=1,3) ENDDO PRINT *,' Normal return from DDE_SOLVER with results' PRINT *," written to the file 'ex4p4p3.dat'." PRINT *,' ' PRINT *,' These results can be accessed in Matlab' PRINT *,' and plotted by' PRINT *,' ' PRINT *," >> [t,y] = ex4p4p3;" PRINT *,' ' PRINT *,' ' PRINT *,' Kind of Event:' DO I = 1,SOL%NE IF(SOL%IE(I) == 1) THEN PRINT *,' Switched from ODEs to DDEs at t = ',SOL%TE(I) END IF END DO PRINT *,' ' ELSE PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',& SOL%FLAG ENDIF STOP END PROGRAM ex4p4p3