simulation.f

多学科优化软件isight培训教程初级pdf有很详细的例子讲解

***************************************************************************
*       10        20        30        40        50        60        70
*23456789012345678901234567890123456789012345678901234567890123456789012
*
       PROGRAM MAIN
       DOUBLE PRECISION Y(100),YP(100),T,TOUT,TEND,STPSZ
       DOUBLE PRECISION RELERR,ABSERR
       INTEGER NEQ,K,IFLAG
*      DOUBLE PRECISION T1,T2,T3(2),T4(2)
       DOUBLE PRECISION KK, CC, MM, FF
       COMMON /PARAMETERS/ KK, CC, MM, FF
       EXTERNAL F
       DATA T,Y,YP / 0.0D0,100*0.0D0,100*0.0D0 /
       DATA TEND,STPSZ,K,NEQ,RELERR,ABSERR / 2*0.0D0,2*0,2*0.0D0 /
       DATA KK, CC, MM, FF / 4*0.0D0 /
*
       OPEN(1,FILE='Inputs.txt')
       OPEN(2,FILE='Outputs.txt')
       read(1,*) T, K, STPSZ, TOUT, TEND, RELERR, ABSERR
       read(1,*)
       read(1,200) NEQ
  200  format(5x,i10)
       read(1,205) MM
  205  format(5x,1e15.6)
       read(1,210) KK
  210  format(5x,1e15.6)
       read(1,215) CC
  215  format(5x,1e15.6)
       read(1,220) FF
  220  format(5x,1e15.6)
       read(1,225) V0
  225  format(12x,1e15.6)
*      print*, T, K, STPSZ, TOUT, TEND, RELERR, ABSERR
*      print*,'       neq = ', NEQ
*      print*,'        MM = ', MM
*      print*,'        KK = ', KK
*      print*,'        CC = ',CC
*      print*,'        FF = ',FF
*      print*,'speed(t=0) = ',V0
*
* Apply initial conditions
*
       Y(1)  = V0
       Y(2)  = 0.0
       YP(1) = V0
       YP(2) = 0.0
*
* ETIME used for timing cpu execution time
*
*      T1 = ETIME(T3)
*
*      T = 0.0
*      K = 10
*      STPSZ = 0.01
*      TOUT = 0.1
*      TEND = 2.0
*
* IFLAG initializes the integrator
*
       IFLAG = 1
*
* tolerances for integration
*
*      RELERR = 0.00001D0
*      ABSERR = 0.00001D0
*
* write output at time zero
*
       WRITE(2,45) T,Y(1),Y(2)
  45   FORMAT(3E15.6)
*
   10  CONTINUE
*
* control output print frequency
*
       PRINT*,'iflag,tout', IFLAG,T
       TOUT = T + K*STPSZ
*
       CALL DE2(F,NEQ,Y,T,TOUT,RELERR,ABSERR,IFLAG)
*
* write output at each output time interval
*
       WRITE(2,45) T,Y(1),Y(2)
*
*
* If time is left in the simulation take the next time step.
*
       IF (T .GE. TEND) GOTO 100
       GOTO 10
*
* Determine the amount of simulation CPU time using ETIME.
*
* 100  T2 = ETIME(T4)
  100  CONTINUE
       PRINT*,'iflag,tout', IFLAG,T
*      PRINT*,'elapsed cpu time (sec):',T2-T1
*
*
       END
****************************************************************************
       SUBROUTINE F(T,Y,YP)
       DOUBLE PRECISION T,Y(100),YP(100)
       DOUBLE PRECISION PI
       DOUBLE PRECISION X1, X2, X1D
       DOUBLE PRECISION KK, CC, MM, FF
       COMMON /PARAMETERS/ KK, CC, MM, FF
       DATA PI / 3.14159265359D0 /
*
* Update working vectors with integrated state variables
*
       X2 = Y(1)              ! displacements
       X1 = Y(2)              ! velocity
*
* RHS
*
       X1D = -(KK*X2 + CC*X1)/MM + FF/MM
*
* initialize the state vector deriviative
*
       Y(1) = X2             ! translational position
       Y(2) = X1             ! translational velocity
*
       YP(1) = X1             ! translational velocity
       YP(2) = X1D            ! translational acceleration
*
       RETURN
*
       END
*
* de2_iprdct.f 
      SUBROUTINE DE2(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG)               DE     1
C                                                                       DE     2
C   SUBROUTINE  DE  INTEGRATES A SYSTEM OF UP TO 20 FIRST ORDER ORDINARYDE     3
C   DIFFERENTIAL EQUATIONS OF THE FORM                                  DE     4
C             DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN))                     DE     5
C             Y(I) GIVEN AT  T .                                        DE     6
C   THE SUBROUTINE INTEGRATES FROM  T  TO  TOUT .  ON RETURN THE        DE     7
C   PARAMETERS IN THE CALL LIST ARE INITIALIZED FOR CONTINUING THE      DE     8
C   INTEGRATION.  THE USER HAS ONLY TO DEFINE A NEW VALUE  TOUT         DE     9
C   AND CALL  DE  AGAIN.                                                DE    10
C                                                                       DE    11
C   DE  CALLS TWO CODES, THE INTEGRATOR  STEP  AND THE INTERPOLATION    DE    12
C   ROUTINE  INTRP .  STEP  USES A MODIFIED DIVIDED DIFFERENCE FORM OF  DE    13
C   THE ADAMS PECE FORMULAS AND LOCAL EXTRAPOLATION.  IT ADJUSTS THE    DE    14
C   ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR.  NORMALLY EACH CALL DE    15
C   TO  STEP  ADVANCES THE SOLUTION ONE STEP IN THE DIRECTION OF  TOUT .DE    16
C   FOR REASONS OF EFFICIENCY  DE  INTEGRATES BEYOND  TOUT  INTERNALLY, DE    17
C   THOUGH NEVER BEYOND T + 10*(TOUT-T), AND CALLS  INTRP  TO           DE    18
C   INTERPOLATE THE SOLUTION AT  TOUT .  AN OPTION IS PROVIDED TO STOP  DE    19
C   THE INTEGRATION AT  TOUT  BUT IT SHOULD BE USED ONLY IF IT IS       DE    20
C   IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND  TOUT .               DE    21
C                                                                       DE    22
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,       DE    23
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS:  THE INITIAL  DE    24
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.                  DE    25
C                                                                       DE    26
C   THE PARAMETERS FOR  DE  ARE:                                        DE    27
C      F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT DE    28
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     DE    29
C      Y(*) -- SOLUTION VECTOR AT  T                                    DE    30
C      T -- INDEPENDENT VARIABLE                                        DE    31
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED                       DE    32
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCALDE    33
C           ERROR TEST.  AT EACH STEP THE CODE REQUIRES                 DE    34
C             ABS(LOCAL ERROR) .LE. ABS(Y)*RELERR + ABSERR              DE    35
C           FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS  DE    36
C      IFLAG -- INDICATES STATUS OF INTEGRATION                         DE    37
C                                                                       DE    38
C   FIRST CALL TO DE --                                                 DE    39
C                                                                       DE    40
C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAY  DE    41
C   IN THE CALL LIST, Y(NEQN), DECLARE  F  IN AN EXTERNAL STATEMENT,    DE    42
C   SUPPLY THE SUBROUTINE F(T,Y,YP) TO EVALUATE                         DE    43
C      DY(I)/DT = YP(I) = F(T,Y(1),Y(2),...,Y(NEQN))                    DE    44
C   AND INITIALIZE THE PARAMETERS:                                      DE    45
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED                     DE    46
C      Y(*) -- VECTOR OF INITIAL CONDITIONS                             DE    47
C      T -- STARTING POINT OF INTEGRATION                               DE    48
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED                       DE    49