CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DRIVER FOR PS63 AT MODIFIED PENDULUM PROBLEM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C "MODIFIED PENDULUM PROBLEM" C E. HAIRER, S.P. NORSET, G. WANNER, C SOLVING ORDINARY DIFFERENTIAL EQUATIONS, C VOL. 1, SPRINGER-VERLAG, SECOND EDITION, 1991 C---------------------------------------------- C LINK dr_Pendulum.f RKPS63.f C------------------------------------------ PROGRAM PENDULUM IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(ND=2,NS=5,NDNS=ND*NS) C---------------------------------------------- C ND : MAXIMAL DIMENSION OF THE PROBLEM C NS : MAXIMAL NUMBER OF INTERNAL STAGES OF C THE EXPLICIT RUNGE-KUTTA METHOD C---------------------------------------------- PARAMETER(LYINT=2*NDNS,LYSOL=ND,LY=ND,LA=NS*(NS-1)/2,LZINT=2*NDNS) DIMENSION Y(LY),YSOL(LYSOL),YINT(LYINT),A(LA),B(NS),C(NS) CHARACTER*20 NOMFIC C---------------------------------------------- C Y : ARRAY FOR THE NUMERICAL SOLUTION AT A GIVEN TIME C YSOL : ARRAY FOR THE EXACT SOLUTION AT A GIVEN TIME C YINT : ARRAY FOR WORKING SPACE C A : ARRAY FOR THE MATRIX MA OF THE EXPLICIT RUNGE-KUTTA METHODS C FOR I>J, MA(I,J)=A((I-1)(I-2)/2+J) C B,C : ARRAYS FOR THE VECTORS B AND C OF THE EXPLICIT RUNGE-KUTTA METHOD C NOMFIC : NAME OF THE RESULT FILE C---------------------------------------------- EXTERNAL FEX,HAMEX,SOLOUT,SOLEX EXTERNAL INTEXP,RKPS63 C---------------------------------------------- C FEX : NAME OF SUBROUTINE COMPUTING THE VALUE OF F(T,Y) C OF THE FORM C SUBROUTINE FEX(NDPB,T,Y,FY) WHERE C NDPB : DIMENSION OF Y AND F C T : A GIVEN TIME C Y : A GIVEN VECTOR C FY : VALUE OF F(T,Y) C HAMEX : NAME OF SUBROUTINE COMPUTING THE VALUE OF THE HAMILTONIAN H(T,Y) C OF THE FORM C SUBROUTINE HAMEX(NDPB,T,Y,HAM) WHERE C NDPB : DIMENSION OF Y C T : A GIVEN TIME C Y : A GIVEN VECTOR C HAM : VALUE OF H(T,Y) C SOLOUT : NAME OF SUBROUTINE PROVIDING THE NUMERICAL SOLUTION DURING C INTEGRATION OF THE FORM C SUBROUTINE SOLOUT(NDPB,T,Y,HAM,ERR,NFIC,NOINT) WHERE C NDPB : DIMENSION OF Y C T : TIME OF THE OUTPUT C Y : NUMERICAL SOLUTION AT TIME T C HAM : VALUE OF THE HAMILTONIAN H(T,Y) C ERR : EUCLIDEAN NORM OF THE VECTOR ERROR : Y_NUMERIC-Y_EXACT C (IF Y_EXACT IS KNOWN) C NFIC : NUMBER OF THE RESULT FILE C NOINT : NUMBER OF PERFORMED INTEGRATIONS AT TIME T C SOLEX : NAME OF SUBROUTINE COMPUTING THE EXACT SOLUTION C (IF THERE EXISTS AN EXPLICIT FORM) OF THE FORM C SUBROUTINE SOLEX(NDPB,T) WHERE C NDPB : DIMENSION OF Y C T : A GIVEN TIME C---------------------------------------------- C INTEXP : NAME OF INTEGRATOR SUBROUTINE C OF THE FORM C SUBROUTINE INTEXP(NDPB,H,NBINT,T0,Y,HAM,FEX,HAMEX, C SOLOUT,SOLEX,A,B,C,NFIC) WHERE C NDPB : DIMENSION OF Y AND F C H : STEP SIZE C NBINT : NUMBER OF INTEGRATIONS C T0 : INITIAL TIME C REMARK : TEND=T0+NBINT*H C Y : INPUT : INITIAL VALUE OF Y C OUTPUT : NUMERICAL SOLUTION AT TIME TEND C HAM : INPUT : INITIAL VALUE OF THE HAMILTONIAN C OUTPUT : NUMERICAL HAMILTONIAN AT TIME TEND C FEX : NAME OF SUBROUTINE COMPUTING F C HAMEX : NAME OF SUBROUTINE COMPUTING H C SOLOUT : NAME OF SUBROUTINE PROVIDING THE NUMERICAL C SOLUTION DURING INTEGRATION C SOLEX : NAME OF SUBROUTINE COMPUTING EXACT SOLUTION C NFIC : NUMBER OF THE RESULT FILE C RKPS63 : NAME OF SUBROUTINE DEFINING THE EXPLICIT C RUNGE-KUTTA METHOD OF ORDER 3 AND C PSEUDO-SYMPLECTICNESS ORDER 6 WITH 5 STAGES, C CALLED PS63. C IT IS OF THE FORM C SUBROUTINE RK-NAME-OF-THE-METHOD(A,B,C) WHERE C A : ARRAY FOR THE MATRIX MA OF THE METHOD C FOR I>J, MA(I,J)=A((I-1)(I-2)/2+J) C B,C : ARRAYS FOR THE VECTORS B AND C OF THE METHOD C---------------------------------------------- COMMON/METHEXP/NSEXP COMMON/TRAVAIL/YINT COMMON/SORTIE/EXACT,YSOL COMMON/STAT/NFCN C---------------------------------------------- C /METHEXP/ : NSEXP : DIMENSION OF THE VECTORS B AND C C OF THE METHOD PS63 C (NUMBER OF INTERNAL STAGES) C /TRAVAIL/ : YINT : ARRAY OF WORKING SPACE C /SORTIE/ : EXACT : IF THERE EXISTS AN EXACT SOLUTION THEN C EXACT<=1 ELSE EXACT>1 C YSOL : VECTOR Y=(P(T),Q(T)) OF THE EXACT SOLUTION C AT A GIVEN TIME T C /STAT/ : NFCN : NUMBER OF F EVALUATIONS C---------------------------------------------- C --- DEFINITION OF THE RESULT FILE WRITE(*,*)'NAME OF THE RESULT FILE?' READ(*,*) NOMFIC WRITE(*,*)'THE RESULT FILE IS : ',NOMFIC NFIC = 10 OPEN(NFIC,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFIC) WRITE(NFIC,101) 101 FORMAT('#',9X,'T',16X,' P',16X,' Q',16X,' HAM',16X,'ERR',10X,'NFCN') C --- DIMENSION OF THE SYSTEM NDPB = 2 C --- INITIAL VALUES T0 = 0.0D0 Y(1) = 0.0D0 Y(2) = DACOS(-0.8D0) C --- INITIAL HAMILTONIAN ("EXACT VALUE" OF THE HAMILTONIAN) HAM0 = 0.8D0 C --- THERE EXISTS NO EXPLICIT FORMULATION OF THE EXACT SOLUTION C --- HENCE EXACT>1 EXACT = 2.0D0 C --- INITIAL STEP SIZE WRITE(*,*)'STEP SIZE H? ' READ(*,*) H C --- NUMBER OF INTEGRATION NBINT OF THE FORM WRITE(*,*)'NUMBER OF INTEGRATION? ' READ(*,*) NBINT C --- CALL OF THE SUBROUTINE RKPS63 CALL RKPS63(A,B,C) C --- CALL OF THE SUBROUTINE INTEXP CALL INTEXP(NDPB,H,NBINT,T0,Y,HAM0,FEX,HAMEX,SOLOUT,SOLEX, + A,B,C,NFIC) C --- CLOSE THE RESULT FILE CLOSE(NFIC) STOP END C---------------------------------------------- SUBROUTINE FEX(NDPB,T,Y,FY) C---------------------------------------------- C INPUT PARAMETERS C NDPB : DIMENSION OF THE PROBLEM (OF Y AND FY) C T : A GIVEN TIME C Y : A GIVEN VECTOR C OUTPUT PARAMETERS C FY : VALUE OF F(T,Y) C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(NDPB),FY(NDPB) U6 = 0.166666666666667 Q = Y(2) P = Y(1) FY(1) = -(1.0D0-U6*P)*DSIN(Q) FY(2) = P + U6*DCOS(Q) RETURN END C---------------------------------------------- SUBROUTINE HAMEX(NDPB,T,Y,HAM) C---------------------------------------------- C INPUT PARAMETERS C NDPB : DIMENSION OF THE PROBLEM (OF Y AND FY) C T : A GIVEN TIME C Y : A GIVEN VECTOR C OUTPUT PARAMETERS C HAM : VALUE OF THE HAMILTONIAN H(T,Y) C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(NDPB) U6 = 0.166666666666667 Q = Y(2) P = Y(1) HAM = 0.50D0*P*P -(1.0D0-P*U6)*DCOS(Q) RETURN END C---------------------------------------------- SUBROUTINE SOLEX(NDPB,T) C---------------------------------------------- C INPUT PARAMETERS C NDPB : DIMENSION OF THE PROBLEM C T : A GIVEN TIME C OUTPUT PARAMETERS C YSOL : EXACT SOLUTION AT TIME T (IF THERE EXISTS AN C EXPLICIT FORM) C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/SORTIE/EXACT,YSOL(2) WRITE(*,*)'NO EXPLICIT FORMULATION FOR THE EXACT SOLUTION' RETURN END C---------------------------------------------- SUBROUTINE SOLOUT(NDPB,T,Y,HAM,ERROR,NFIC,NOINT) C---------------------------------------------- C INPUT PARAMETERS C NDPB : DIMENSION OF Y C T : TIME OF THE OUTPUT C Y : NUMERICAL SOLUTION AT TIME T C HAM : VALUE OF THE HAMILTONIAN H(T,Y) C ERROR : EUCLIDEAN NORM OF THE VECTOR ERROR : Y_NUMERIC-Y_EXACT C (IF Y_EXACT IS KNOWN) C NFIC : NUMBER OF THE RESULT FILE C NOINT : NUMBER OF PERFORMED INTEGRATIONS AT TIME T C PURPOSE : WRITE ON THE RESULT FILE NUMBER "NFIC" THE TIME T, TWO C CHOSEN COMPONENTS OF THE NUMERICAL SOLUTION AT TIME T, C THE NUMERICAL HAMILTONIAN, THE ERROR AND THE NUMBER OF C F EVALUATIONS C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/STAT/NFCN DIMENSION Y(NDPB) C --- FOR AN OUTPUT AT ALL THE 100 INTEGRATIONS : NR = MOD(NOINT,100) IF((NOINT.EQ.0).OR.(NR.EQ.0))THEN WRITE(NFIC,110)T,Y(1),Y(2),HAM,ERROR,NFCN ENDIF 110 FORMAT(1X,5E18.10,I10) RETURN END