The following source text for the main version of the computer program, which integrated the Lorenz model, is written in Fortran 77 and uses the standard library Runge-Kutta order 4 method. The source text is commented. C------------------------------------------------------------------------- PROGRAM LORENZ IMPLICIT NONE REAL*8 So(3), S(3), EST1(3), EST2(3), EST(3), DSDT1(9), DSDT2(9) REAL*8 H, H_NEW, T, T_COPY, T_UP, T_END REAL*8 TOLERANCE, ERROR REAL*8 SIGMA, R, B, AX, AY, AZ !parameter: Ai = forcing REAL*8 COUNTER, T_PRINT, NR_DATA INTEGER STP, STP_UP CHARACTER LOG1, LOG2, LOG3, LOG4, LOG5 LOGICAL LOG_A COMMON SIGMA, R, B, AX, AY, AZ EXTERNAL RHS EXTERNAL DRKSTP C------------------------ INPUT : ---------------------------------- WRITE(*,*)'__________________________________________' WRITE(*,*)'ENTER TOLERANCE FOR ADAPTIVE STEPSIZE :' READ(*,*) TOLERANCE WRITE(*,*)'ENTER OVERALL INTEGRATION TIME ( T_END ) :' READ(*,*) T_END WRITE(*,*)'ENTER TIME ( T_UP ) UP WHICH DATA SHALL' WRITE(*,*)'BE PRINTED TO THE DATA FILE(S) :' READ(*,*) T_UP WRITE(*,*)'ENTER NUMBER ( NR_DATA ) OF DATA POINTS' WRITE(*,*)'PRINTED TO FILE :' READ(*,*) NR_DATA WRITE(*,*)'ENTER INITIAL Xo :' READ(*,*) S(1) WRITE(*,*)'ENTER INITIAL Yo :' READ(*,*) S(2) WRITE(*,*)'ENTER INITIAL Zo :' READ(*,*) S(3) C WRITE(*,*)'ENTER PARAMETER "SIGMA" :' C READ(*,*) SIGMA WRITE(*,*)'ENTER PARAMETER "R" :' READ(*,*) R C WRITE(*,*)'ENTER PARAMETER "B" :' C READ(*,*) B WRITE(*,*)'ENTER AX, AY, AZ, ( I.E. FORCING) :' READ(*,*) AX READ(*,*) AY READ(*,*) AZ SIGMA = 10 B = 2.666666666666666666 C-------------------- SOME VALUES : --------------------- H = 1.0d-10 T = 0.0d0 STP = 0 T_PRINT = ( (T_END - T_UP) / NR_DATA ) COUNTER = 0.0D0 LOG_A = .TRUE. So = S C----------- DATA TO FILE TO BE CHOSEN : ------------------ WRITE(*,*)'PRINT X-Y DATA TO FILE "xy.dat" ? :' READ(*,*) LOG1 If(LOG1.NE.'n') OPEN(UNIT=10, FILE='xy.dat') WRITE(*,*)'PRINT X-Z DATA TO FILE "xz.dat" ? :' READ(*,*) LOG2 If(LOG2.NE.'n') OPEN(UNIT=20, FILE='xz.dat') WRITE(*,*)'PRINT Y-Z DATA TO FILE "yz.dat" ? :' READ(*,*) LOG3 If(LOG3.NE.'n') OPEN(UNIT=30, FILE='yz.dat') WRITE(*,*)'PRINT X-Y-Z DATA TO FILE "xyz.dat" ? :' READ(*,*) LOG4 If(LOG4.NE.'n') OPEN(UNIT=40, FILE='xyz.dat') WRITE(*,*)'PRINT X-T DATA TO FILE "xt.dat" ? :' READ(*,*) LOG5 IF(LOG5.NE.'n') OPEN(UNIT=50, FILE='xt.dat') WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)' program is running ...' C---------------- MAINPART : --------------------------- DO WHILE(T.LE.T_END) STP = STP + 1 IF(LOG_A.EQ..TRUE.) THEN IF(T.GT.T_UP) THEN STP_UP = STP LOG_A = .FALSE. END IF END IF C--------------- WRITING DATAS TO FILE : --------------------- IF (T.GE.(COUNTER * T_PRINT + T_UP)) THEN IF(LOG1.NE.'n') WRITE(10,*) S(1),S(2) IF(LOG2.NE.'n') WRITE(20,*) S(1),S(3) IF(LOG3.NE.'n') WRITE(30,*) S(2),S(3) IF(LOG4.NE.'n') WRITE(40,*) S(1),S(2),S(3) IF(LOG5.NE.'n') WRITE(50,*) T, S(1) COUNTER = COUNTER + 1 WRITE(*,*) COUNTER END IF C------------------- INTEGRATION SCHEME : -------------------- T_COPY = T EST1 = S EST2 = S CALL DRKSTP(3,H,T,EST1,RHS,DSDT1) T = T_COPY CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) C---------- COMPUTE ESTIMATES AND RESULTING ERROR : --------------- EST(1) = DABS( EST1(1)-EST2(1) ) EST(2) = DABS( EST1(2)-EST2(2) ) EST(3) = DABS( EST1(3)-EST2(3) ) ERROR = EST(1) IF (EST(2).GT.ERROR) ERROR = EST(2) IF (EST(3).GT.ERROR) ERROR = EST(3) C------------ ADAPTIVE STEPSIZE : ------------------------ C-------- DETERMING NEW STEPSIZE H : --------------------- IF (ERROR.GT.TOLERANCE) THEN H_NEW = H * (TOLERANCE/ERROR)**.25d0 T = T - H END IF IF (ERROR.LE.TOLERANCE) THEN IF (ERROR.EQ.0.0d0) THEN H_NEW = 5.0d0 * H ELSE H_NEW = H * (TOLERANCE/ERROR)**.2d0 END IF S = (16.0d0 * EST2 - EST1) / 15.0d0 END IF IF (H_NEW.LT.(H/5.0d0)) then H_NEW = (H/5.0d0) END IF IF (H_NEW.GT.(5.0d0*H)) THEN H_NEW = (5.0d0*H) END IF IF (H_NEW.LT.(1.0d-15)) THEN H_NEW = 1.0d-15 END IF H = H_NEW END DO PRINT*,'...end of running !!!' C------------------ LAST SCREEN OUTPUT : ----------------------- WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)'YOUR CHOSEN CONFIGURATION WAS :' WRITE(*,*)' ' WRITE(*,*)'# TOLERANCE = ',TOLERANCE WRITE(*,*)' ' WRITE(*,*)'# DATA ARE PRINTED TO FILE BETWEEN :' WRITE(*,*)' TIME T-UP = ',T_UP,' ( STEP-NR ',STP_UP,' )' WRITE(*,*)' AND' WRITE(*,*)' TIME T_END = ',T_END,' ( STEP-NR ',STP,' )' WRITE(*,*)'# NUMBER OF WRITTEN POINTS = ',COUNTER,' OF ',NR_DATA WRITE(*,*)' ( DISTANCE OF DATA POINTS = CONST = ',T_PRINT,' )' WRITE(*,*)' ' WRITE(*,*)' # PARAMETER : SIGMA = ',SIGMA WRITE(*,*)' R = ',R WRITE(*,*)' B = ',B WRITE(*,*)' AX = ',AX WRITE(*,*)' AY = ',AY WRITE(*,*)' AZ = ',AZ WRITE(*,*)'# INITIAL CONFIGURATION : Xo = ',So(1) WRITE(*,*)' Yo = ',So(2) WRITE(*,*)' Zo = ',So(3) WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' STOP END C-------------------------------------------------------------------------- c****************************************************** c*** SUBROUTINE RHS *** c*** COMPUTES THE RIGHT HAND SIDE OF THE EQUATIONS *** c****************************************************** SUBROUTINE RHS(T, S, DSDT) IMPLICIT NONE REAL*8 S(3), DSDT(3) REAL*8 SIGMA, R, B, AX, AY, AZ, T COMMON SIGMA, R, B, AX, AY, AZ DSDT(1) = - SIGMA * S(1) + SIGMA * S(2) + AX DSDT(2) = - S(1) * S(3) + R * S(1) - S(2) + AY DSDT(3) = S(1) * S(2) - B * S(3) + AZ RETURN END C---------------------------------------------------------------------------- VARIATIONS/ EXTENSION OF THE ABOVE PROGRAM: BIF.f C------------------------------------------------------------------------- PROGRAM LORENZ1 IMPLICIT NONE REAL*8 So(3), S(3), EST1(3), EST2(3), EST(3), DSDT1(9), DSDT2(9) REAL*8 H, H_NEW, T, T_COPY, T_UP, T_END REAL*8 TOLERANCE, ERROR REAL*8 SIGMA, R, B, AX, AY, AZ ,DA REAL*8 COUNTER, T_PRINT, NR_DATA INTEGER STP, STP_UP CHARACTER LOG1, LOG2, LOG3, LOG4, LOG5 LOGICAL LOG_A COMMON SIGMA, R, B, AX, AY, AZ EXTERNAL RHS EXTERNAL DRKSTP C------------------------ INPUT : ---------------------------------- WRITE(*,*)'__________________________________________' WRITE(*,*)'ENTER TOLERANCE FOR ADAPTIVE STEPSIZE :' READ(*,*) TOLERANCE WRITE(*,*)'ENTER OVERALL INTEGRATION TIME ( T_END ) :' READ(*,*) T_END WRITE(*,*)'ENTER TIME ( T_UP ) UP WHICH DATA SHALL' WRITE(*,*)'BE PRINTED TO THE DATA FILE(S) :' READ(*,*) T_UP WRITE(*,*)'ENTER NUMBER ( NR_DATA ) OF DATA POINTS' WRITE(*,*)'PRINTED TO FILE :' READ(*,*) NR_DATA WRITE(*,*)'ENTER INITIAL Xo :' READ(*,*) S(1) WRITE(*,*)'ENTER INITIAL Yo :' READ(*,*) S(2) WRITE(*,*)'ENTER INITIAL Zo :' READ(*,*) S(3) C WRITE(*,*)'ENTER PARAMETER "SIGMA" :' C READ(*,*) SIGMA WRITE(*,*)'ENTER PARAMETER "R" :' READ(*,*) R C WRITE(*,*)'ENTER PARAMETER "B" :' C READ(*,*) B c WRITE(*,*)'ENTER AX, AY, AZ, ( I.E. FORCING) :' c READ(*,*) AX c READ(*,*) AY c READ(*,*) AZ SIGMA = 10 B = 2.666666666666666666 C-------------------- SOME VALUES : --------------------- H = 1.0d-10 T = 0.0d0 STP = 0 T_PRINT = ( (T_END - T_UP) / NR_DATA ) COUNTER = 0.0D0 LOG_A = .TRUE. So = S R = 28 Az = 0.0d0 C----------- DATA TO FILE TO BE CHOSEN : ------------------ OPEN(UNIT=10, FILE='bifx.dat') open(unit=20,file='bify.dat') 10 WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)' program is running ...' ax = 4.0 C---------------- MAINPART : --------------------------- DA = 0.2 DO WHILE(Ax.LE.30) Ay = Ax S = SO T = 0 COUNTER = 0 DO WHILE(T.LE.T_END) C--------------- WRITING DATAS TO FILE : --------------------- IF (T.GE.(COUNTER * T_PRINT + T_UP)) THEN WRITE(10,*) Ax, S(1) write(20,*) Ax, s(2) COUNTER = COUNTER + 1 WRITE(*,*) COUNTER,Ax END IF C------------------- INTEGRATION SCHEME : -------------------- T_COPY = T EST1 = S EST2 = S CALL DRKSTP(3,H,T,EST1,RHS,DSDT1) T = T_COPY CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) C---------- COMPUTE ESTIMATES AND RESULTING ERROR : --------------- EST(1) = DABS( EST1(1)-EST2(1) ) EST(2) = DABS( EST1(2)-EST2(2) ) EST(3) = DABS( EST1(3)-EST2(3) ) ERROR = EST(1) IF (EST(2).GT.ERROR) ERROR = EST(2) IF (EST(3).GT.ERROR) ERROR = EST(3) C------------ ADAPTIVE STEPSIZE : ------------------------ C-------- DETERMING NEW STEPSIZE H : --------------------- IF (ERROR.GT.TOLERANCE) THEN H_NEW = H * (TOLERANCE/ERROR)**.25d0 T = T - H END IF IF (ERROR.LE.TOLERANCE) THEN IF (ERROR.EQ.0.0d0) THEN H_NEW = 5.0d0 * H ELSE H_NEW = H * (TOLERANCE/ERROR)**.2d0 END IF S = (16.0d0 * EST2 - EST1) / 15.0d0 END IF IF (H_NEW.LT.(H/5.0d0)) then H_NEW = (H/5.0d0) END IF IF (H_NEW.GT.(5.0d0*H)) THEN H_NEW = (5.0d0*H) END IF IF (H_NEW.LT.(1.0d-15)) THEN H_NEW = 1.0d-15 END IF H = H_NEW END DO Ax = Ax + DA END DO if (log_a) then SO = - SO GOTO 10 log_a = .false. end if PRINT*,'...end of running !!!' C------------------ LAST SCREEN OUTPUT : ----------------------- WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)'YOUR CHOSEN CONFIGURATION WAS :' WRITE(*,*)' ' WRITE(*,*)'# TOLERANCE = ',TOLERANCE WRITE(*,*)' ' WRITE(*,*)'# DATA ARE PRINTED TO FILE BETWEEN :' WRITE(*,*)' TIME T-UP = ',T_UP,' ( STEP-NR ',STP_UP,' )' WRITE(*,*)' AND' WRITE(*,*)' TIME T_END = ',T_END,' ( STEP-NR ',STP,' )' WRITE(*,*)'# NUMBER OF WRITTEN POINTS = ',COUNTER,' OF ',NR_DATA WRITE(*,*)' ( DISTANCE OF DATA POINTS = CONST = ',T_PRINT,' )' WRITE(*,*)' ' WRITE(*,*)' # PARAMETER : SIGMA = ',SIGMA WRITE(*,*)' R = ',R WRITE(*,*)' B = ',B WRITE(*,*)' AX = ',AX WRITE(*,*)' AY = ',AY WRITE(*,*)' AZ = ',AZ WRITE(*,*)'# INITIAL CONFIGURATION : Xo = ',So(1) WRITE(*,*)' Yo = ',So(2) WRITE(*,*)' Zo = ',So(3) WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' STOP END C-------------------------------------------------------------------------- c****************************************************** c*** SUBROUTINE RHS *** c*** COMPUTES THE RIGHT HAND SIDE OF THE EQUATIONS *** c****************************************************** SUBROUTINE RHS(T, S, DSDT) IMPLICIT NONE REAL*8 S(3), DSDT(3) REAL*8 SIGMA, R, B, AX, AY, AZ, T COMMON SIGMA, R, B, AX, AY, AZ DSDT(1) = - SIGMA * S(1) + SIGMA * S(2) + AX DSDT(2) = - S(1) * S(3) + R * S(1) - S(2) + AY DSDT(3) = S(1) * S(2) - B * S(3) + AZ RETURN END C---------------------------------------------------------------------------- INITIAL.f C------------------------------------------------------------------------- PROGRAM LORENZ1 IMPLICIT NONE REAL*8 So(3), S(3), EST1(3), EST2(3), EST(3), DSDT1(9), DSDT2(9) REAL*8 H, H_NEW, T, T_COPY, T_UP, T_END REAL*8 TOLERANCE, ERROR REAL*8 SIGMA, R, B, AX, AY, AZ REAL*8 COUNTER, T_PRINT, NR_DATA INTEGER STP, STP_UP CHARACTER LOG1, LOG2, LOG3, LOG4, LOG5 LOGICAL LOG_A COMMON SIGMA, R, B, AX, AY, AZ EXTERNAL RHS EXTERNAL DRKSTP C------------------------ INPUT : ---------------------------------- WRITE(*,*)'__________________________________________' c WRITE(*,*)'ENTER TOLERANCE FOR ADAPTIVE STEPSIZE :' c READ(*,*) TOLERANCE WRITE(*,*)'ENTER OVERALL INTEGRATION TIME ( T_END ) :' READ(*,*) T_END WRITE(*,*)'ENTER TIME ( T_UP ) UP WHICH DATA SHALL' WRITE(*,*)'BE PRINTED TO THE DATA FILE(S) :' READ(*,*) T_UP WRITE(*,*)'ENTER NUMBER ( NR_DATA ) OF DATA POINTS' WRITE(*,*)'PRINTED TO FILE :' READ(*,*) NR_DATA WRITE(*,*)'ENTER INITIAL Xo :' READ(*,*) S(1) WRITE(*,*)'ENTER INITIAL Yo :' READ(*,*) S(2) WRITE(*,*)'ENTER INITIAL Zo :' READ(*,*) S(3) C WRITE(*,*)'ENTER PARAMETER "SIGMA" :' C READ(*,*) SIGMA c WRITE(*,*)'ENTER PARAMETER "R" :' c READ(*,*) R C WRITE(*,*)'ENTER PARAMETER "B" :' C READ(*,*) B WRITE(*,*)'ENTER AX, AY, AZ, ( I.E. FORCING) :' READ(*,*) AX READ(*,*) AY READ(*,*) AZ SIGMA = 10 B = 2.666666666666666666 r = 28.0 C-------------------- SOME VALUES : --------------------- H = 1.0d-10 T = 0.0d0 STP = 0 T_PRINT = ( (T_END - T_UP) / NR_DATA ) COUNTER = 0.0D0 LOG_A = .TRUE. So = S tolerance = 1e-14 C----------- DATA TO FILE TO BE CHOSEN : ------------------ c WRITE(*,*)'PRINT X-Y DATA TO FILE "xy.dat" ? :' c READ(*,*) LOG1 OPEN(UNIT=10, FILE='xy.dat') c WRITE(*,*)'PRINT X-Z DATA TO FILE "xz.dat" ? :' c READ(*,*) LOG2 OPEN(UNIT=20, FILE='xz.dat') c WRITE(*,*)'PRINT Y-Z DATA TO FILE "yz.dat" ? :' c READ(*,*) LOG3 OPEN(UNIT=30, FILE='yz.dat') c WRITE(*,*)'PRINT X-Y-Z DATA TO FILE "xyz.dat" ? :' c READ(*,*) LOG4 c If(LOG4.NE.'n') OPEN(UNIT=40, FILE='xyz.dat') c WRITE(*,*)'PRINT X-T DATA TO FILE "xt.dat" ? :' c READ(*,*) LOG5 c IF(LOG5.NE.'n') OPEN(UNIT=50, FILE='xt.dat') WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)' program is running ...' C---------------- MAINPART : --------------------------- DO WHILE(T.LE.T_END) STP = STP + 1 IF(LOG_A.EQ..TRUE.) THEN IF(T.GT.T_UP) THEN STP_UP = STP LOG_A = .FALSE. END IF END IF C--------------- WRITING DATAS TO FILE : --------------------- IF (T.GE.(COUNTER * T_PRINT + T_UP)) THEN IF(LOG1.NE.'n') WRITE(10,*) S(1),S(2) IF(LOG2.NE.'n') WRITE(20,*) S(1),S(3) IF(LOG3.NE.'n') WRITE(30,*) S(2),S(3) IF(LOG4.NE.'n') WRITE(40,*) S(1),S(2),S(3) IF(LOG5.NE.'n') WRITE(50,*) T, S(1) COUNTER = COUNTER + 1 c WRITE(*,*) COUNTER END IF C------------------- INTEGRATION SCHEME : -------------------- T_COPY = T EST1 = S EST2 = S CALL DRKSTP(3,H,T,EST1,RHS,DSDT1) T = T_COPY CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) C---------- COMPUTE ESTIMATES AND RESULTING ERROR : --------------- EST(1) = DABS( EST1(1)-EST2(1) ) EST(2) = DABS( EST1(2)-EST2(2) ) EST(3) = DABS( EST1(3)-EST2(3) ) ERROR = EST(1) IF (EST(2).GT.ERROR) ERROR = EST(2) IF (EST(3).GT.ERROR) ERROR = EST(3) C------------ ADAPTIVE STEPSIZE : ------------------------ C-------- DETERMING NEW STEPSIZE H : --------------------- IF (ERROR.GT.TOLERANCE) THEN H_NEW = H * (TOLERANCE/ERROR)**.25d0 T = T - H END IF IF (ERROR.LE.TOLERANCE) THEN IF (ERROR.EQ.0.0d0) THEN H_NEW = 5.0d0 * H ELSE H_NEW = H * (TOLERANCE/ERROR)**.2d0 END IF S = (16.0d0 * EST2 - EST1) / 15.0d0 END IF IF (H_NEW.LT.(H/5.0d0)) then H_NEW = (H/5.0d0) END IF IF (H_NEW.GT.(5.0d0*H)) THEN H_NEW = (5.0d0*H) END IF IF (H_NEW.LT.(1.0d-15)) THEN H_NEW = 1.0d-15 END IF H = H_NEW END DO PRINT*,'...end of running !!!' C------------------ LAST SCREEN OUTPUT : ----------------------- WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)'YOUR CHOSEN CONFIGURATION WAS :' WRITE(*,*)' ' WRITE(*,*)'# TOLERANCE = ',TOLERANCE WRITE(*,*)' ' WRITE(*,*)'# DATA ARE PRINTED TO FILE BETWEEN :' WRITE(*,*)' TIME T-UP = ',T_UP,' ( STEP-NR ',STP_UP,' )' WRITE(*,*)' AND' WRITE(*,*)' TIME T_END = ',T_END,' ( STEP-NR ',STP,' )' WRITE(*,*)'# NUMBER OF WRITTEN POINTS = ',COUNTER,' OF ',NR_DATA WRITE(*,*)' ( DISTANCE OF DATA POINTS = CONST = ',T_PRINT,' )' WRITE(*,*)' ' WRITE(*,*)' # PARAMETER : SIGMA = ',SIGMA WRITE(*,*)' R = ',R WRITE(*,*)' B = ',B WRITE(*,*)' AX = ',AX WRITE(*,*)' AY = ',AY WRITE(*,*)' AZ = ',AZ WRITE(*,*)'# INITIAL CONFIGURATION : Xo = ',So(1) WRITE(*,*)' Yo = ',So(2) WRITE(*,*)' Zo = ',So(3) WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' STOP END C-------------------------------------------------------------------------- c****************************************************** c*** SUBROUTINE RHS *** c*** COMPUTES THE RIGHT HAND SIDE OF THE EQUATIONS *** c****************************************************** SUBROUTINE RHS(T, S, DSDT) IMPLICIT NONE REAL*8 S(3), DSDT(3) REAL*8 SIGMA, R, B, AX, AY, AZ, T COMMON SIGMA, R, B, AX, AY, AZ DSDT(1) = - SIGMA * S(1) + SIGMA * S(2) + AX DSDT(2) = - S(1) * S(3) + R * S(1) - S(2) + AY DSDT(3) = S(1) * S(2) - B * S(3) + AZ RETURN END C---------------------------------------------------------------------------- C--------------------------------------------------------------------- PROGRAM LOR_DENSI IMPLICIT NONE REAL*8 So(3), S(3), EST1(3), EST2(3), EST(3), DSDT1(9), DSDT2(9) REAL*8 H, H_NEW, T, T_COPY, T_UP, T_END, H_SUM, H_COUNT REAL*8 TOLERANCE, ERROR REAL*8 SIGMA, R, B, AX, AY, AZ REAL*8 COUNTER, T_PRINT, NR_DATA, rel, PASSED INTEGER X, Y, Z INTEGER DENSITY(-20:20,-30:30,0:50) INTEGER DENSI_Z(-20:20,-30:30), DENSI_Y(-20:20,0:50) INTEGER DENSI_X(-30:30,0:50) real*8 Xp, Xm, Xz, Yp, Ym, Yz, sum_y, sum_x INTEGER PART CHARACTER LOG_Z, LOG_Y, LOG_X, LOG_XYZ LOGICAL LOG_A COMMON SIGMA, R, B, AX, AY, AZ EXTERNAL RHS EXTERNAL DRKSTP open(unit=90,file='abs_xp-xm.dat') open(unit=80,file='rel_xp-xm.dat') open(unit=99,file='t-h.dat') C---------------------- INPUT : ------------------------------------- WRITE(*,*)'__________________________________________' WRITE(*,*)'ENTER TOLERANCE FOR ADAPTIVE STEPSIZE :' READ(*,*) TOLERANCE WRITE(*,*)'ENTER OVERALL INTEGRATION TIME ( T_END ) :' READ(*,*) T_END WRITE(*,*)'ENTER TIME ( T_UP ) UP WHICH DATA SHALL' WRITE(*,*)'BE PRINTED TO THE DATA FILE(S) :' READ(*,*) T_UP WRITE(*,*)'ENTER NUMBER ( NR_DATA ) OF DATA POINTS' WRITE(*,*)'PRINTED TO FILE :' READ(*,*) NR_DATA WRITE(*,*)'ENTER INITIAL Xo :' READ(*,*) S(1) WRITE(*,*)'ENTER INITIAL Yo :' READ(*,*) S(2) WRITE(*,*)'ENTER INITIAL Zo :' READ(*,*) S(3) C WRITE(*,*)'ENTER PARAMETER "SIGMA" :' C READ(*,*) SIGMA WRITE(*,*)'ENTER PARAMETER "R" :' READ(*,*) R C WRITE(*,*)'ENTER PARAMETER "B" :' C READ(*,*) B WRITE(*,*)'ENTER "AX", "AY", "AZ", ( I.E. FORCING ) :' READ(*,*) AX READ(*,*) AY READ(*,*) AZ SIGMA = 10 B = 2.666666666666666 C----------- FILE FOR DATAS TO BE CHOSEN : ------------------- WRITE(*,*) ' ' WRITE(*,*) 'WISH TO PLOT DENSITY IN XY-PLANE ?' READ(*,*) LOG_Z WRITE(*,*) 'WISH TO PLOT DENSITY IN XZ-PLANE ?' READ(*,*) LOG_Y WRITE(*,*) 'WISH TO PLOT DENSITY IN YZ-PLANE ?' READ(*,*) LOG_X WRITE(*,*) 'WISH TO PLOT DENSITY IN XYZ-SPACE ?' READ(*,*) LOG_XYZ WRITE(*,*) 'program is running...' C------------------------ SOME VALUES : --------------- H = tolerance ! first "h" to get started T = 0.0d0 ! time starts at 0 T_PRINT = ( (T_END - T_UP) / NR_DATA ) ! interval for printing ! computed by inserted COUNTER = 0.0d0 ! values LOG_A = .TRUE. PASSED = 0.1 So = S ! intial values for last screen output Xp = 0 Xm = 0 Xz = 0 Yp = 0 Ym = 0 Yz = 0 DENSITY = 0 C------------------- MAINPART : ------------------------ DO WHILE(T.LE.T_END) IF (T.GT.1) THEN IF (MOD(NINT(T),NINT(T_END*PASSED)).EQ.0) THEN WRITE(*,*)' ', NINT(PASSED*100),' % ' PASSED = PASSED + 0.1 END IF END IF C------------- WRITING TO ARRAY : ---------------- IF (T.GE.(COUNTER * T_PRINT + T_UP)) THEN X = NINT(S(1)) Y = NINT(S(2)) Z = NINT(S(3)) DENSITY(X,Y,Z) = DENSITY(X,Y,Z) + 1 C---------- COUNTING POSITIVE AND NEGATIVE VALUES ---------- IF (S(1).GT.0) Xp = Xp + 1.0 IF (S(1).EQ.0) Xz = Xz + 1.0 IF (S(1).LT.0) Xm = Xm + 1.0 IF (S(2).GT.0) Yp = Yp + 1.0 IF (S(2).EQ.0) Yz = Yz + 1.0 IF (S(2).LT.0) Ym = Ym + 1.0 rel = (xp - xm)/(xp+xm) c write(*,*) xp,' dabs = ',dabs(xp-xm),' / = ',xp/(dabs(xp - xm)) c & ,\ =',(dabs(xp - xm))/xp write(90,*)T, (Xp - Xm) write(80,*)T, rel COUNTER = COUNTER + 1.0d0 c WRITE(*,*) COUNTER END IF C------- INTEGRATION SCHEME : ----------------------- T_COPY = T EST1 = S EST2 = S CALL DRKSTP(3,H,T,EST1,RHS,DSDT1) T = T_COPY CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) CALL DRKSTP(3,H/2.0d0,T,EST2,RHS,DSDT2) C----- COMPUTE ESTIMATES AND RESULTING ERROR : ------------ EST(1) = DABS( EST1(1)-EST2(1) ) EST(2) = DABS( EST1(2)-EST2(2) ) EST(3) = DABS( EST1(3)-EST2(3) ) ERROR = EST(1) IF (EST(2).GT.ERROR) ERROR = EST(2) IF (EST(3).GT.ERROR) ERROR = EST(3) C----------------- ADAPTIVE STEPSIZE : ------------------- C------------- DETERMING NEW STEPSIZE H : ---------------- IF (ERROR.GT.TOLERANCE) THEN H_NEW = H * (TOLERANCE/ERROR)**.25d0 T = T - H END IF IF (ERROR.LE.TOLERANCE) THEN IF (ERROR.EQ.0.0d0) THEN H_NEW = 5.0d0 * H ELSE H_NEW = H * (TOLERANCE/ERROR)**.2d0 END IF S = (16.0d0 * EST2 - EST1) / 15.0d0 END IF IF (H_NEW.LT.(H/5.0d0)) then H_NEW = (H/5.0d0) END IF IF (H_NEW.GT.(5.0d0*H)) THEN H_NEW = (5.0d0*H) END IF IF (H_NEW.LT.(1.0d-15)) THEN H_NEW = 1.0d-15 END IF H = H_NEW if (mod(h_count,1000).eq.0) then write(99,*) t, error end if H_SUM = H_SUM + H H_COUNT = H_COUNT + 1.0 END DO WRITE(*,*) '100 % ...end of running !!!' WRITE(*,*) ' ' WRITE(*,*)'AVERAGE H = ', H_SUM / H_COUNT IF (LOG_Z.NE.'N'.OR.LOG_Y.NE.'N'.OR.LOG_X.NE.'N') THEN WRITE(*,*)'printing density-data to file ...' WRITE(*,*)'Please be patient !' END IF WRITE(*,*)'______________________________________________' C---------- PRODUCE 2-DIM ARRAY(S) FROM THE 3-DIM ARRAY ---------- C---------- ( I.E. 3 POSSIBILITIES: X,Y; X,Z; Y,Z ) TO PLOT ------- C---------- A PROJECTION OF DENSITY(X,Y,Z) AS CHOSEN ABOVE : ------- IF (LOG_Z.NE.'n') THEN OPEN(UNIT=10,FILE='proba-xy.dat') DENSI_Z = 0 DO Z=0,50 DENSI_Z = DENSI_Z + DENSITY(:,:,Z) END DO DO X=-20,20 DO Y=-30,30 WRITE(10,*) X, Y, DENSI_Z(X,Y) / COUNTER END DO END DO CLOSE(10) END IF IF (LOG_Y.NE.'n') THEN OPEN(UNIT=20,FILE='proba-xz.dat') DENSI_Y = 0 DO Y=-30,30 DENSI_Y = DENSI_Y + DENSITY(:,Y,:) END DO sum_y = 0 do x=-20,20 do z=0,50 sum_y = sum_y + densi_y(x,z) end do end do DO X=-20,20 DO Z=0,50 WRITE(20,*) X, Z, DENSI_Y(X,Z)/sum_y END DO END DO CLOSE(20) END IF IF (LOG_X.NE.'n') THEN OPEN(UNIT=30,FILE='proba-yz.dat') DENSI_X = 0 DO X=-20,20 DENSI_X = DENSI_X + DENSITY(X,:,:) END DO sum_x = 0 do y=-30,30 do z =0,50 sum_x = sum_x + densi_x(y,z) end do end do DO Y=-30,30 DO Z=0,50 WRITE(30,*) Y, Z, DENSI_X(Y,Z)/sum_x END DO END DO CLOSE(30) END IF IF (LOG_XYZ.NE.'n') THEN OPEN(UNIT=40,FILE='proba-xyz.dat') PART = 0 DO X=-20,20 DO Y=-30,30 DO Z=0,50 PART = PART+1 IF (DENSITY(X,Y,Z).EQ.0) THEN IF(MOD(PART,40).EQ.0) THEN WRITE(40,*) X, Y, Z, DENSITY(X,Y,Z)/COUNTER END IF ELSE WRITE(40,*) X, Y, Z, DENSITY(X,Y,Z) / COUNTER END IF END DO END DO END DO CLOSE(40) END IF C-------------------- LAST SCREEN OUTPUT : ------------------------ WRITE(*,*)counter,' data are printed to file(s) !!!' WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)'YOUR CHOSEN CONFIGURATION WAS :' WRITE(*,*)' ' WRITE(*,*)'# TOLERANCE = ',TOLERANCE WRITE(*,*)' ' WRITE(*,*)'# DATA ARE PRINTED TO FILE(S) BETWEEN :' WRITE(*,*)' TIME T-UP = ',T_UP WRITE(*,*)' AND' WRITE(*,*)' TIME T_END = ',T_END WRITE(*,*)'# NUMBER OF WRITTEN DATA = ',COUNTER,' OF ',NR_DATA WRITE(*,*)' (TIME POINTS SHOULD BE = ',T_PRINT,' )' WRITE(*,*)' ' WRITE(*,*)' # PARAMETER : SIGMA = ',SIGMA WRITE(*,*)' R = ',R WRITE(*,*)' B = ',B WRITE(*,*)' AX = ',AX WRITE(*,*)' AY = ',AY WRITE(*,*)' AZ = ',AZ WRITE(*,*)'# INITIAL CONFIGURATION : Xo = ',So(1) WRITE(*,*)' Yo = ',So(2) WRITE(*,*)' Zo = ',So(3) WRITE(*,*)'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' WRITE(*,*)'Xp = ',Xp,' Xz = ',Xz,' Xm = ',Xm WRITE(*,*)'Xp + Xz + Xm = ',Xp + Xz + Xm WRITE(*,*)'Yp = ',Yp,' Yz = ',Yz,' Ym = ',Ym WRITE(*,*)'Y = ',Ym + Yz + Yp write(*,*)'c ',counter,' sx ',sum_x,' sy ', sum_y STOP END C--------------------------------------------------------------------------- c****************************************************** c*** SUBROUTINE RHS *** c*** COMPUTES THE RIGHT HAND SIDE OF THE EQUATIONS *** c****************************************************** SUBROUTINE RHS(T, S, DSDT) IMPLICIT NONE REAL*8 S(3), DSDT(3) REAL*8 SIGMA, R, B, AX, AY, AZ, T COMMON SIGMA, R, B, AX, AY, AZ DSDT(1) = - SIGMA * S(1) + SIGMA * S(2) + AX DSDT(2) = - S(1) * S(3) + R * S(1) - S(2) + AY DSDT(3) = S(1) * S(2) - B * S(3) + AZ RETURN END C----------------------------------------------------------------------------