(original) (raw)

gnicodes/ 40755 12122 12120 0 7624126575 11324 5ustar hairermathgnicodes/Readme100644 12122 12120 4411 7561743075 12541 0ustar hairermath This folder contains some programs that are used for computations and illustrations of the monograph GEOMETRIC NUMERICAL INTEGRATION Structure-Preserving Algorithms for Ordinary Differential Equations E. Hairer, Ch. Lubich and G. Wanner Springer Series in Computational Mathematics 31, Springer-Verlag, 2002 ISBN 3-540-43003-2 A MATLAB version of these codes is also available. Its description is given in the publication: E. Hairer and M. Hairer, GNICODE - Matlab Programs for Geometric Numerical Integration in: Frontiers in Numerical Analysis (Durham 2002), Springer, Berlin, 2003. REMARK: Some programs use the GGG graphics library. The Fortran program GGG.f is available from "http://www.unige.ch/math/folks/hairer" CONTENTS: gni_irk2.f implicit Runge-Kutta code (Gauss methods) for second order differential equations gni_comp.f symmetric composition methods based on second order symmetric methods gni_lmm2.f symmetric linear multistep methods for second order differential equations dr_irk2.f driver for gni_irk2.f dr_comp.f driver for gni_comp.f dr_lmm2.f driver for gni_lmm2.f problem.f data and equations for the problems: Kepler problem, harmonic oscillator, pendulum, and outer solar system FOLDER comparison for Kepler problem with 200 periods, outer solar system all three methods are of order 8 cal_irk2.f computation for Gauss methot cal_comp.f computation for composition method cal_lmm2.f computation for linear multistep method compar_dess.f makes the figure (uses GGG library) FOLDER examples rigidbody.f gives an application of composition methods two basic methods are presented: "quater" splitting and "ratori" rattle twobdsph.f two body problem on the sphere rattle as basic method for composition use of "last" is illustrated henon.f computation of the Poincar� section for Henon-Heiles problem gnicodes/examples/ 40755 12122 12120 0 7571432012 13125 5ustar hairermathgnicodes/examples/rigidbody.f100644 12122 12120 15570 7571430053 15400 0ustar hairermathccfeh rigidbody include '../gni_comp.f' IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Q(4),AM(3),QI(4),PI(4),RPAR(20),IPAR(20) REAL TIME0,TIME1 EXTERNAL QUATER,RATORI,POTENP,SOLFIX C--- AI1=0.5D0 AI2=0.9D0 AI3=1.0D0 RPAR(11)=AI1 RPAR(12)=AI2 RPAR(13)=AI3 IOUT=1 C--- XEND=1.0D3 H=0.6D0 DO IH=1,5 H=H/2.0d0 NSTEP=XEND/H N=4 C --- INITIAL VALUES X=0.d0 Q(1)=0.4D0 Q(2)=0.2D0 Q(3)=0.4D0 Q(4)=SQRT(1.0D0-Q(1)**2-Q(2)**2-Q(3)**2) AM(1)=0.2D0 AM(2)=1.0D0 AM(3)=0.4D0 C --- METHC=45 C DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO CALL CPU_TIME(TIME0) C --- COMPOSITION WITH SPLITTING AS BSIC C CALL GNI_COMP(QUATER,N,POTENP,NSTEP,X,AM,Q,XEND, C & METHC,SOLFIX,IOUT,RPAR,IPAR) C --- COMPOSITION WITH RATTLE AS BSIC CALL GNI_COMP(RATORI,N,POTENP,NSTEP,X,AM,Q,XEND, & METHC,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- WRITE (6,*) ' HAMILT = ',RPAR(14) WRITE (6,*) ' H = ',H,' ERR HAM = ',RPAR(15)/RPAR(14), & ' TIME = ',TIME1-TIME0 C --- PRINT FINAL SOLUTION END DO STOP END C SUBROUTINE HAMIL (Q,AM,HAM,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(4),AM(3) DIMENSION IPAR(*),RPAR(*) HAM=AM(1)**2/RPAR(11)+AM(2)**2/RPAR(12)+AM(3)**2/RPAR(13) CALL POTEN(Q,POT,RPAR,IPAR) HAM=HAM/2.0D0+POT RETURN END c SUBROUTINE POTEN(Q,POT,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(4) DIMENSION IPAR(*),RPAR(*) POT=Q(1)**2-Q(2)**2-Q(3)**2+Q(4)**2 C POT=0.0D0 RETURN END c SUBROUTINE POTENP(Q,POTP,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(4),POTP(3) DIMENSION IPAR(*),RPAR(*) POTP(1)=-2*(Q(1)*Q(2)+Q(3)*Q(4)) POTP(2)=-2*(Q(1)*Q(3)-Q(2)*Q(4)) C POTP(1)=0.0d0 C POTP(2)=0.0d0 POTP(3)=0.0d0 RETURN END C SUBROUTINE SOLFIX (NR,XOLD,X,AM,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),AM(N) DIMENSION IPAR(20),RPAR(20) CALL HAMIL (Q,AM,HAM,RPAR,IPAR) IF(NR.EQ.0) THEN RPAR(14)=HAM RPAR(15)=0.0D0 ELSE RPAR(15)=MAX(ABS(HAM-RPAR(14)),RPAR(15)) IF (ABS(HAM-RPAR(14)).GE.1.D0) THEN WRITE (6,*) ' EXIT IN SOLFIX ' IRTRN=-1 END IF END IF RETURN END C SUBROUTINE QUATER (N,X,AM,Q,EAM,EQ,POTENP,HA,HB,HC, & FIRST,LAST,RPAR,IPAR) C--- SPLITTING METHOD FOR RIGID BODY SIMULATIONS C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),EAM(3),EQ(4),POTP(3) DIMENSION IPAR(20),RPAR(20) COMMON /INTERN/FAC01,FAC02,FAC03 LOGICAL FIRST,LAST AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) IF (FIRST) THEN CALL POTENP(Q,POTP,RPAR,IPAR) C --- HALF-STEP POTENTIAL AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) FAC01=(AI2-AI3)/(AI2*AI3) FAC02=0.5D0*(AI2-AI1)/(AI2*AI1) FAC03=1.0D0/AI2 ELSE AM1=AM(1) AM2=AM(2) AM3=AM(3) END IF HD2=HB/2.0D0 FAC1=HD2*FAC01 FAC2=HD2*FAC02 FAC3=HD2*FAC03 C --- HALF-STEP "R" ALPH=FAC2*AM1 SA=SIN(ALPH) CA=COS(ALPH) CAP=1-2*SA**2 SAP=2*CA*SA SAVE=CAP*AM2+SAP*AM3 AM3=-SAP*AM2+CAP*AM3 AM2=SAVE C --- Q0=CA*Q(1)-SA*Q(2) Q1=CA*Q(2)+SA*Q(1) Q2=CA*Q(3)+SA*Q(4) Q3=CA*Q(4)-SA*Q(3) C --- STEP "S" BETA=FAC1*AM3 ANORM=SQRT(AM1**2+AM2**2+AM3**2) GAM=ANORM*FAC3 CA=COS(GAM) SA=SIN(GAM) CB=COS(BETA) SB=SIN(BETA) CBP=1-2*SB**2 SBP=2*CB*SB FAC=SA/ANORM V1=AM1*FAC V2=AM2*FAC V3=AM3*FAC C --- QP0=CA*Q0-V1*Q1-V2*Q2-V3*Q3 QP1=CA*Q1+V1*Q0+V3*Q2-V2*Q3 QP2=CA*Q2+V2*Q0+V1*Q3-V3*Q1 QP3=CA*Q3+V3*Q0+V2*Q1-V1*Q2 C --- Q0=CB*QP0-SB*QP3 Q1=CB*QP1+SB*QP2 Q2=CB*QP2-SB*QP1 Q3=CB*QP3+SB*QP0 C --- SAVE=CBP*AM1+SBP*AM2 AM2=-SBP*AM1+CBP*AM2 AM1=SAVE C --- HALF-STEP "R" ALPH=FAC2*AM1 SA=SIN(ALPH) CA=COS(ALPH) CAP=1-2*SA**2 SAP=2*CA*SA SAVE=CAP*AM2+SAP*AM3 AM3=-SAP*AM2+CAP*AM3 AM2=SAVE C --- Q(2)=CA*Q1+SA*Q0 Q(3)=CA*Q2+SA*Q3 Q(4)=CA*Q3-SA*Q2 Q1=CA*Q0-SA*Q1 C --- PROJECTION SUM=Q(2)**2+Q(3)**2+Q(4)**2 Q(1)=(Q1**2-SUM+1.0D0)/(2*Q1) C --- HALF-STEP POTENTIAL CALL POTENP(Q,POTP,RPAR,IPAR) AM(1)=AM1-HC*POTP(1) AM(2)=AM2-HC*POTP(2) AM(3)=AM3-HC*POTP(3) RETURN END C SUBROUTINE RATORI (N,X,AM,Q,EAM,EQ,POTENP,HA,HB,HC, & FIRST,LAST,RPAR,IPAR) C--- RATTLE APPLIED TO FORMULATION WITH P AND Q (MATRICES) C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),EAM(3),EQ(4),POTP(3) DIMENSION IPAR(20),RPAR(20) LOGICAL FIRST,LAST AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) IF (FIRST) THEN CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) ELSE AM1=AM(1) AM2=AM(2) AM3=AM(3) END IF EPS=1.D-16/HB HD2=HB/2.0D0 VDH=4.0D0/HB FAC1=(AI2-AI3)/AI1 FAC2=(AI3-AI1)/AI2 FAC3=(AI1-AI2)/AI3 FAD1=(AI2-AI3)*VDH FAD2=(AI3-AI1)*VDH FAD3=(AI1-AI2)*VDH C --- SOLVE FOR INTERNAL STAGE BM1=AM1*HD2/AI1 BM2=AM2*HD2/AI2 BM3=AM3*HD2/AI3 SUM=BM1**2+BM2**2+BM3**2 CM0=1.0D0-SUM/2.0D0 CM1=(BM1+FAC1*BM2*BM3)/CM0 CM2=(BM2+FAC2*CM1*BM3)/CM0 CM3=(BM3+FAC3*CM1*CM2)/CM0 DO I=1,10 SUM=CM1**2+CM2**2+CM3**2 STE=SUM+CM0**2-1.0D0 CM0=(CM0**2-SUM+1.0D0)/(2*CM0) CM1=(BM1+FAC1*CM2*CM3)/CM0 CM2=(BM2+FAC2*CM1*CM3)/CM0 CM3=(BM3+FAC3*CM1*CM2)/CM0 IF (ABS(STE).LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q1=Q(1) Q2=Q(2) Q3=Q(3) Q4=Q(4) Q(2)=CM0*Q2+CM1*Q1+CM3*Q3-CM2*Q4 Q(3)=CM0*Q3+CM2*Q1+CM1*Q4-CM3*Q2 Q(4)=CM0*Q4+CM3*Q1+CM2*Q2-CM1*Q3 Q1=CM0*Q1-CM1*Q2-CM2*Q3-CM3*Q4 C --- PROJECTION SUM=Q(2)**2+Q(3)**2+Q(4)**2 Q(1)=(Q1**2-SUM+1.0D0)/(2*Q1) C --- UPDATE M CALL POTENP(Q,POTP,RPAR,IPAR) AM(1)=AM1-HC*POTP(1)+FAD1*CM2*CM3 AM(2)=AM2-HC*POTP(2)+FAD2*CM1*CM3 AM(3)=AM3-HC*POTP(3)+FAD3*CM1*CM2 RETURN END gnicodes/examples/twobdsph.f100644 12122 12120 10200 7535336354 15250 0ustar hairermathccfeh twobdsph include '../gni_comp.f' IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM) DIMENSION RPAR(20),IPAR(20) EXTERNAL GRPOT,SOLFIX,RATTWO METHC=817 WRITE (6,*) ' METHOD ',METHC H=0.002D0 XEND=2.0D1 NSTEP=XEND/H C --- INITIAL VALUES CALL INDATA(X,XEND,NDIM,NN,Q,P,QEX,PEX,RPAR,IPAR) IOUT=1 CALL GNI_COMP(RATTWO,NN,GRPOT,NSTEP,X,P,Q,XEND, & METHC,SOLFIX,IOUT,RPAR,IPAR) WRITE (6,*) ' ERROR IN HAMIL ',RPAR(12) STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N),RPAR(20),IPAR(20) COMMON /INTERN/HAMIL0 CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(11)=HAMIL WRITE (6,*) ' HAMILTONIAN ',HAMIL RPAR(12)=0.0D0 ELSE RPAR(12)=MAX(RPAR(12),ABS(HAMIL-RPAR(11))) PROD=Q(1)*Q(4)+Q(2)*Q(5)+Q(3)*Q(6) C WRITE (6,*) NR,HAMAX,PROD END IF RETURN END C C SUBROUTINE INDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(NDIM),P(NDIM),QEX(NDIM),PEX(NDIM) N=6 X=0.D0 PHI1=1.3D0 THE1=2.1D0 Q(1)=COS(PHI1)*SIN(THE1) Q(2)=SIN(PHI1)*SIN(THE1) Q(3)=COS(THE1) PHI2=-2.1D0 THE2=-1.1D0 Q(4)=COS(PHI2)*SIN(THE2) Q(5)=SIN(PHI2)*SIN(THE2) Q(6)=COS(THE2) DPHI1=1.2D0 DTHE1=0.1D0 P(1)=-DPHI1*SIN(PHI1)*SIN(THE1)+DTHE1*COS(PHI1)*COS(THE1) P(2)=DPHI1*COS(PHI1)*SIN(THE1)+DTHE1*SIN(PHI1)*COS(THE1) P(3)=-DTHE1*SIN(THE1) DPHI2=0.1D0 DTHE2=-0.5D0 P(4)=-DPHI2*SIN(PHI2)*SIN(THE2)+DTHE2*COS(PHI2)*COS(THE2) P(5)=DPHI2*COS(PHI2)*SIN(THE2)+DTHE2*SIN(PHI2)*COS(THE2) P(6)=-DTHE2*SIN(THE2) RETURN END C SUBROUTINE HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Q(N),P(N) HAMIL=0.0D0 DO I=1,N HAMIL=HAMIL+P(I)**2 END DO CALL POTEN (N,X,Q,POT,RPAR,IPAR) HAMIL=HAMIL*0.5D0+POT RETURN END C SUBROUTINE POTEN (N,X,Q,POT,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION Q(N) PROD=0.0D0 DO I=1,3 PROD=PROD+Q(I)*Q(I+3) END DO POT=-1.0D0/ACOS(PROD) C POT=Q(3)+Q(6) RETURN END C SUBROUTINE GRPOT(N,X,Q,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),F(N) PROD=0.0D0 DO I=1,3 PROD=PROD+Q(I)*Q(I+3) END DO FAC=ACOS(PROD)**2*SQRT(1.0D0-PROD**2) DO I=1,3 II=I+3 F(I)=-Q(II)/FAC F(II)=-Q(I)/FAC END DO RETURN END C SUBROUTINE RATTWO (N,X,P,Q,EP,EQ,FCN,HA,HB,HC, & FIRST,LAST,RPAR,IPAR) C--- RATTLE FOR THE TWO-BODY PROBLEM ON THE SPHERE C--- FCN(N,X,Q,F,...) COMPUTES THE GRADIENT OF THE POTENTIAL IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N),Q(N),EP(N),EQ(N),F(100) LOGICAL FIRST,LAST C CALL FCN(N,X,Q,F,RPAR,IPAR) DO I=1,N EP(I)=P(I)-HA*F(I) EQ(I)=Q(I)+HB*EP(I) END DO EE1=0.0D0 EQ1=0.0D0 EE2=0.0D0 EQ2=0.0D0 DO I=1,3 EE1=EE1+EQ(I)**2 EQ1=EQ1+EQ(I)*Q(I) II=I+3 EE2=EE2+EQ(II)**2 EQ2=EQ2+EQ(II)*Q(II) END DO BET1=1.0D0-EE1 ALAM1=-BET1/(HB*(EQ1+SQRT(BET1+EQ1**2))) BET2=1.0D0-EE2 ALAM2=-BET2/(HB*(EQ2+SQRT(BET2+EQ2**2))) DO I=1,3 II=I+3 P(I)=EP(I)-ALAM1*Q(I) P(II)=EP(II)-ALAM2*Q(II) Q(I)=Q(I)+HB*P(I) Q(II)=Q(II)+HB*P(II) END DO C IF (LAST) THEN CALL FCN(N,X,Q,F,RPAR,IPAR) AMU1=0.0D0 AMU2=0.0D0 DO I=1,3 II=I+3 P(I)=P(I)-HC*F(I) P(II)=P(II)-HC*F(II) AMU1=AMU1+P(I)*Q(I) AMU2=AMU2+P(II)*Q(II) END DO DO I=1,3 II=I+3 P(I)=P(I)-AMU1*Q(I) P(II)=P(II)-AMU2*Q(II) END DO END IF RETURN END gnicodes/examples/henon.f100644 12122 12120 7036 7535351065 14516 0ustar hairermathC * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR DOPRI5 ON VAN DER POL'S EQUATION C * * * * * * * * * * * * * * * * * * * * * * * * * ccfeh henon include '../gni_irk2.f' include '../gni_lmm2.f' include '../gni_comp.f' IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDGL=2) DIMENSION P(NDGL),Q(NDGL),IPAR(20),RPAR(20) EXTERNAL HENON,SOLFIX,STVERL DO IGNI=1,3 C --- DIMENSION OF THE SYSTEM N=2 C --- INITIAL VALUES X=0.0D0 VAL=0.18D0 P(1)=VAL P(2)=VAL Q(1)=VAL Q(2)=VAL U=0.5D0*(Q(1)**2+Q(2)**2)+Q(1)**2*Q(2)-Q(2)**3/3.0D0 T=0.5D0*(P(1)**2+P(2)**2) WRITE (6,*) ' HAMIL ',T+U C --- ENDPOINT OF INTEGRATION X=0.0D0 XEND=20.0D0 IOUT=1 IF (IGNI.EQ.1) H=0.22D0 IF (IGNI.EQ.2) H=1.5D0 IF (IGNI.EQ.3) H=1.2D0 NSTEP=(XEND-X)/H DO I=1,20 IPAR(I)=0 RPAR(I)=0.0D0 END DO WRITE (6,*) ' NSTEP ',NSTEP IF (IGNI.EQ.1) THEN METH=803 CALL GNI_LMM2(N,HENON,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) END IF IF (IGNI.EQ.2) THEN METH=4 CALL GNI_IRK2(N,HENON,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) END IF IF (IGNI.EQ.3) THEN METH=817 CALL GNI_COMP(STVERL,N,HENON,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) END IF C --- PRINT FINAL SOLUTION WRITE (6,99) X,Q(1),Q(2) 99 FORMAT(1X,'X =',F10.2,' Q =',2E18.10) END DO STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N),IPAR(20),RPAR(20) REAL XD,YD COMMON /INTERN/XOUT,EVOLD,HAM0,HAMMAX,QOLD(10),POLD(10) IF (NR.EQ.0) EVOLD=Q(1) EVEN=Q(1) IF (EVOLD*EVEN.LT.0.0D0) THEN CALL ZEROEV(EVOLD,EVEN,EVENT,N,P,Q,POLD,QOLD,XOLD,X,XP) CALL POLINT(N,P,Q,POLD,QOLD,XOLD,X,XP,2,POLQ,POLP) C --- (POLP,POLQ) IS THE POINT OF THE POINCARE SECTION WRITE (6,88) XP,POLP,POLQ 88 FORMAT(1X,3F15.10) END IF EVOLD=EVEN DO I=1,N QOLD(I)=Q(I) POLD(I)=P(I) END DO RETURN END C SUBROUTINE ZEROEV(EVOLD,EVEN,EVENT,N,P,Q,POLD,QOLD,X0,X1,XM) C --- ROOT FINDING BY BISECTION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N),QOLD(10),POLD(10) XM0=X0 XM1=X1 33 XM=(XM1+XM0)/2.0D0 CALL POLINT(N,P,Q,POLD,QOLD,X0,X1,XM,1,POLQ,POLP) IF (EVOLD*POLQ.LT.0.0D0) THEN XM1=XM ELSE XM0=XM END IF IF (XM1-XM0.GE.1.0D-10) GOTO 33 RETURN END C SUBROUTINE POLINT(N,P,Q,POLD,QOLD,X0,X1,X,I,POLQ,POLP) C --- HERMITE INTERPOLATION IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N),QOLD(10),POLD(10) H=X1-X0 H2=H*H DIFF=(Q(I)-QOLD(I))/H XM0=X-X0 XM1=X-X1 P0DIF=POLD(I)-DIFF P1DIF=P(I)-DIFF POLQ=QOLD(I)+XM0*(DIFF+XM1*(P1DIF*XM0+P0DIF*XM1)/H2) POLP=DIFF+(XM0*(2*XM1+XM0)*P1DIF+XM1*(2*XM0+XM1)*P0DIF)/H2 RETURN END C SUBROUTINE HENON(N,X,Q,F,RPAR,IPAR) C --- RIGHT-HAND SIDE OF THE HENON-HEILES PROBLEM IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),F(N) F(1)=-Q(1)*(1+2*Q(2)) F(2)=-Q(2)*(1-Q(2))-Q(1)**2 RETURN END gnicodes/gni_comp.f100644 12122 12120 24526 7624125747 13414 0ustar hairermathC----------------------------------------------------------------------- SUBROUTINE GNI_COMP(BASIC,N,FCN,NSTEP,X,P,Q,XEND, & METHC,SOLFIX,IOUT,RPAR,IPAR) C----------------------------------------------------------------------- C VERSION OF SEPTEMBER 4,2002 C E-MAIL CONTACT ADDRESS : Ernst.Hairer@math.unige.ch C----------------------------------------------------------------------- C SYMMETRIC COMPOSITION OF SYMMETRIC METHODS C DESCRIBED IN SECTIONS II.4, V.3 OF THE BOOK: C C E. HAIRER, C. LUBICH, G. WANNER, GEOMETRIC NUMERICAL INTEGRATION, C STRUCTURE-PRESERVING ALGORITHMS FOR ODES. C SPRINGER SERIES IN COMPUT. MATH. 31, SPRINGER 2002. C C AND IN THE PUBLICATION C C E. HAIRER, M. HAIRER, GNI-CODES - MATLAB PROGRAMS FOR C GEOMETRIC NUMERICAL INTEGRATION. C C INPUT.. C BASIC NAME (EXTERNAL) OF BASIC INTEGRATION METHOD: C BASIC (N,X,P,Q,EP,EQ,FCN,HA,HB,HC,FIRST,LAST,RPAR,IPAR) C REAL*8 P(N),Q(N),EP(N),EQ(N) C LOGICAL FIRST,LAST C "STVERL" SOLVES SECOND ORDER ODES Q'' = F(X,Q) C C N DIMENSION OF Q AND F(Q) C C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING F(X,Q): C SUBROUTINE FCN(N,X,Q,F,RPAR,IPAR) C REAL*8 Q(N),F(N) C F(1)=... ETC. C C NSTEP NUMBER OF INTEGRATION STEPS C CONSTANT STEP SIZE, H=(XEND-X)/NSTEP C C X INITIAL X-VALUE C P(N) INITIAL VELOCITY VECTOR C Q(N) INITIAL POSITION VECTOR C XEND FINAL X-VALUE C C METHC ORDER AND STAGES OF THE COMPOSITION METHOD C E.G., 817 IS A METHOD OF ORDER 8 WITH 17 STAGES C POSSIBLE VALUES : 21,43,45,67,69,815,817,1033 C C SOLFIX NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE C NUMERICAL SOLUTION DURING INTEGRATION. C IF IOUT=1, IT IS CALLED AFTER EVERY STEP. C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) C DOUBLE PRECISION X,Y(N),CONT(LRC) C .... C SOLFIX FURNISHES THE SOLUTION "Q,P" AT THE NR-TH C GRID-POINT "X" (INITIAL VALUE FOR NR=0). C "XOLD" IS THE PRECEEDING GRID-POINT. C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN C IS SET <0, RETURN TO THE CALLING PROGRAM. C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLFIX: C IOUT=0: SUBROUTINE IS NEVER CALLED C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. C C RPAR(LR) REAL PARAMETER ARRAY; LR MUST BE AT LEAST LR=10 C RPAR(1),...,RPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C IPAR(LI) INTEGER PARAMETER ARRAY; LI MUST BE AT LEAST LI=10 C IPAR(1),...,IPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C C OUTPUT.. C P(N) SOLUTION (VELOCITY) AT XEND C Q(N) SOLUTION (POSITION) AT XEND C----------------------------------------------------------------------- C SOPHISTICATED SETTING OF PARAMETERS C----------------------------------------------------------------------- C NONE C----------------------------------------------------------------------- PARAMETER (NDGL=500,NSD=100) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N),Q(N),EP(NDGL),EQ(NDGL) DIMENSION GAM(NSD),HG(NSD),HGP(NSD) DIMENSION IPAR(*),RPAR(*) LOGICAL FIRST,LAST EXTERNAL FCN HH=(XEND-X)/NSTEP C --- PREPARE THE COEFFICIENTS IF (METHC.NE.21.AND.METHC.NE.43.AND.METHC.NE.45.AND.METHC.NE.67 & .AND.METHC.NE.69.AND.METHC.NE.815.AND.METHC.NE.817 & .AND.METHC.NE.1033.AND.METHC.NE.1034) THEN WRITE (6,*) ' NO METHOD AVAILABLE FOR METHC = ',METHC RETURN END IF CALL COMPCO (METHC,NGAM,GAM) DO IGAM=1,NGAM HG(IGAM)=HH*GAM(IGAM) END DO HGP(1)=HG(1)/2.0D0 DO IGAM=2,NGAM HGP(IGAM)=(HG(IGAM)+HG(IGAM-1))/2.0D0 END DO HGP(NGAM+1)=HG(NGAM)/2.0D0 HFL=HGP(1)+HGP(NGAM+1) IF (IOUT.NE.0) CALL SOLFIX (0,X,X,P,Q,N,IRTRN,RPAR,IPAR) FIRST=.TRUE. LAST=.FALSE. C --- THE BIG LOOP DO ISTEP=1,NSTEP X=X+HH DO IGAM=1,NGAM-1 CALL BASIC (N,X,P,Q,EP,EQ,FCN,HGP(IGAM),HG(IGAM), & HGP(IGAM+1),FIRST,LAST,RPAR,IPAR) FIRST=.FALSE. END DO IF (IOUT.NE.0) THEN LAST=.TRUE. CALL BASIC (N,X,P,Q,EP,EQ,FCN,HGP(NGAM),HG(NGAM), & HGP(NGAM+1),FIRST,LAST,RPAR,IPAR) CALL SOLFIX (ISTEP,X-HH,X,P,Q,N,IRTRN,RPAR,IPAR) IF (ISTEP.EQ.NSTEP.OR.IRTRN.LT.0) RETURN FIRST=.TRUE. LAST=.FALSE. ELSE IF (ISTEP.EQ.NSTEP) THEN LAST=.TRUE. CALL BASIC (N,X,P,Q,EP,EQ,FCN,HGP(NGAM),HG(NGAM), & HGP(NGAM+1),FIRST,LAST,RPAR,IPAR) RETURN ELSE CALL BASIC (N,X,P,Q,EP,EQ,FCN,HGP(NGAM),HG(NGAM), & HFL,FIRST,LAST,RPAR,IPAR) FIRST=.FALSE. END IF END IF END DO RETURN END C SUBROUTINE STVERL (N,X,P,Q,EP,EQ,FCN,HA,HB,HC, & FIRST,LAST,RPAR,IPAR) C--- STOERMER/VERLET FOR SECOND ORDER EQUATIONS C--- A SUBROUTINE FCN(N,X,Q,F) HAS TO PROVIDED PARAMETER (NDGL=500) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N),Q(N),EP(N),EQ(N),F(NDGL) DIMENSION IPAR(*),RPAR(*) LOGICAL FIRST,LAST IF (FIRST) THEN DO I=1,N TEMP=Q(I) EQ(I)=HA*P(I) Q(I)=TEMP+EQ(I) EQ(I)=EQ(I)+(TEMP-Q(I)) END DO END IF CALL FCN(N,X,Q,F,RPAR,IPAR) DO I=1,N TEMP=P(I) EP(I)=EP(I)+HB*F(I) P(I)=TEMP+EP(I) EP(I)=EP(I)+(TEMP-P(I)) END DO DO I=1,N TEMP=Q(I) EQ(I)=EQ(I)+HC*P(I) Q(I)=TEMP+EQ(I) EQ(I)=EQ(I)+(TEMP-Q(I)) END DO RETURN END C SUBROUTINE COMPCO (METHC,NGAM,GAM) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 GAM(100) IF (METHC.EQ.21) THEN C --- BASIC METHOD, NO COMPOSITION NGAM=1 GAM(1)=1.0D0 RETURN END IF IF (METHC.EQ.43) THEN C --- YOSHIDA (1990) NGAM=3 QR2=2.0D0**(1.0D0/3.0D0) GAM(1)=1.0D0/(2.0D0-QR2) GAM(2)=-QR2*GAM(1) END IF IF (METHC.EQ.45) THEN C --- SUZUKI (1990) NGAM=5 QR4=4.0D0**(1.0D0/3.0D0) GAM(1)=1.0D0/(4.0D0-QR4) GAM(2)=GAM(1) GAM(3)=-QR4*GAM(1) END IF IF (METHC.EQ.67) THEN C --- YOSHIDA (1990) NGAM=7 GAM(1)= 0.78451361047755726381949763D0 GAM(2)= 0.23557321335935813368479318D0 GAM(3)=-1.17767998417887100694641568D0 GAM(4)= 1.31518632068391121888424973D0 END IF IF (METHC.EQ.69) THEN C --- KAHAN & LI (1997) NGAM=9 GAM(1)= 0.39215846226558135393725811D0 GAM(2)= 0.33260288849639382848984234D0 GAM(3)=-0.70624607359230088534941023D0 GAM(4)= 0.08221272729243641010472095D0 GAM(5)= 0.79854399107577858563517766D0 END IF IF (METHC.EQ.815) THEN C --- SUZUKI & UMENO (1993) NGAM=15 GAM(1)= 0.74167036435061295344822780D0 GAM(2)=-0.40910082580003159399730010D0 GAM(3)= 0.19075471029623837995387626D0 GAM(4)=-0.57386247111608226665638773D0 GAM(5)= 0.29906418130365592384446354D0 GAM(6)= 0.33462491824529818378495798D0 GAM(7)= 0.31529309239676659663205666D0 GAM(8)=-0.79688793935291635401978884D0 END IF IF (METHC.EQ.817) THEN C --- KAHAN & LI (1997) NGAM=17 GAM(1)= 0.13020248308889008087881763D0 GAM(2)= 0.56116298177510838456196441D0 GAM(3)=-0.38947496264484728640807860D0 GAM(4)= 0.15884190655515560089621075D0 GAM(5)=-0.39590389413323757733623154D0 GAM(6)= 0.18453964097831570709183254D0 GAM(7)= 0.25837438768632204729397911D0 GAM(8)= 0.29501172360931029887096624D0 GAM(9)=-0.60550853383003451169892108D0 END IF IF (METHC.EQ.1033) THEN C --- HAIRER, LUBICH & WANNER (2002) NGAM=33 GAM( 1)= 0.09040619368607278492161150D0 GAM( 2)= 0.53591815953030120213784983D0 GAM( 3)= 0.35123257547493978187517736D0 GAM( 4)=-0.31116802097815835426086544D0 GAM( 5)=-0.52556314194263510431065549D0 GAM( 6)= 0.14447909410225247647345695D0 GAM( 7)= 0.02983588609748235818064083D0 GAM( 8)= 0.17786179923739805133592238D0 GAM( 9)= 0.09826906939341637652532377D0 GAM(10)= 0.46179986210411860873242126D0 GAM(11)=-0.33377845599881851314531820D0 GAM(12)= 0.07095684836524793621031152D0 GAM(13)= 0.23666960070126868771909819D0 GAM(14)=-0.49725977950660985445028388D0 GAM(15)=-0.30399616617237257346546356D0 GAM(16)= 0.05246957188100069574521612D0 GAM(17)= 0.44373380805019087955111365D0 END IF IF (METHC.EQ.1034) THEN C --- MARK SOFRONIOU & GIULIA SPALETTA (2003) NGAM=33 GAM( 1)= 0.070711261586085399079302771810203D0 GAM( 2)= 0.090342080937267568345577914389234D0 GAM( 3)= 0.14103133297152486103524322594476D0 GAM( 4)= 0.40206004554029767537357060971803D0 GAM( 5)=-0.30239722849131075165735249848238D0 GAM( 6)=-0.22462355658241460137093154363351D0 GAM( 7)= 0.061496988956063121940380707068411D0 GAM( 8)= 0.11346847775740802675296685287062D0 GAM( 9)= 0.23654672241381781124636015203490D0 GAM(10)= 0.27211409645898977643699556260890D0 GAM(11)= 0.076129418470277234386530906651024D0 GAM(12)=-0.18543093454238185309165565783301D0 GAM(13)=-0.46495110925607623804616342747217D0 GAM(14)= 0.10423014962104084592532590279051D0 GAM(15)= 0.13621181452383232935841998116651D0 GAM(16)=-0.27010275720513252644976102610064D0 GAM(17)= 0.48632639368142264147037913293721D0 END IF DO I=1,NGAM/2 GAM(NGAM+1-I)=GAM(I) END DO RETURN END gnicodes/comparison/ 40755 12122 12120 0 7535343655 13477 5ustar hairermathgnicodes/comparison/compar_dess.f100644 12122 12120 12306 7467660667 16275 0ustar hairermathcfggg compar_dess C this program uses GGG.f DIMENSION YDES0(1000),XDESF0(1000),XDEST0(1000) DIMENSION YDES4(1000),XDESF4(1000),XDEST4(1000) DIMENSION YDES6(1000),XDESF6(1000),XDEST6(1000) common/sizes/dim(4) data dim/5.1,5.0,0.8,0.6/ c call begin_ggg('compar_dess') C C --- OPEN DATA FILES ---- OPEN(11,FILE='kep_comp.dat') OPEN(12,FILE='kep_irk2.dat') OPEN(13,FILE='kep_lmm2.dat') OPEN(14,FILE='sol_comp.dat') OPEN(15,FILE='sol_irk2.dat') OPEN(16,FILE='sol_lmm2.dat') c --- READ DATA ----- DO IPROB=1,4,3 write(6,*)' iprob=',iprob IFILE=IPROB+11 DO I=1,500 READ(IFILE,*,END=21)YDES0(I),XDESF0(I),XDEST0(I) NPT0=I END DO 21 CONTINUE IF (IPROB.EQ.1) NPT0=NPT0 IF (IPROB.EQ.4) NPT0=NPT0-10 write (6,*) 'npt0 ',npt0 IFILE=IPROB+12 DO I=1,500 READ(IFILE,*,END=22)YDES4(I),XDESF4(I),XDEST4(I) NPT4=I END DO 22 CONTINUE IF (IPROB.EQ.1) NPT4=NPT4 IF (IPROB.EQ.4) NPT4=NPT4-11 write (6,*) 'npt4 ',npt4 IFILE=IPROB+10 DO I=1,500 READ(IFILE,*,END=23)YDES6(I),XDESF6(I),XDEST6(I) NPT6=I END DO 23 CONTINUE IF (IPROB.EQ.1) NPT6=NPT6 IF (IPROB.EQ.4) NPT6=NPT6 write (6,*) 'npt6 ',npt6 C ---- LOOP PICTURES ----- DO IPIC=1,2 call scale_char(0.75) call thick_pixel(3) C --- KEPLER --- IF(IPROB.EQ.1)THEN ymax=1.0 ymin=-12.2 IF (IPIC.EQ.1) THEN xmin=4.3 xmax=6.9 END IF IF (IPIC.EQ.2) THEN xmin=-1.5 xmax=0.5 END IF END IF C ---- SOLAR ---- IF(IPROB.EQ.4)THEN ymax=-1. ymin=-11.9 IF (IPIC.EQ.1) THEN xmin=3.5 xmax=5.5 END IF IF (IPIC.EQ.2) THEN xmin=-1.1 xmax= 0.6 END IF END IF call masog(xmin,xmax,ymin,ymax) call framg call logax_x(0.,YMIN,xmin,xmax,1,1) call logax_y(xmin,0.,ymin,ymax,3,1) C ---- DRAW LINES ---- call thick_pixel(7) IF(IPIC.EQ.1)THEN call select_line_type(2) CALL LING(XDESF0,YDES0,NPT0) call select_line_type(7) CALL LING(XDESF4,YDES4,NPT4) CALL LING(XDESF6,YDES6,NPT6) call thick_pixel(3) AMM=0.6 DO I=1,NPT0 CALL BUBBLE(XDESF0(I),YDES0(I),AMM,0.9,'circle') END DO DO I=1,NPT4 CALL BUBBLE(XDESF4(I),YDES4(I),AMM,0.5,'square') END DO DO I=1,NPT6 CALL BUBBLE(XDESF6(I),YDES6(I),AMM,0.5,'mhex') END DO END IF IF(IPIC.EQ.2)THEN call select_line_type(2) CALL LING(XDEST0,YDES0,NPT0) call select_line_type(7) CALL LING(XDEST4,YDES4,NPT4) CALL LING(XDEST6,YDES6,NPT6) call thick_pixel(3) AMM=0.6 DO I=1,NPT0 CALL BUBBLE(XDEST0(I),YDES0(I),AMM,0.9,'circle') END DO DO I=1,NPT4 CALL BUBBLE(XDEST4(I),YDES4(I),AMM,0.5,'square') END DO DO I=1,NPT6 CALL BUBBLE(XDEST6(I),YDES6(I),AMM,0.5,'mhex') END DO END IF C ---- TEXTES ------ CALL TEXT_LATEX_RG(XMIN,(YMAX+2*YMIN)/3.,0.5,3.0,90.,1. &,'error') IF(IPROB.EQ.1)THEN CALL TEXT_LATEX((XMIN+XMAX)/2.,YMAX,0.3,1.4 &,'\normalsize Kepler') END IF IF(IPROB.EQ.4)THEN CALL TEXT_LATEX((XMIN+XMAX)/2.,YMAX,-0.4,1.4 &,'\normalsize Solar') END IF IF(IPIC.EQ.1)THEN CALL TEXT_LATEX((2*XMIN+XMAX)/3.,YMIN,0.8,-1.1, &'fcn.\,eval.') END IF IF(IPIC.EQ.2)THEN CALL TEXT_LATEX((2*XMIN+XMAX)/3.,YMIN,0.8,-0.6, &'cpu time') END IF C ----- FIRST PICTURE --- IF(IPROB.EQ.1.AND.IPIC.EQ.1)THEN CALL TEXT_RELAD_LATEX(0.66,0.33,0.5,0.5,'comp') CALL TEXT_RELAD_LATEX(0.4,0.2,0.5,0.5,'irk2') CALL TEXT_RELAD_LATEX(0.25,0.45,0.5,0.5,'lmm2') END IF C ----- SECOND PICTURE --- IF(IPROB.EQ.1.AND.IPIC.EQ.2)THEN CALL TEXT_RELAD_LATEX(0.6,0.56,0.5,0.5,'comp') CALL TEXT_RELAD_LATEX(0.4,0.36,0.5,0.5,'irk2') CALL TEXT_RELAD_LATEX(0.15,0.6,0.5,0.5,'lmm2') END IF C ----- THIRD PICTURE --- IF(IPROB.EQ.4.AND.IPIC.EQ.1)THEN CALL TEXT_RELAD_LATEX(0.82,0.40,0.5,0.5,'comp') CALL TEXT_RELAD_LATEX(0.52,0.36,0.5,0.5,'irk2') CALL TEXT_RELAD_LATEX(0.32,0.6,0.5,0.5,'lmm2') END IF C ----- FOURTH PICTURE --- IF(IPROB.EQ.4.AND.IPIC.EQ.2)THEN CALL TEXT_RELAD_LATEX(0.44,0.58,0.5,0.5,'comp') CALL TEXT_RELAD_LATEX(0.68,0.75,0.5,0.5,'irk2') CALL TEXT_RELAD_LATEX(0.32,0.25,0.5,0.5,'lmm2') END IF END DO END DO call greater_boundingbox (0.7,-0.1,0.3,0.) call end_ggg STOP END gnicodes/comparison/cal_comp.f100644 12122 12120 4036 7467215116 15515 0ustar hairermathccfehop cal_comp include '../gni_comp.f' include '../problem.f' IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX,STVERL C --- OPEN DATA FILES ---- OPEN(12,FILE='kep_comp.dat') OPEN(15,FILE='sol_comp.dat') C --- DO IPROB=1,4,3 IFILE=IPROB+11 REWIND (IFILE) METH=817 WRITE (6,*) ' METHOD, PROBLEM ',METH,IPROB IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 H=800.0D0 IF (IPROB.EQ.1) H=0.5D0 DO IH=1,40 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) IF (IPROB.EQ.1) XEND=200*XEND NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=0 CALL CPU_TIME(TIME0) CALL GNI_COMP(STVERL,N,EQUA,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR 90 FORMAT (1X,I8,F7.2,1X,2E24.16) c XDAT=LOG10(max(1.d-16,ERR)) ANFCN=NFCN YDATF=LOG10(ANFCN) YDATT=LOG10(max(1.e-16,TIME1-TIME0)) WRITE(IFILE,956)XDAT,YDATF,YDATT c WRITE(6,957)XDAT,YDATF,YDATT,NSTEP 956 FORMAT(3F11.6) 957 FORMAT(3F11.6,I10) c H=H*0.9D0 END DO CLOSE (IFILE) END DO STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF RETURN END gnicodes/comparison/cal_irk2.f100644 12122 12120 4017 7467214710 15424 0ustar hairermathccfehop cal_irk2 include '../gni_irk2.f' include '../problem.f' IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX C --- OPEN DATA FILES ---- OPEN(12,FILE='kep_irk2.dat') OPEN(15,FILE='sol_irk2.dat') DO IPROB=1,4,3 IFILE=IPROB+11 REWIND (IFILE) METH=4 WRITE (6,*) ' METHOD, PROBLEM ',METH,IPROB IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 H=800.0D0 IF (IPROB.EQ.1) H=0.4D0 DO IH=1,40 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) IF (IPROB.EQ.1) XEND=200*XEND NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=0 CALL CPU_TIME(TIME0) CALL GNI_IRK2(N,EQUA,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR 90 FORMAT (1X,I8,F7.2,1X,2E24.16) c XDAT=LOG10(max(1.d-16,ERR)) ANFCN=NFCN YDATF=LOG10(ANFCN) YDATT=LOG10(max(1.e-16,TIME1-TIME0)) WRITE(IFILE,956)XDAT,YDATF,YDATT c WRITE(6,957)XDAT,YDATF,YDATT,NSTEP 956 FORMAT(3F11.6) 957 FORMAT(3F11.6,I10) c H=H*0.9D0 END DO CLOSE (IFILE) END DO STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF RETURN END gnicodes/comparison/cal_lmm2.f100644 12122 12120 4060 7467215140 15420 0ustar hairermathccfehop cal_lmm2 include '../gni_irk2.f' include '../gni_lmm2.f' include '../problem.f' IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX C --- OPEN DATA FILES ---- OPEN(12,FILE='kep_lmm2.dat') OPEN(15,FILE='sol_lmm2.dat') DO IPROB=1,4,3 IFILE=IPROB+11 REWIND (IFILE) METH=803 WRITE (6,*) ' METHOD, PROBLEM ',METH,IPROB IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 H=200.0D0 IF (IPROB.EQ.1) H=0.05D0 DO IH=1,40 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) IF (IPROB.EQ.1) XEND=200*XEND NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=0 CALL CPU_TIME(TIME0) CALL GNI_LMM2(N,EQUA,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR 90 FORMAT (1X,I8,F7.2,1X,2E24.16) c XDAT=LOG10(max(1.d-16,ERR)) ANFCN=NFCN YDATF=LOG10(ANFCN) YDATT=LOG10(max(1.e-16,TIME1-TIME0)) WRITE(IFILE,956)XDAT,YDATF,YDATT c WRITE(6,957)XDAT,YDATF,YDATT,NSTEP 956 FORMAT(3F11.6) 957 FORMAT(3F11.6,I10) c H=H*0.9D0 END DO CLOSE (IFILE) END DO STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF RETURN END gnicodes/gni_irk2.f100644 12122 12120 43226 7570427524 13321 0ustar hairermathC----------------------------------------------------------------------- SUBROUTINE GNI_IRK2(N,FCN,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) C----------------------------------------------------------------------- C VERSION OF SEPTEMBER 4,2002 C E-MAIL CONTACT ADDRESS : Ernst.Hairer@math.unige.ch C----------------------------------------------------------------------- C SOLVES SECOND ORDER ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C Q'' = F(X,Q) C BASED ON THE SYMPLECTIC AND SYMMETRIC GAUSS (IRK) METHODS C DESCRIBED IN SECTIONS II.1, VIII.6 OF THE BOOK: C C E. HAIRER, C. LUBICH, G. WANNER, GEOMETRIC NUMERICAL INTEGRATION, C STRUCTURE-PRESERVING ALGORITHMS FOR ODES. C SPRINGER SERIES IN COMPUT. MATH. 31, SPRINGER 2002. C C AND IN THE PUBLICATION C C E. HAIRER, M. HAIRER, GNI-CODES - MATLAB PROGRAMS FOR C GEOMETRIC NUMERICAL INTEGRATION. C C INPUT.. C N DIMENSION OF Q AND F(X,Q) C C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING F(X,Q): C SUBROUTINE FCN(N,X,Q,F,RPAR,IPAR) C REAL*8 Q(N),F(N) C F(1)=... ETC. C C NSTEP NUMBER OF INTEGRATION STEPS C CONSTANT STEP SIZE, H=(XEND-X)/NSTEP C C X INITIAL X-VALUE C P(N) INITIAL VELOCITY VECTOR C Q(N) INITIAL POSITION VECTOR C XEND FINAL X-VALUE C C METH NUMBER OF STAGES OF THE GAUSS METHOD C FOR THE MOMENT ONLY POSSIBLE VALUES: 2,4,6. C C SOLFIX NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE C NUMERICAL SOLUTION DURING INTEGRATION. C IF IOUT=1, IT IS CALLED AFTER EVERY STEP. C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) C DOUBLE PRECISION X,Y(N),CONT(LRC) C .... C SOLFIX FURNISHES THE SOLUTION "Q,P" AT THE NR-TH C GRID-POINT "X" (INITIAL VALUE FOR NR=0). C "XOLD" IS THE PRECEEDING GRID-POINT. C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN C IS SET <0, RETURN TO THE CALLING PROGRAM. C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLFIX: C IOUT=0: SUBROUTINE IS NEVER CALLED C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. C C RPAR(LR) REAL PARAMETER ARRAY; LR MUST BE AT LEAST LR=10 C RPAR(1),...,RPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C IPAR(LI) INTEGER PARAMETER ARRAY; LI MUST BE AT LEAST LI=10 C IPAR(1),...,IPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C C OUTPUT.. C P(N) SOLUTION (VELOCITY) AT XEND C Q(N) SOLUTION (POSITION) AT XEND C----------------------------------------------------------------------- C SOPHISTICATED SETTING OF PARAMETERS C----------------------------------------------------------------------- C RPAR(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. C IPAR(1) NITMAX, MAXIMAL NUMER OF FIXED POINT ITERAT., DEFAULT 50 C----------------------------------------------------------------------- PARAMETER (NDGL=500,NSD=6,NMD=3) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(NDGL*NSD),EQ(NDGL),EP(NDGL),YH(NDGL),QQ(NDGL) DIMENSION C(NSD),AA(NSD,NSD),E(NSD,NSD+NMD),B(NSD),BC(NSD) DIMENSION SM(NMD),AMF(NMD,NSD+NMD),AM(NSD+NMD) DIMENSION Q(N),P(N),FS(NDGL),PS(NDGL),ZQ(NDGL,NSD) DIMENSION IPAR(*),RPAR(*) EXTERNAL FCN C -------- UROUND SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0 IF (RPAR(1).EQ.0.0D0) THEN UROUND=1.0D-16 ELSE UROUND=RPAR(1) END IF C -------- NITMAX, MAXIMAL NUMER OF FIXED POINT ITERATIONS IF (IPAR(1).EQ.0) THEN NITMAX=50 ELSE NITMAX=IPAR(1) END IF C -------- NS=METH H=(XEND-X)/NSTEP CALL COEFG(NS,C,B,BC,NSD,AA,E,NM,SM,NMD,AM,H) IF (IOUT.NE.0) CALL SOLFIX (0,X,X,P,Q,N,IRTRN,RPAR,IPAR) CALL FCN(N,X,Q,FS,RPAR,IPAR) DO IS=1,NS FAC=C(IS)**2/2 DO I=1,N ZQ(I,IS)=C(IS)*P(I)+FAC*FS(I) END DO END DO DO I=1,N PS(I)=P(I) EQ(I)=0.0D0 EP(I)=0.0D0 END DO C --- LOOP FOR THE ITERATIONS DO ISTEP=1,NSTEP IF (ISTEP.GT.1) CALL STARTB (FCN,N,X,P,Q,NS,NDGL,FS,PS, & ZQ,NSD,E,YH,NM,SM,AM,F,C,RPAR,IPAR) C --- FIXED POINT ITERATION NITER=0 DYNOLD=0.0D0 40 CONTINUE CALL RKNITE(FCN,N,NS,X,Q,P,NSD,AA,C,NDGL,QQ,ZQ,F,DYNO,RPAR,IPAR) NITER=NITER+1 IF (DYNOLD.LT.DYNO.AND.DYNO.LT.10*UROUND) GOTO 50 IF (NITER.GE.NITMAX) THEN WRITE (6,*) 'NO CONVERGENCE OF ITERATION',DYNO IF (DYNO.GT.UROUND*1.D6) RETURN END IF IF (DYNO.GT.0.1D0*UROUND) GOTO 40 50 CONTINUE C --- UPDATE OF THE SOLUTION X=X+H DO I=1,N EQI=EQ(I) QI=Q(I) AY=QI EQI=EQI+H*P(I) QI=AY+EQI EQI=EQI+(AY-QI) DO IS=1,NS AY=QI EQI=EQI+F(I+(IS-1)*N)*BC(IS) QI=AY+EQI EQI=EQI+(AY-QI) END DO AY=Q(I) Q(I)=QI EQ(I)=EQI EPI=EP(I) PI=P(I) DO IS=1,NS/2 AY=PI EPI=EPI+(F(I+(IS-1)*N)+F(I+(NS-IS)*N))*B(IS) PI=AY+EPI EPI=EPI+(AY-PI) END DO P(I)=PI EP(I)=EPI END DO IF (IOUT.NE.0) CALL SOLFIX (ISTEP,X-H,X,P,Q,N,IRTRN,RPAR,IPAR) END DO RETURN END C SUBROUTINE STARTB (FCN,N,X,P,Q,NS,NDGL,FS,PS,ZQ,NSD,E,YH, & NM,SM,AM,F,C,RPAR,IPAR) C ---------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ZQ(NDGL,NS),E(NSD,NS+NM),PS(N),F(N*NS),C(NS) REAL*8 AM(NS+NM),SM(NM),P(N),Q(N),YH(NDGL),FS(NDGL) DIMENSION IPAR(*),RPAR(*) NS1=NS+1 NS2=NS+2 NSM=NS+NM DO I=1,N SAV=0.0D0 DO JS=1,NS SAV=SAV+AM(JS)*ZQ(I,JS) END DO YH(I)=SAV+AM(NS1)*PS(I)+AM(NS2)*P(I)+Q(I) DO IS=1,NS SAV=0.0D0 DO JS=1,NS SAV=SAV+E(IS,JS)*F(I+(JS-1)*N) END DO ZQ(I,IS)=SAV+E(IS,NS1)*FS(I) END DO END DO C CALL FCN(N,X+H,Q,FS,RPAR,IPAR) CALL FCN(N,X+H*SM(NM),YH,F,RPAR,IPAR) DO I=1,N PS(I)=P(I) DO IS=1,NS ZQIIS=ZQ(I,IS)+E(IS,NS2)*FS(I)+E(IS,NSM)*F(I) ZQ(I,IS)=ZQIIS+C(IS)*P(I) END DO END DO RETURN END C SUBROUTINE RKNITE(FCN,N,NS,X,Q,P,NSD,AA,C,NDGL,QQ,ZQ, & F,DYNO,RPAR,IPAR) C ---------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) REAL*8 ZQ(NDGL,NS),C(NS),Q(N),P(N) REAL*8 AA(NSD,NS),F(NDGL*NS),QQ(NDGL) DIMENSION IPAR(*),RPAR(*) EXTERNAL FCN C --- DO JS=1,NS DO J=1,N QQ(J)=Q(J)+ZQ(J,JS) END DO CALL FCN(N,X+C(JS),QQ,F(1+(JS-1)*N),RPAR,IPAR) END DO C --- DYNO=0.D0 DO I=1,N DNOM=MAX(1.0D-1,ABS(Q(I))) DO IS=1,NS SUM=C(IS)*P(I) DO JS=1,NS SUM=SUM+AA(IS,JS)*F(I+(JS-1)*N) END DO DYNO=DYNO+((SUM-ZQ(I,IS))/DNOM)**2 ZQ(I,IS)=SUM END DO END DO DYNO=DSQRT(DYNO/(NS*N)) C --- RETURN END C SUBROUTINE COEFG(NS,C,B,BC,NSD,AA,E,NM,SM,NMD,AM,HSTEP) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 C(NS),AA(NSD,NS),E(NSD,NS+NM),B(NS) REAL*8 AM(NS+NM),SM(NM),BC(NS) NM=3 IF (NS.EQ.2) THEN C(1)= 0.21132486540518711775D+00 C(2)= 0.78867513459481288225D+00 B(1)= 0.50000000000000000000D+00 B(2)= 0.50000000000000000000D+00 BC(1)= 0.39433756729740644113D+00 BC(2)= 0.10566243270259355887D+00 AA(1,1)= 0.41666666666666666667D-01 AA(1,2)=-0.19337567297406441127D-01 AA(2,1)= 0.26933756729740644113D+00 AA(2,2)= 0.41666666666666666667D-01 E(1,1)=-0.28457905077110526160D-02 E(1,2)=-0.63850024471784160410D-01 E(1,3)= 0.48526095198694517563D-02 E(1,4)= 0.11305688530429939012D+00 E(1,5)=-0.28884580475413403312D-01 E(2,1)= 0.41122751744511433137D-01 E(2,2)=-0.18654814888622834132D+00 E(2,3)=-0.18110185277445209332D-01 E(2,4)= 0.36674109449368040786D+00 E(2,5)= 0.10779872188955481745D+00 SM(1)= 0.00000000000000000000D+00 SM(2)= 0.10000000000000000000D+01 SM(3)= 0.16000000000000000000D+01 AM(1)= 0.25279583039343438291D+02 AM(2)=-0.86907830393434382912D+01 AM(3)=-0.80640000000000000000D+00 AM(4)= 0.29184000000000000000D+01 AM(5)= 0.00000000000000000000D+00 END IF IF (NS.EQ.4) THEN C(1)= 0.69431844202973712388D-01 C(2)= 0.33000947820757186760D+00 C(3)= 0.66999052179242813240D+00 C(4)= 0.93056815579702628761D+00 B(1)= 0.17392742256872692869D+00 B(2)= 0.32607257743127307131D+00 B(3)= 0.32607257743127307131D+00 B(4)= 0.17392742256872692869D+00 BC(1)= 0.16185132086231030665D+00 BC(2)= 0.21846553629538057030D+00 BC(3)= 0.10760704113589250101D+00 BC(4)= 0.12076101706416622036D-01 AA(1,1)= 0.40381914508467311298D-02 AA(1,2)=-0.32958609449446961650D-02 AA(1,3)= 0.26447829520668538006D-02 AA(1,4)=-0.97672296325588161023D-03 AA(2,1)= 0.43563580902396261254D-01 AA(2,2)= 0.13818951406296126013D-01 AA(2,3)=-0.43401341944349953440D-02 AA(2,4)= 0.14107297391595337720D-02 AA(3,1)= 0.10586435263357640763D+00 AA(3,2)= 0.10651836096505307395D+00 AA(3,3)= 0.13818951406296126013D-01 AA(3,4)=-0.17580153590805494993D-02 AA(4,1)= 0.14879849619263780300D+00 AA(4,2)= 0.19847049885237718995D+00 AA(4,3)= 0.81671359795877570687D-01 AA(4,4)= 0.40381914508467311298D-02 E(1,1)=-0.21272768296134340207D-01 E(1,2)= 0.11059138674756969912D-01 E(1,3)= 0.38999255049973564023D-02 E(1,4)=-0.43986226789008967612D-01 E(1,5)= 0.13581590305438849621D-01 E(1,6)= 0.39922421675314269059D-01 E(1,7)=-0.79369058065113002021D-03 E(2,1)=-0.75671119283734809953D-02 E(2,2)= 0.10209394000843457002D-01 E(2,3)=-0.12880197817980892596D-01 E(2,4)=-0.56381316813776501277D-01 E(2,5)= 0.37440782682669799960D-02 E(2,6)= 0.11522469441011273193D+00 E(2,7)= 0.21035877343246316334D-02 E(3,1)=-0.39890571772473709759D+00 E(3,2)= 0.26819725655216894347D+00 E(3,3)=-0.82551711648854471247D-01 E(3,4)=-0.85516559106259630212D+00 E(3,5)= 0.24433810515772642570D+00 E(3,6)= 0.10234155624049009806D+01 E(3,7)= 0.25115745967236579242D-01 E(4,1)=-0.40964796048052939224D+00 E(4,2)= 0.29949323098224574487D+00 E(4,3)=-0.13867460566101912494D+00 E(4,4)=-0.98859300714628940382D+00 E(4,5)= 0.24671351779481625627D+00 E(4,6)= 0.12912760231350872304D+01 E(4,7)= 0.13241134766742798418D+00 SM(1)= 0.00000000000000000000D+00 SM(2)= 0.10000000000000000000D+01 SM(3)= 0.16500000000000000000D+01 AM(1)= 0.10806374869244001787D+04 AM(2)=-0.66008818661284690206D+03 AM(3)= 0.61810154357557529566D+03 AM(4)=-0.31341427826212857229D+03 AM(5)=-0.10187174765625000000D+02 AM(6)= 0.31173050390625000000D+02 AM(7)= 0.00000000000000000000D+00 END IF IF (NS.EQ.6) THEN C(1)= 0.33765242898423986094D-01 C(2)= 0.16939530676686774317D+00 C(3)= 0.38069040695840154568D+00 C(4)= 0.61930959304159845432D+00 C(5)= 0.83060469323313225683D+00 C(6)= 0.96623475710157601391D+00 B(1)= 0.85662246189585172520D-01 B(2)= 0.18038078652406930378D+00 B(3)= 0.23395696728634552369D+00 B(4)= 0.23395696728634552369D+00 B(5)= 0.18038078652406930378D+00 B(6)= 0.85662246189585172520D-01 BC(1)= 0.82769839639769234611D-01 BC(2)= 0.14982512785597570103D+00 BC(3)= 0.14489179419935320895D+00 BC(4)= 0.89065173086992314743D-01 BC(5)= 0.30555658668093602753D-01 BC(6)= 0.28924065498159379092D-02 AA(1,1)= 0.90625420195651151857D-03 AA(1,2)=-0.72859711612531400024D-03 AA(1,3)= 0.79102695861167691135D-03 AA(1,4)=-0.70675390218535384182D-03 AA(1,5)= 0.45647714224056921122D-03 AA(1,6)=-0.14836147050330408643D-03 AA(2,1)= 0.11272367531794365387D-01 AA(2,2)= 0.39083482447840698486D-02 AA(2,3)=-0.14724868010943911900D-02 AA(2,4)= 0.10992669056588431310D-02 AA(2,5)=-0.67689040729401428165D-03 AA(2,6)= 0.21677950347174141516D-03 AA(3,1)= 0.30008019623627547434D-01 AA(3,2)= 0.36978289259468146662D-01 AA(3,3)= 0.65490339168957822692D-02 AA(3,4)=-0.16615098173008262274D-02 AA(3,5)= 0.84753461862041607649D-03 AA(3,6)=-0.25877462623437421721D-03 AA(4,1)= 0.49900269650650898941D-01 AA(4,2)= 0.82003427445271620462D-01 AA(4,3)= 0.54165111295060067982D-01 AA(4,4)= 0.65490339168957822692D-02 AA(4,5)=-0.11352871017627472322D-02 AA(4,6)= 0.28963081055952389031D-03 AA(5,1)= 0.68475836671617248304D-01 AA(5,2)= 0.11859257878058808400D+00 AA(5,3)= 0.10635984886129551097D+00 AA(5,4)= 0.47961474042181382443D-01 AA(5,5)= 0.39083482447840698486D-02 AA(5,6)=-0.34600839001342442657D-03 AA(6,1)= 0.79729071619449992615D-01 AA(6,2)= 0.14419100392702230613D+00 AA(6,3)= 0.13628542646896576408D+00 AA(6,4)= 0.81956586217401900627D-01 AA(6,5)= 0.23736460480774324642D-01 AA(6,6)= 0.90625420195651151857D-03 E(1,1)=-0.16761132335280609813D-01 E(1,2)= 0.10201050166615899799D-01 E(1,3)=-0.58593121685075943100D-02 E(1,4)=-0.11907383391366998251D-03 E(1,5)= 0.10615611118132982241D-01 E(1,6)=-0.30692054230989138447D-01 E(1,7)= 0.10615182045216224925D-01 E(1,8)= 0.22586707045496892369D-01 E(1,9)=-0.16931992776201068110D-04 E(2,1)= 0.10671755276327262128D-01 E(2,2)=-0.51098203653251450913D-02 E(2,3)= 0.16062647299186369205D-03 E(2,4)= 0.64818802653621866868D-02 E(2,5)=-0.12132386914873895089D-01 E(2,6)=-0.99709737725909584834D-02 E(2,7)=-0.70287093442894942752D-02 E(2,8)= 0.31243249755879001843D-01 E(2,9)= 0.31763603839792897936D-04 E(3,1)=-0.40875203230945019464D+00 E(3,2)= 0.28214948905763253599D+00 E(3,3)=-0.22612660499718519054D+00 E(3,4)= 0.13640993962985420478D+00 E(3,5)= 0.15888529591697266925D+00 E(3,6)=-0.11667863471317749710D+01 E(3,7)= 0.25224964119340060668D+00 E(3,8)= 0.10440940643938620983D+01 E(3,9)= 0.33914722176493324285D-03 E(4,1)=-0.29437531285359759661D+01 E(4,2)= 0.20017220470127690267D+01 E(4,3)=-0.15383035791443948798D+01 E(4,4)= 0.78114323215109899716D+00 E(4,5)= 0.13930345104184182146D+01 E(4,6)=-0.75958281612589849630D+01 E(4,7)= 0.18220129530415584951D+01 E(4,8)= 0.62663163493155487560D+01 E(4,9)= 0.54279630166374655267D-02 E(5,1)=-0.79572842006457093076D+01 E(5,2)= 0.53527892762707449170D+01 E(5,3)=-0.40049139768467199697D+01 E(5,4)= 0.18326058141135591515D+01 E(5,5)= 0.39753886181058367500D+01 E(5,6)=-0.19423696478604790213D+02 E(5,7)= 0.49362128400107292627D+01 E(5,8)= 0.15601708062381928560D+02 E(5,9)= 0.32142123424873719685D-01 E(6,1)=-0.78463118056075171475D+01 E(6,2)= 0.53580869574441241664D+01 E(6,3)=-0.41476905275607763365D+01 E(6,4)= 0.21275912797813913113D+01 E(6,5)= 0.37642416878253538582D+01 E(6,6)=-0.20329681631523484613D+02 E(6,7)= 0.48515418060343387549D+01 E(6,8)= 0.16604467346259915039D+02 E(6,9)= 0.84559690262225766975D-01 SM(1)= 0.00000000000000000000D+00 SM(2)= 0.10000000000000000000D+01 SM(3)= 0.17500000000000000000D+01 AM(1)= 0.58080578375796358720D+05 AM(2)=-0.33214989339522861968D+05 AM(3)= 0.28376088288312020853D+05 AM(4)=-0.27923430684614999462D+05 AM(5)= 0.29743005589491042677D+05 AM(6)=-0.15525927919158826444D+05 AM(7)=-0.27700591278076171875D+03 AM(8)= 0.73086943817138671875D+03 AM(9)= 0.00000000000000000000D+00 END IF C HSTEP2=HSTEP*HSTEP DO IS=1,NS B(IS)=HSTEP*B(IS) BC(IS)=HSTEP2*BC(IS) C(IS)=HSTEP*C(IS) DO JS=1,NS AA(IS,JS)=HSTEP2*AA(IS,JS) E(IS,JS)=HSTEP2*E(IS,JS) END DO END DO DO IM=1,NM DO IS=1,NS E(IS,NS+IM)=HSTEP2*E(IS,NS+IM) END DO AM(NS+IM)=HSTEP*AM(NS+IM) END DO C RETURN END gnicodes/gni_lmm2.f100644 12122 12120 24211 7561742233 13307 0ustar hairermathC----------------------------------------------------------------------- SUBROUTINE GNI_LMM2(N,FCN,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) C----------------------------------------------------------------------- C VERSION OF NOVEMBER 3,2002 C E-MAIL CONTACT ADDRESS : Ernst.Hairer@math.unige.ch C----------------------------------------------------------------------- C SOLVES SECOND ORDER ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C Q'' = F(X,Q) C BASED ON SYMMETRIC LINEAR MULTISTEP METHODS C DESCRIBED IN CHAPTER XIV OF THE BOOK: C C E. HAIRER, C. LUBICH, G. WANNER, GEOMETRIC NUMERICAL INTEGRATION, C STRUCTURE-PRESERVING ALGORITHMS FOR ODES. C SPRINGER SERIES IN COMPUT. MATH. 31, SPRINGER 2002. C C AND IN THE PUBLICATION C C E. HAIRER, M. HAIRER, GNI-CODES - MATLAB PROGRAMS FOR C GEOMETRIC NUMERICAL INTEGRATION. C C INPUT.. C N DIMENSION OF Q AND F(X,Q) C C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING F(X,Q): C SUBROUTINE FCN(N,X,Q,F,RPAR,IPAR) C REAL*8 Q(N),F(N) C F(1)=... ETC. C C NSTEP NUMBER OF INTEGRATION STEPS C CONSTANT STEP SIZE, H=(XEND-X)/NSTEP C C X INITIAL X-VALUE C P(N) INITIAL VELOCITY VECTOR C Q(N) INITIAL POSITION VECTOR C XEND FINAL X-VALUE C C METH NUMBER OF METHOD (ORDER = METH/100) C ORDER 2 : 201 (STOERMER/VERLET, MOD.INIT.VELOCITY) C ORDER 4 : 401 C ORDER 8 : 801,802,803 C C SOLFIX NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE C NUMERICAL SOLUTION DURING INTEGRATION. C IF IOUT=1, IT IS CALLED AFTER EVERY STEP. C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) C DOUBLE PRECISION X,P(N),Q(N) C .... C SOLFIX FURNISHES THE SOLUTION "Q,P" AT THE NR-TH C GRID-POINT "X" (INITIAL VALUE FOR NR=0). C "XOLD" IS THE PRECEEDING GRID-POINT. C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN C IS SET <0, RETURN TO THE CALLING PROGRAM. C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLFIX: C IOUT=0: SUBROUTINE IS NEVER CALLED C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. C C RPAR(LR) REAL PARAMETER ARRAY; LR MUST BE AT LEAST LR=10 C RPAR(1),...,RPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C IPAR(LI) INTEGER PARAMETER ARRAY; LI MUST BE AT LEAST LI=10 C IPAR(1),...,IPAR(10) SERVE AS PARAMETERS FOR C THE CODE. FURTHER VALUES CAN BE USED FOR DEFINING C PARAMETERS IN THE PROBLEM C C OUTPUT.. C P(N) SOLUTION (VELOCITY) AT XEND C Q(N) SOLUTION (POSITION) AT XEND C----------------------------------------------------------------------- C SOPHISTICATED SETTING OF PARAMETERS C----------------------------------------------------------------------- C NONE C----------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (ND=500) DOUBLE PRECISION Q0(ND),U(ND),Z0(ND),P(N),Q(N) DOUBLE PRECISION F1(ND),F2(ND),F3(ND),F4(ND),F5(ND) DOUBLE PRECISION F6(ND),F7(ND),F8(ND),F9(ND) DOUBLE PRECISION Q1(ND),Q2(ND),Q3(ND),Q4(ND),Q5(ND) DOUBLE PRECISION Q6(ND),Q7(ND),Q8(ND),Q9(ND) DOUBLE PRECISION Z1(ND),Z2(ND),Z3(ND),Z4(ND),Z5(ND) DOUBLE PRECISION Z6(ND),Z7(ND),Z8(ND),Z9(ND) DIMENSION IPAR(*),RPAR(*) EXTERNAL FCN,SOLFIX H=(XEND-X)/NSTEP KK=METH/100 CALL SYCOE(METH,A1,A2,A3,A4,B1,B2,B3,B4,H,HM) IF (METH.LE.0) THEN METH=-METH WRITE (6,*) ' NO COEFFICIENTS AVAILABLE FOR METHOD ',METH RETURN END IF NR=0 IRTRN=0 DO I=1,N Q0(I)=Q(I) END DO IF (IOUT.GE.1) CALL SOLFIX (NR,X,X,P,Q,N,IRTRN,RPAR,IPAR) XRK=X C DO IK=1,KK-1 XDEND=XRK+H NSTRK=1 METHRK=4 IOUTRK=0 CALL GNI_IRK2(N,FCN,NSTRK,XRK,P,Q,XDEND, & METHRK,SOLFIX,IOUTRK,RPAR,IPAR) XRK=XDEND IF (IK.EQ.1) CALL PREPST(N,FCN,H,Q,Q0,Q1,Z0,F1,RPAR,IPAR) IF (IK.EQ.2) CALL PREPST(N,FCN,H,Q,Q1,Q2,Z1,F2,RPAR,IPAR) IF (IK.EQ.3) CALL PREPST(N,FCN,H,Q,Q2,Q3,Z2,F3,RPAR,IPAR) IF (IK.EQ.4) CALL PREPST(N,FCN,H,Q,Q3,Q4,Z3,F4,RPAR,IPAR) IF (IK.EQ.5) CALL PREPST(N,FCN,H,Q,Q4,Q5,Z4,F5,RPAR,IPAR) IF (IK.EQ.6) CALL PREPST(N,FCN,H,Q,Q5,Q6,Z5,F6,RPAR,IPAR) IF (IK.EQ.7) CALL PREPST(N,FCN,H,Q,Q6,Q7,Z6,F7,RPAR,IPAR) IF (IK.EQ.8) CALL PREPST(N,FCN,H,Q,Q7,Q8,Z7,F8,RPAR,IPAR) IF (IK.EQ.9) CALL PREPST(N,FCN,H,Q,Q8,Q9,Z8,F9,RPAR,IPAR) IF (IK.LT.KK/2) THEN NR=NR+1 XOLD=X X=XDEND IF (IOUT.GE.1) CALL SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) END IF END DO C IF (KK.EQ.2) THEN 22 CONTINUE CALL FCN(N,X,Q1,F1,RPAR,IPAR) DO I=1,N Z0(I)=Z0(I)+H*F1(I) Q2(I)=Q1(I)+H*Z0(I) END DO DO I=1,N Q(I)=Q1(I) P(I)=(Q2(I)-Q0(I))/(2*H) Q0(I)=Q1(I) Q1(I)=Q2(I) END DO XOLD=X X=X+H NR=NR+1 IF (IOUT.GE.1) CALL SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IF (NR.LT.NSTEP.AND.IRTRN.GE.0) GOTO 22 END IF C IF (KK.EQ.4) THEN DO I=1,N U(I)=A1*(Z0(I)+Z2(I))+A2*Z1(I) END DO 44 CONTINUE CALL FCN(N,X,Q3,F3,RPAR,IPAR) DO I=1,N U(I)=U(I)+(B1*(F1(I)+F3(I))+B2*F2(I)) SUM=A1*Z1(I)+A2*Z2(I) Z3(I)=U(I)-SUM Q4(I)=Q3(I)+H*Z3(I) END DO DO I=1,N Q(I)=Q2(I) P(I)=(8*(Q3(I)-Q1(I))-(Q4(I)-Q0(I)))/(12*H) F1(I)=F2(I) F2(I)=F3(I) Z1(I)=Z2(I) Z2(I)=Z3(I) Q0(I)=Q1(I) Q1(I)=Q2(I) Q2(I)=Q3(I) Q3(I)=Q4(I) END DO XOLD=X X=X+H NR=NR+1 IF (IOUT.GE.1) CALL SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IF (NR.LT.NSTEP.AND.IRTRN.GE.0) GOTO 44 END IF C IF (KK.EQ.6) THEN DO I=1,N U(I)=A1*(Z0(I)+Z4(I))+A2*(Z1(I)+Z3(I))+A3*Z2(I) END DO 66 CONTINUE CALL FCN(N,X,Q5,F5,RPAR,IPAR) DO I=1,N U(I)=U(I)+(B1*(F1(I)+F5(I))+B2*(F2(I)+F4(I))+B3*F3(I)) SUM=A1*Z1(I)+A2*(Z2(I)+Z4(I))+A3*Z3(I) Z5(I)=U(I)-SUM Q6(I)=Q5(I)+H*Z5(I) END DO DO I=1,N Q(I)=Q3(I) P(I)=(45*(Q4(I)-Q2(I))-9*(Q5(I)-Q1(I))+(Q6(I)-Q0(I)))/(60*H) F1(I)=F2(I) F2(I)=F3(I) F3(I)=F4(I) F4(I)=F5(I) Z1(I)=Z2(I) Z2(I)=Z3(I) Z3(I)=Z4(I) Z4(I)=Z5(I) Q0(I)=Q1(I) Q1(I)=Q2(I) Q2(I)=Q3(I) Q3(I)=Q4(I) Q4(I)=Q5(I) Q5(I)=Q6(I) END DO XOLD=X X=X+H NR=NR+1 IF (IOUT.GE.1) CALL SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IF (NR.LT.NSTEP.AND.IRTRN.GE.0) GOTO 66 END IF C IF (KK.EQ.8) THEN DO I=1,N U(I)=A1*(Z0(I)+Z6(I))+A2*(Z1(I)+Z5(I)) & +A3*(Z2(I)+Z4(I))+A4*Z3(I) END DO 88 CONTINUE CALL FCN(N,X,Q7,F7,RPAR,IPAR) DO I=1,N U(I)=U(I)+(B1*(F1(I)+F7(I))+B2*(F2(I)+F6(I)) & +B3*(F3(I)+F5(I))+B4*F4(I)) SUM=A1*Z1(I)+A2*(Z2(I)+Z6(I))+A3*(Z3(I)+Z5(I))+A4*Z4(I) Z7(I)=U(I)-SUM Q8(I)=Q7(I)+H*Z7(I) END DO DO I=1,N Q(I)=Q4(I) P(I)=(672*(Q5(I)-Q3(I))-168*(Q6(I)-Q2(I)) & +32*(Q7(I)-Q1(I))-3*(Q8(I)-Q0(I)))/(840*H) F1(I)=F2(I) F2(I)=F3(I) F3(I)=F4(I) F4(I)=F5(I) F5(I)=F6(I) F6(I)=F7(I) Z1(I)=Z2(I) Z2(I)=Z3(I) Z3(I)=Z4(I) Z4(I)=Z5(I) Z5(I)=Z6(I) Z6(I)=Z7(I) Q0(I)=Q1(I) Q1(I)=Q2(I) Q2(I)=Q3(I) Q3(I)=Q4(I) Q4(I)=Q5(I) Q5(I)=Q6(I) Q6(I)=Q7(I) Q7(I)=Q8(I) END DO XOLD=X X=X+H NR=NR+1 IF (IOUT.GE.1) CALL SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IF (NR.LT.NSTEP.AND.IRTRN.GE.0) GOTO 88 END IF RETURN END C SUBROUTINE PREPST(N,FCN,H,Q,QB,QE,Z,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),QB(N),QE(N),Z(N),F(N) DIMENSION IPAR(*),RPAR(*) DO I=1,N QE(I)=Q(I) Z(I)=(Q(I)-QB(I))/H END DO CALL FCN(N,X,QE,F,RPAR,IPAR) RETURN END C SUBROUTINE SYCOE(METH,A1,A2,A3,A4,B1,B2,B3,B4,H,HM) IMPLICIT REAL*8 (A-H,O-Z) IF (METH.EQ.201) RETURN IF (METH.EQ.401) THEN C --- METHOD ORDER 4 A1=1.0D0 A2=1.0D0 HM=H/4.0D0 B1=5*HM B2=2*HM RETURN END IF IF (METH.EQ.601) THEN C --- METHOD ORDER 6 A1=1.0D0 A2=1.0D0 A3=1.0D0 HM=H/48.0D0 B1=67*HM B2=-8*HM B3=122*HM RETURN END IF IF (METH.EQ.801) THEN C --- METHOD SY8 (QUINLAN & TREMAINE 1990) A1=1.0D0 A2=0.0D0 A3=1.0D0 A4=1.0D0 HM=H/12096.0D0 B1=17671*HM B2=-23622*HM B3=61449*HM B4=-50516*HM RETURN END IF IF (METH.EQ.802) THEN C --- METHOD SY8B A1=1.0D0 A2=2.0D0 A3=3.0D0 A4=3.5D0 HM=H/120960.0D0 B1=192481*HM B2=6582*HM B3=816783*HM B4=-156812*HM RETURN END IF IF (METH.EQ.803) THEN C --- METHOD SY8C A1=1.0D0 A2=1.0D0 A3=1.0D0 A4=1.0D0 HM=H/8640.0D0 B1=13207*HM B2=-8934*HM B3=42873*HM B4=-33812*HM RETURN END IF METH=-METH RETURN END gnicodes/problem.f100644 12122 12120 13522 7467146630 13252 0ustar hairermathC C IPROB = 1 : KEPLER PROBLEM, ECCENTRICITY IN RPAR(1) C IPROB = 2 : HARMONIC OSCILLATOR C IPROB = 3 : PENDULUM C IPROB = 4 : OUTER SOLAR SYSTEM C SUBROUTINE PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(NDIM),P(NDIM),QEX(NDIM),PEX(NDIM) DIMENSION IPAR(*),RPAR(*) IPROB=IPAR(11) IF (IPROB.EQ.1) THEN N=2 PI=4*ATAN(1.D0) ECCENT=RPAR(11) X=0.D0 XEND=2*PI Q(1)=1-ECCENT Q(2)=0.d0 P(1)=0.d0 P(2)=SQRT((1+ECCENT)/(1-ECCENT)) QEX(1)=1-ECCENT QEX(2)=0.d0 PEX(1)=0.d0 PEX(2)=SQRT((1+ECCENT)/(1-ECCENT)) END IF IF (IPROB.EQ.2) THEN N=1 PI=4*ATAN(1.D0) X=0.D0 XEND=2*PI Q(1)=0.0D0 P(1)=1.0D0 QEX(1)=0.0D0 PEX(1)=1.0D0 END IF IF (IPROB.EQ.3) THEN N=1 PI=4*ATAN(1.D0) X=0.D0 XEND=2*PI Q(1)=0.0D0 P(1)=1.0D0 QEX(1)=-0.443944662290259D0 PEX(1)= 0.897846803312944D0 END IF IF (IPROB.EQ.4) THEN N=18 X=0.D0 XEND=500000.D0 Q(1)=-3.5023653D0 Q(2)=-3.8169847D0 Q(3)=-1.5507963D0 Q(4)=9.0755314D0 Q(5)=-3.0458353D0 Q(6)=-1.6483708D0 Q(7)=8.3101420D0 Q(8)=-16.2901086D0 Q(9)=-7.2521278D0 Q(10)=11.4707666D0 Q(11)=-25.7294829D0 Q(12)=-10.8169456D0 Q(13)=-15.5387357D0 Q(14)=-25.2225594D0 Q(15)=-3.1902382D0 Q(16)=0.0D0 Q(17)=0.0D0 Q(18)=0.0D0 P(1)=0.00565429D0 P(2)=-0.00412490D0 P(3)=-0.00190589D0 P(4)=0.00168318D0 P(5)=0.00483525D0 P(6)=0.00192462D0 P(7)=0.00354178D0 P(8)=0.00137102D0 P(9)=0.00055029D0 P(10)=0.00288930D0 P(11)=0.00114527D0 P(12)=0.00039677D0 P(13)=0.00276725D0 P(14)=-0.00170702D0 P(15)=-0.00136504D0 P(16)=0.0D0 P(17)=0.0D0 P(18)=0.0D0 QEX( 1)= 0.7766584086800482D+01 QEX( 2)= 0.2531065754551048D+00 QEX( 3)= -0.9410571402013185D-01 QEX( 4)= -0.5564967162844037D+01 QEX( 5)= 0.1674849740822012D+01 QEX( 6)= 0.9767232069533176D+00 QEX( 7)= 0.1963899572895227D+02 QEX( 8)= 0.8958504552286460D+01 QEX( 9)= 0.3611839157057347D+01 QEX(10)= 0.2493570870305177D+02 QEX(11)= 0.1769518676153705D+02 QEX(12)= 0.6583785164549242D+01 QEX(13)= 0.3178592511375764D+02 QEX(14)= 0.3863618958160644D+02 QEX(15)= 0.3192794169732889D+01 QEX(16)= 0.3084118473380683D+01 QEX(17)= -0.1227726356581642D+01 QEX(18)= -0.6162537634647217D+00 PEX( 1)= -0.2495503201917009D-02 PEX( 2)= 0.6896467194473328D-02 PEX( 3)= 0.3007950247474123D-02 PEX( 4)= -0.2255335935351989D-02 PEX( 5)= -0.4905913854771086D-02 PEX( 6)= -0.1938473641716708D-02 PEX( 7)= -0.2186170231167942D-02 PEX( 8)= 0.2817177012110666D-02 PEX( 9)= 0.1262882639181183D-02 PEX(10)= -0.2148728705895163D-02 PEX(11)= 0.2128650077635786D-02 PEX(12)= 0.9248501411662923D-03 PEX(13)= -0.1675173186229401D-02 PEX(14)= 0.1011833320388655D-02 PEX(15)= 0.8231800038576520D-03 PEX(16)= 0.9417379703028725D-05 PEX(17)= -0.7855256238249194D-05 PEX(18)= -0.3646926313230521D-05 END IF RETURN END c SUBROUTINE EQUA(N,X,Q,F,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION D(6,6),K,M(6) DIMENSION Q(N),F(N) DIMENSION IPAR(*),RPAR(*) IPROB=IPAR(11) IPAR(12)=IPAR(12)+1 IF (IPROB.EQ.1) THEN RAD=Q(1)*Q(1)+Q(2)*Q(2) RAD=RAD*SQRT(RAD) F(1)=-Q(1)/RAD F(2)=-Q(2)/RAD END IF IF (IPROB.EQ.2) THEN F(1)=-Q(1) END IF IF (IPROB.EQ.3) THEN F(1)=-SIN(Q(1)) END IF IF (IPROB.EQ.4) THEN K=2.95912208286D-4 M(1)=0.000954786104043D0 M(2)=0.000285583733151D0 M(3)=0.0000437273164546D0 M(4)=0.0000517759138449D0 M(5)=1.0D0/1.3D8 M(6)=1.00000597682D0 DO I=1,5 I1=3*(I-1)+1 DO J=I+1,6 J1=3*(J-1)+1 D(I,J)=(SQRT((Q(I1)-Q(J1))**2+(Q(I1+1)-Q(J1+1))**2+ * (Q(I1+2)-Q(J1+2))**2))**3 D(J,I)=D(I,J) END DO END DO DO I=1,6 I1=3*(I-1)+1 F(I1)=0.0D0 F(I1+1)=0.0D0 F(I1+2)=0.0D0 DO J=1,6 IF (J.NE.I) THEN J1=3*(J-1)+1 F(I1)=F(I1)+M(J)*(Q(J1)-Q(I1))/D(I,J) F(I1+1)=F(I1+1)+M(J)*(Q(J1+1)-Q(I1+1))/D(I,J) F(I1+2)=F(I1+2)+M(J)*(Q(J1+2)-Q(I1+2))/D(I,J) END IF END DO F(I1)=K*F(I1) F(I1+1)=K*F(I1+1) F(I1+2)=K*F(I1+2) END DO END IF RETURN END c SUBROUTINE HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION D(6,6),K,M(6),Q(N),P(N),Y(100) DIMENSION IPAR(*),RPAR(*) IPROB=IPAR(11) IF (IPROB.EQ.1) THEN HAMIL=0.5d0*(P(1)**2+P(2)**2)-1/SQRT(Q(1)**2+Q(2)**2) END IF IF (IPROB.EQ.2) THEN HAMIL=0.5d0*(P(1)**2+Q(1)**2) END IF IF (IPROB.EQ.3) THEN HAMIL=0.5d0*P(1)**2-COS(Q(1)) END IF IF (IPROB.EQ.4) THEN DO I=1,N Y(I)=Q(I) Y(I+N)=P(I) END DO K=2.95912208286D-4 M(1)=0.000954786104043D0 M(2)=0.000285583733151D0 M(3)=0.0000437273164546D0 M(4)=0.0000517759138449D0 M(5)=1.0D0/1.3D8 M(6)=1.00000597682D0 DO I=1,5 I1=3*(I-1)+1 DO J=I+1,6 J1=3*(J-1)+1 D(I,J)=SQRT((Y(I1)-Y(J1))**2+(Y(I1+1)-Y(J1+1))**2+ * (Y(I1+2)-Y(J1+2))**2) D(J,I)=D(I,J) END DO END DO HAMIL=0.0D0 DO I=1,6 I1=N+3*(I-1)+1 HAMIL=HAMIL+M(I)*(Y(I1)**2+Y(I1+1)**2+Y(I1+2)**2) END DO HAMIL=HAMIL/2.0D0 POT=0.0D0 DO I=2,6 DO J=1,I-1 POT=POT+M(I)*M(J)/D(I,J) END DO END DO HAMIL=HAMIL-K*POT END IF RETURN END gnicodes/dr_irk2.f100644 12122 12120 3121 7535337725 13123 0ustar hairermathccfeh dr_irk2 include 'gni_irk2.f' include 'problem.f' C --- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX C --- CHOOSE THE PROBLEM C IPROB = 1 : KEPLER PROBLEM, ECCENTRICITY IN RPAR(1) C IPROB = 2 : HARMONIC OSCILLATOR C IPROB = 3 : PENDULUM C IPROB = 4 : OUTER SOLAR SYSTEM IPROB=4 IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) C --- CHOOSE THE THE METHOD C --- GAUSS METHOD OF ORDER 2*METH METH=6 WRITE (6,*) ' METHOD, PROBLEM ',METH,IPROB H=1.0D2 NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=1 CALL CPU_TIME(TIME0) CALL GNI_IRK2(N,EQUA,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR,RPAR(13) 90 FORMAT (1X,I8,F7.2,1X,2E24.16) STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF RETURN END gnicodes/dr_comp.f100644 12122 12120 3203 7624126561 13204 0ustar hairermathccfeh dr_comp include 'gni_comp.f' include 'problem.f' C --- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX,STVERL C --- CHOOSE THE PROBLEM C IPROB = 1 : KEPLER PROBLEM, ECCENTRICITY IN RPAR(1) C IPROB = 2 : HARMONIC OSCILLATOR C IPROB = 3 : PENDULUM C IPROB = 4 : OUTER SOLAR SYSTEM IPROB=1 IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) C --- CHOOSE THE THE COMPOSITION METHOD C --- ORDER-STAGES : 21,43,45,67,69,815,817,1033,1034 METHC=817 WRITE (6,*) ' METHOD, PROBLEM ',METHC,IPROB H=0.1D0 NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=1 CALL CPU_TIME(TIME0) CALL GNI_COMP(STVERL,N,EQUA,NSTEP,X,P,Q,XEND, & METHC,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR,RPAR(13) 90 FORMAT (1X,I8,F7.2,1X,2E24.16) STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF RETURN END gnicodes/dr_lmm2.f100644 12122 12120 3322 7561742672 13126 0ustar hairermathccfeh dr_lmm2 include 'gni_lmm2.f' include 'gni_irk2.f' include 'problem.f' C --- IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=100) DIMENSION P(NDIM),Q(NDIM),PEX(NDIM),QEX(NDIM),IPAR(20),RPAR(20) REAL TIME0,TIME1 EXTERNAL EQUA,SOLFIX C --- CHOOSE THE PROBLEM C IPROB = 1 : KEPLER PROBLEM, ECCENTRICITY IN RPAR(1) C IPROB = 2 : HARMONIC OSCILLATOR C IPROB = 3 : PENDULUM C IPROB = 4 : OUTER SOLAR SYSTEM IPROB=3 IPAR(11)=IPROB IF (IPROB.EQ.1) RPAR(11)=0.6D0 CALL PDATA(X,XEND,NDIM,N,Q,P,QEX,PEX,RPAR,IPAR) C --- CHOOSE THE METHOD C --- SYMMETRIC LINEAR MULTISTEP METHOD (ORDER 2,4,6,8) C --- FOR THE MOMENT : 201,401,601,801,802,803. METH=801 WRITE (6,*) ' METHOD, PROBLEM ',METH,IPROB H=1.0D-1 NSTEP=(XEND-X)/H C --- CALL OF THE METHOD DO I=1,10 RPAR(I)=0.0D0 IPAR(I)=0 END DO IPAR(12)=0 IOUT=1 CALL CPU_TIME(TIME0) CALL GNI_LMM2(N,EQUA,NSTEP,X,P,Q,XEND, & METH,SOLFIX,IOUT,RPAR,IPAR) CALL CPU_TIME(TIME1) C --- STATISTICS NFCN=IPAR(12) ERR=0.0D0 DO I=1,N ERR=ERR+(Q(I)-QEX(I))**2+(P(I)-PEX(I))**2 END DO ERR=SQRT(ERR) WRITE (6,90) NFCN,TIME1-TIME0,ERR,RPAR(13) 90 FORMAT (1X,I8,F7.2,1X,2E24.16) STOP END C SUBROUTINE SOLFIX (NR,XOLD,X,P,Q,N,IRTRN,RPAR,IPAR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N),P(N) DIMENSION IPAR(20),RPAR(20) CALL HAMILTON (N,X,Q,P,HAMIL,RPAR,IPAR) IF (NR.EQ.0) THEN RPAR(12)=HAMIL RPAR(13)=0.0D0 ELSE RPAR(13)=MAX(RPAR(13),ABS(RPAR(12)-HAMIL)) END IF C write (6,*) nr,x,q(1),p(1) RETURN END