C The software PS63.f in its version 1.1 of the 1st October 1997, C hereinafter referred to as « the software » has been designed and C produced by Anne Aubry and Philippe Chartier, C the researchers of the project ALADIN, a research project of the National C Computer and Automatics Institute (INRIA). C INRIA C Domaine de Voluceau C Rocquencourt C 78153 Le Chesnay Cedex. C C This Software is C copyright INRIA . 1997. C C INRIA holds all the ownership rights on the Software. The scientific C community is asked to use the SOFTWARE in order to test and evaluate it. C C INRIA freely grants the right to use the Software. Any use or C reproduction of this Software to obtain profit or for commercial ends C being subject to obtaining the prior express authorisation of INRIA. C C INRIA authorises any reproduction of this Software C C - in limits defined in clauses 9 and 10 of the Berne agreement for the C protection of literary and artistic works respectively specify in their C paragraphs 2 and 3 authorising only the reproduction and quoting of C works on the condition that : C - "this reproduction does not adversely affect the normal C exploitation of the work or cause any unjustified prejudice to the C legitimate interests of the author". C - that the quotations given by way of illustration and/or tuition C conform to the proper uses and that it mentions the source and name of C the author if this name features in the source", C C - under the condition that this file is included with any reproduction. C C Any commercial use made without obtaining the prior express agreement C of INRIA would therefore constitute a fraudulent imitation. C C the Software being currently developed, INRIA is assuming no liability, C and should not be responsible, in any manner or any case, for any direct C or indirect damages sustained by the user. C C ******************************************* C C---------------------------------------------- SUBROUTINE RKPS63(A,B,C) C---------------------------------------------- C OUTPUT PARAMETER C A : ARRAY FOR THE MATRIX MA OF THE EXPLICIT RUNGE-KUTTA METHOD C OF ORDER 3 AND PSEUDO-SYMPLECTICNESS ORDER 6 CALLED PS63 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 (A. AUBRY AND P. CHARTIER, PSEUDO-SYMPLECTIC RUNGE-KUTTA METHODS 1997) C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/METHEXP/NSEXP DIMENSION A(*),B(*),C(*) C---------------------------------------------- C /METHEXP/ : NSEXP : DIMENSION OF THE VECTORS B AND C C (NUMBER OF INTERNAL STAGES) C---------------------------------------------- NSEXP = 5 A(1) = 0.13502027922908531467870770059890109845029771948745 A(2) = -0.47268213605236986919157018872698934707969134911740 A(3) = 1.0598025041541896819941478977026920654585670044207 A(4) = -1.2165046059568853893460276983775969452327650543924 A(5) = 2.1621763021675253301167955602335235500247340216386 A(6) = -0.37234592426536003038396583048267719962403454576156 A(7) = 0.3327444303638736757817999651407277733302922435900 A(8) = -0.2088266829658723128357049049438063817099773643970 A(9) = 1.8786561773792085608958844927526700974335019616735 A(10) = -1.0025739247772099238419795529495914890538168408667 B(1) = 0.041138944570917691837239341448087001008370597280928 B(2) = 0.26732123194413937348518250968177776549784236278735 B(3) = 0.86700906289954518480142449994305187262185646586241 B(4) = -0.30547139552035758861867759375782586500945637201934 B(5) = 0.13000215610575533849483124268490922588138694608995 C(1) = 0.0D0 C(2) = 0.13502027922908531467870770059890109845029771948745 C(3) = 0.58712036810181981280257770897570271837887565530330 C(4) = 0.57332577194527991038680203137324940516793442148468 C(5) = 1.0D0 RETURN END C---------------------------------------------- SUBROUTINE INTEXP(NDPB,H,NBINT,T0,Y,HAM0,FCN,HAMCN,SOLOUT,SOLEX, + A,B,C,NFIC) C---------------------------------------------- C NUMERICAL SOLUTION OF A HAMILTONIAN PROBLEM: C Y'(T)=F(T,Y(T)) C USING AN EXPLICIT RUNGE-KUTTA METHOD C DEFINED IN A SUBROUTINE OF THE FORM C SUBROUTINE RK-NAME-OF-THE-METHOD(A,B,C) C WHERE A IS AN ARRAY FOR THE MATRIX MA OF THE EXPLICIT C RUNGE-KUTTA METHOD SUCH AS FOR I>J, C MA(I,J)=A((I-1)(I-2)/2+J) C AND B,C ARE ARRAYS FOR THE VECTORS B AND C OF THE METHOD. C THE INTEGRATION IS PERFORMED WITH CONSTANT STEP SIZE H. C F IS OF THE FORM F=-J GRADIENT(H) WHERE C J IS THE FOLLOWING MATRIX: C J = | I_D 0 | C | 0 I_D | C WHERE I_D IS THE D-DIMENSIONAL IDENTITY MATRIX C F IS DEFINED IN THE SUBROUTINE FEX C H IS DEFINED IN THE SUBROUTINE HAMEX C C INPUT PARAMETERS C ----------------- 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 : INITIAL VALUE OF Y C HAM : INITIAL VALUE OF THE HAMILTONIAN C FCN : NAME OF SUBROUTINE COMPUTING F C HAMCN : 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 A,B,C : ARRAYS FOR THE COMPONENTS OF THE EXPLICIT METHOD C NFIC : NUMBER OF THE RESULT FILE C OUTPUT PARAMETERS C ------------------ C Y : NUMERICAL SOLUTION AT TIME TEND C HAM : NUMERICAL HAMILTONIAN AT TIME TEND C---------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/METHEXP/NSEXP COMMON/TRAVAIL/YINT(1) COMMON/SORTIE/EXACT,YSOL(1) COMMON/STAT/NFCN DIMENSION A(NSEXP*(NSEXP-1)/2),B(NSEXP),C(NSEXP),Y(NDPB) EXTERNAL FCN,HAMCN,SOLOUT,SOLEX C---------------------------------------------- C /METHEXP/ : NSEXP : DIMENSION OF THE VECTORS B AND C C (NUMBER OF INTERNAL STAGES) C /TRAVAIL/ : YINT : ARRAY OF WORKING SPACE C YINT(1) IS THE FIRST COMPONENT OF C THE TRANSPOSE VECTOR : C (Y1,Y2,...,YS) C LET IKDEB BE IKDEB = NSEXP*NDPB+1 C YINT(IKDEB) IS THE FIRST COMPONENT OF C THE TRANSPOSE VECTOR : C (F(T,Y1),F(T,Y2),...,F(T,YS)) 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--------------------------------------------- IKDEB = NSEXP*NDPB+1 NFCN = 0 NOINT = 0 T = T0 ERROR = 0.0D0 HAMS = 0.0D0 CALL SOLOUT(NDPB,T0,Y,HAMS,ERROR,NFIC,NOINT) DO 100 NOINT=1,NBINT INDY = 1 DO 10 I=1,NSEXP C --- INITIALIZATION YI=Y_{NOINT-1} DO 11 J=1,NDPB 11 YINT(INDY+J-1) = Y(J) INDY=INDY+NDPB 10 CONTINUE INDY = 1 INDK = IKDEB INDA = 1 DO 12 I=2,NSEXP C --- COMPUTATION OF K_{I-1}=F(T,Y_{I-1}) TI = T+C(I-1)*H CALL FCN(NDPB,TI,YINT(INDY),YINT(INDK)) NFCN = NFCN+1 C --- COMPUTATION OF Y_{NOINT} DO 13 J=1,NDPB Y(J) = Y(J) + H*B(I-1)*YINT(INDK+J-1) 13 CONTINUE C --- COMPUTATION OF Y_{NI} INDYI=INDY+NDPB DO 16 NI=I,NSEXP INDA=(NI-2)*(NI-1)/2+I-1 DO 15 J=1,NDPB YINT(INDYI+J-1) = YINT(INDYI+J-1) + H*A(INDA)*YINT(INDK+J-1) 15 CONTINUE INDYI=INDYI+NDPB 16 CONTINUE INDY = INDY+NDPB INDK = INDK + NDPB 12 CONTINUE C --- COMPUTATION OF K_S=F(T,YS) TI = T+C(NSEXP)*H CALL FCN(NDPB,TI,YINT(INDY),YINT(INDK)) NFCN = NFCN + 1 C --- COMPUTATION DE Y_{NOINT} DO 14 J=1,NDPB 14 Y(J) = Y(J) + H*B(NSEXP)*YINT(INDK+J-1) C --- COMPUTATION OF THE HAMILTONIAN T = T+H CALL HAMCN(NDPB,T,Y,HAM) C --- TREATMENTS OF THE RESULTS IF (EXACT.LE.1) THEN CALL SOLEX(NDPB,T) ERROR = 0.0D0 DO 20 J=1,NDPB AUX = YSOL(J)-Y(J) ERROR = ERROR+ AUX*AUX 20 CONTINUE ERROR=DSQRT(ERROR) ENDIF HAMS = HAM-HAM0 CALL SOLOUT(NDPB,T,Y,HAMS,ERROR,NFIC,NOINT) 100 CONTINUE RETURN END C----------------------------------------------