CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DRIVER FOR PS63 AT GALACTIC PROBLEM CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C "GALACTIC PROBLEM" C E. HAIRER, S.P. NORSET, G. WANNER, C SOLVING ORDINARY DIFFERENTIAL EQUATIONS, VOL. 1, C SPRINGER-VERLAG, 1991, SECOND EDITION, SECTION II.16 C---------------------------------------------- C LINK dr_Galactic.f RKPS63.f C------------------------------------------ PROGRAM GALACTIC IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(ND=6,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/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2 COMMON/METHEXP/NSEXP COMMON/TRAVAIL/YINT COMMON/SORTIE/EXACT,YSOL COMMON/STAT/NFCN C---------------------------------------------- C /GALAXY/ : CONSTANTS FOR THE "GALACTIC PROBLEM" 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 --- INITIALIZATION OF CONSTANTS SPECIFIC TO THE GALACTIC PROBLEM PAG = 1.25D0 PBG = 1.0D0 PCG = 0.75D0 PAG2 = PAG*PAG PBG2 = PBG*PBG PCG2 = PCG*PCG AG = 1.0D0 CG = 1.0D0 OMEGAG = 0.25D0 C --- DIMENSION OF THE SYSTEM NDPB = 6 C --- INITIAL VALUES T0 = 0.0D0 Y(1) = 0.0D0 Y(2) = 1.688837005904475591 Y(3) = 0.2D0 Y(4) = 2.5D0 Y(5) = 0.0D0 Y(6) = 0.0D0 C --- INITIAL HAMILTONIAN ("EXACT VALUE" OF THE HAMILTONIAN) HAM0 = 2.0D0 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) COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2 P1 = Y(1) P2 = Y(2) P3 = Y(3) Q1 = Y(4) Q2 = Y(5) Q3 = Y(6) QNR = CG + Q1*Q1/PAG2 + Q2*Q2/PBG2 + Q3*Q3/PCG2 RIQNR = 1.0D0/QNR PNR = P1*P1 + P2*P2 + P3*P3 FY(1) = OMEGAG*P2 - 2.0D0*AG*Q1*RIQNR/PAG2 FY(2) = -OMEGAG*P1 - 2.0D0*AG*Q2*RIQNR/PBG2 FY(3) = -2.0D0*AG*Q3*RIQNR/PCG2 FY(4) = P1 + OMEGAG*Q2 FY(5) = P2 - OMEGAG*Q1 FY(6) = P3 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) COMMON/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2 P1 = Y(1) P2 = Y(2) P3 = Y(3) Q1 = Y(4) Q2 = Y(5) Q3 = Y(6) QNR = CG + Q1*Q1/PAG2 + Q2*Q2/PBG2 + Q3*Q3/PCG2 RIQNR = 1.0D0/QNR PNR = P1*P1 + P2*P2 + P3*P3 HAM = 0.50D0*PNR + OMEGAG*(P1*Q2-P2*Q1) + AG*DLOG(QNR) 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/GALAXY/PAG,PBG,PCG,AG,CG,OMEGAG,PAG2,PBG2,PCG2 COMMON/SORTIE/EXACT,YSOL(6) 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