(original) (raw)

SUBROUTINE ODEX2(N,FCN,X,Y,YP,XEND,H, & RTOL,ATOL,ITOL, & SOLOUT,IOUT, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C ---------------------------------------------------------- C NUMERICAL SOLUTION OF A SYSTEM OF SECOND 0RDER C ORDINARY DIFFERENTIAL EQUATIONS Y''=F(X,Y). C THIS IS AN EXTRAPOLATION-ALGORITHM, BASED ON C THE STOERMER RULE (WITH STEPSIZE CONTROL C ORDER SELECTION AND DENSE OUTPUT). C C AUTHORS: E. HAIRER AND G. WANNER C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES C CH-1211 GENEVE 24, SWITZERLAND C E-MAIL: Ernst.Hairer@math.unige.ch C Gerhard.Wanner@math.unige.ch C C THIS CODE IS DESCRIBED IN SECTION II.14 OF THE BOOK: C E. HAIRER, S.P. NORSETT AND G. WANNER, SOLVING ORDINARY C DIFFERENTIAL EQUATIONS I. NONSTIFF PROBLEMS. 2ND EDITION. C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, C SPRINGER-VERLAG (1993) C C VERSION SEPTEMBER 30, 1995 C SMALL CORRECTIONS ON JUNE 11, 1999 C C INPUT PARAMETERS C ---------------- C N DIMENSION OF THE SYSTEM C C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE C VALUE OF F(X,Y): C SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) C DOUBLE PRECISION X,Y(N),F(N) C F(1)=... ETC. C C X INITIAL X-VALUE C C Y(N) INITIAL VALUES FOR Y C C YP(N) INITIAL VALUES FOR Y' C C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) C C H INITIAL STEP SIZE GUESS; C USUALLY 1.D-1 OR 1.D-3, IS GOOD. THIS CHOICE IS NOT C VERY IMPORTANT, THE CODE QUICKLY ADAPTS ITS STEPSIZE. C WHEN YOU ARE NOT SURE, THEN STUDY THE CHOSEN VALUES C FOR A FEW STEPS IN SUBROUTINE "SOLOUT". C (IF H=0.D0, THE CODE PUTS H=1.D-4). C C RTOL,ATOL RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY CAN C BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH 2*N. C C ITOL SWITCH FOR RTOL AND ATOL: C ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF C Y(I) BELOW RTOL*ABS(Y(I))+ATOL C YP(I) BELOW RTOL*ABS(YP(I))+ATOL C ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. C THE CODE KEEPS THE LOCAL ERROR OF C Y(I) BELOW RTOL(I)*ABS(Y(I))+ATOL(I). C YP(I) BELOW RTOL(I+N)*ABS(YP(I))+ATOL(I+N). C C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE C NUMERICAL SOLUTION DURING INTEGRATION. C IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. C IT MUST HAVE THE FORM C SUBROUTINE SOLOUT (NR,XOLD,X,Y,YP,N,CON,NCON,ICOMP,ND, C RPAR,IPAR,IRTRN) C DIMENSION Y(N),YP(N),CON(NCON),ICOMP(ND) C .... C SOLOUT FURNISHES THE SOLUTIONS "Y, YP" AT THE NR-TH C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS C THE FIRST GRID-POINT). C "XOLD" IS THE PRECEEDING GRID-POINT. C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN C IS SET <0, ODEX2 WILL RETURN TO THE CALLING PROGRAM. C C ----- CONTINUOUS OUTPUT (IF IOUT=2): ----- C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH C THE DOUBLE PRECISION FUNCTION C >>> CONTX2(I,S,CON,NCON,ICOMP,ND) <<< C WHICH PROVIDES AN APPROXIMATION TO THE I-TH C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE C S SHOULD LIE IN THE INTERVAL [XOLD,X]. C C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLOUT: C IOUT=0: SUBROUTINE IS NEVER CALLED C IOUT=1: SUBROUTINE IS USED FOR OUTPUT C IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT C C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". C SERVES AS WORKING SPACE FOR ALL VECTORS. C "LWORK" MUST BE AT LEAST C N*(2*KM+6)+5*KM+20+(KM*(2*KM+5)+6)*NRDENS C WHERE NRDENS=IWORK(8) (SEE BELOW) AND C KM=9 IF IWORK(2)=0 C KM=IWORK(2) IF IWORK(2).GT.0 C WORK(1),...,WORK(20) SERVE AS PARAMETERS C FOR THE CODE. FOR STANDARD USE, SET THESE C PARAMETERS TO ZERO BEFORE CALLING. C C LWORK DECLARED LENGTH OF ARRAY "WORK". C C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". C "LIWORK" MUST BE AT LEAST C KM+20+NRDENS C IWORK(1),...,IWORK(20) SERVE AS PARAMETERS C FOR THE CODE. FOR STANDARD USE, SET THESE C PARAMETERS TO ZERO BEFORE CALLING. C C LIWORK DECLARED LENGTH OF ARRAY "IWORK". C C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH C CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING C PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. C C----------------------------------------------------------------------- C C SOPHISTICATED SETTING OF PARAMETERS C ----------------------------------- C SEVERAL PARAMETERS (WORK(1),...,IWORK(1),...) ALLOW C TO ADAPT THE CODE TO THE PROBLEM AND TO THE NEEDS OF C THE USER. FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES. C C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 2.3D-16. C C WORK(2) MAXIMAL STEP SIZE, DEFAULT XEND-X. C C WORK(3) STEP SIZE IS REDUCED BY FACTOR WORK(3), IF DURING THE C COMPUTATION OF THE EXTRAPOLATION TABLEAU DIVERGENCE C IS OBSERVED; DEFAULT 0.5. C C WORK(4), WORK(5) PARAMETERS FOR STEP SIZE SELECTION C THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS C CHOSEN SUBJECT TO THE RESTRICTION C FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN C WHERE FACMIN=WORK(4)**(1/(2*J-1)) C DEFAULT VALUES: WORK(4)=0.02D0, WORK(5)=4.D0 C C WORK(6), WORK(7) PARAMETERS FOR THE ORDER SELECTION C STEP SIZE IS DECREASED IF W(K-1) <= W(K)*WORK(6) C STEP SIZE IS INCREASED IF W(K) <= W(K-1)*WORK(7) C DEFAULT VALUES: WORK(6)=0.8D0, WORK(7)=0.9D0 C C WORK(8), WORK(9) SAFETY FACTORS FOR STEP CONTROL ALGORITHM C HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1)) C DEFAULT VALUES: WORK(8)=0.65D0, C WORK(9)=0.94D0 IF "HOPE FOR CONVERGENCE" C C IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. C THE DEFAULT VALUE (FOR IWORK(1)=0) IS 10000. C C IWORK(2) THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION C TABLE. THE DEFAULT VALUE (FOR IWORK(2)=0) IS 9. C IF IWORK(2).NE.0 THEN IWORK(2) SHOULD BE .GE.3. C C IWORK(3) SWITCH FOR THE STEP SIZE SEQUENCE (EVEN NUMBERS ONLY) C IF IWORK(3).EQ.1 THEN 2,4,6,8,10,12,14,16,... C IF IWORK(3).EQ.2 THEN 2,4,8,12,16,20,24,28,... C IF IWORK(3).EQ.3 THEN 2,4,6,8,12,16,24,32,... C IF IWORK(3).EQ.4 THEN 2,6,10,14,18,22,26,30,... C THE DEFAULT VALUE IS IWORK(3)=1 IF IOUT.LE.1; C THE DEFAULT VALUE IS IWORK(3)=4 IF IOUT.GE.2. C C IWORK(6) IF IWORK(6)=0 ERROR ESTIMATOR IN THE DENSE C OUTPUT FORMULA IS ACTIVATED. IT CAN BE SUPPRESSED C BY PUTTING IWORK(6)=1. C DEFAULT IWORK(6)=0 (IF IOUT.GE.2). C C IWORK(7) DETERMINES THE DEGREE OF INTERPOLATION FORMULA C MU = 2 * KAPPA - IWORK(7) + 1 C IWORK(7) SHOULD LIE BETWEEN 1 AND 8 C DEFAULT IWORK(7)=6 (IF IWORK(7)=0). C C IWORK(8) = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT C IS REQUIRED C C IWORK(21),...,IWORK(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH C DENSE OUTPUT IS REQUIRED. THEY NEED NOT BE SET IF C NRDENS = N . C C----------------------------------------------------------------------C C OUTPUT PARAMETERS C ----------------- C X X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED C (AFTER SUCCESSFUL RETURN X=XEND). C C Y(N) NUMERICAL SOLUTION AT X C C YP(N) NUMERICAL DERIVATIVE OF SOLUTION AT X C C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP C C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: C IDID=1 COMPUTATION SUCCESSFUL, C IDID=-1 COMPUTATION UNSUCCESSFUL. C C IWORK(17) NFCN NUMBER OF FUNCTION EVALUATIONS C IWORK(18) NSTEP NUMBER OF COMPUTED STEPS C IWORK(19) NACCPT NUMBER OF ACCEPTED STEPS C IWORK(20) NREJCT NUMBER OF REJECTED STEPS (DUE TO ERROR TEST), C (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED) C----------------------------------------------------------------------- C *** *** *** *** *** *** *** *** *** *** *** *** *** C DECLARATIONS C *** *** *** *** *** *** *** *** *** *** *** *** *** IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),YP(N),WORK(LWORK),IWORK(LIWORK) DIMENSION ATOL(*),RTOL(*),RPAR(*),IPAR(*) LOGICAL ARRET EXTERNAL FCN,SOLOUT C *** *** *** *** *** *** *** C SETTING THE PARAMETERS C *** *** *** *** *** *** *** NFCN=0 NSTEP=0 NACCPT=0 NREJCT=0 ARRET=.FALSE. C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- IF(IWORK(1).EQ.0)THEN NMAX=10000 ELSE NMAX=IWORK(1) IF(NMAX.LE.0)THEN WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1) ARRET=.TRUE. END IF END IF C -------- KM MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION IF(IWORK(2).EQ.0)THEN KM=9 ELSE KM=IWORK(2) IF(KM.LE.2)THEN WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2) ARRET=.TRUE. END IF END IF C -------- NSEQU CHOICE OF STEP SIZE SEQUENCE NSEQU=IWORK(3) IF(IWORK(3).EQ.0.AND.IOUT.LE.1) NSEQU=1 IF(IWORK(3).EQ.0.AND.IOUT.GE.2) NSEQU=4 IF(NSEQU.LE.0.OR.NSEQU.GE.5)THEN WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3) ARRET=.TRUE. END IF IF (NSEQU.LE.3.AND.IOUT.GE.2) THEN WRITE(6,*)' IWORK(3) NOT COMPATIBLE WITH IOUT' ARRET=.TRUE. END IF C -------- IDERR PARAMETER FOR ERROR ESTIMATION IN DENSE OUTPUT IF(IWORK(6).EQ.0)THEN IF(IOUT.LE.1) IDERR=1 IF(IOUT.GE.2) IDERR=0 ELSE IDERR=IWORK(6) IF(IOUT.LE.1)THEN WRITE(6,*)' ERROR ESTIMATION IN DENSE OUTPUT', & ' NOT POSSIBLE, WRONG IWORK(6)=',IWORK(6) ARRET=.TRUE. END IF END IF C -------- MUDIF IF(IWORK(7).EQ.0)THEN MUDIF=6 ELSE MUDIF=IWORK(7) IF(MUDIF.LE.0.OR.MUDIF.GE.9)THEN WRITE(6,*)' WRONG INPUT IWORK(7)=',IWORK(7) ARRET=.TRUE. END IF END IF C -------- NRDENS NUMBER OF DENSE OUTPUT COMPONENTS NRDENS=IWORK(8) IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN WRITE(6,*)' CURIOUS INPUT IWORK(8)=',IWORK(8) ARRET=.TRUE. END IF IF (NRDENS.EQ.N) THEN DO 17 I=1,NRDENS 17 IWORK(20+I)=I END IF C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 IF(WORK(1).EQ.0.D0)THEN UROUND=2.3D-16 ELSE UROUND=WORK(1) IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN WRITE(6,*)' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:' & ,WORK(1) ARRET=.TRUE. END IF END IF C -------- MAXIMAL STEP SIZE IF(WORK(2).EQ.0.D0)THEN HMAX=XEND-X ELSE HMAX=ABS(WORK(2)) END IF C -------- STEP SIZE REDUCTION FACTOR IF(WORK(3).EQ.0.D0)THEN SAFE3=0.5D0 ELSE SAFE3=WORK(3) IF(SAFE3.LE.UROUND.OR.SAFE3.GE.1.D0)THEN WRITE(6,*)' CURIOUS INPUT WORK(3)=',WORK(3) ARRET=.TRUE. END IF END IF C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION IF(WORK(4).EQ.0.D0)THEN FAC1=0.02D0 ELSE FAC1=WORK(4) END IF IF(WORK(5).EQ.0.D0)THEN FAC2=4.0D0 ELSE FAC2=WORK(5) END IF C ------- FAC3, FAC4 PARAMETERS FOR THE ORDER SELECTION IF(WORK(6).EQ.0.D0)THEN FAC3=0.8D0 ELSE FAC3=WORK(6) END IF IF(WORK(7).EQ.0.D0)THEN FAC4=0.9D0 ELSE FAC4=WORK(7) END IF C ------- SAFE1, SAFE2 SAFETY FACTORS FOR STEP SIZE PREDICTION IF(WORK(8).EQ.0.D0)THEN SAFE1=0.65D0 ELSE SAFE1=WORK(8) END IF IF(WORK(9).EQ.0.D0)THEN SAFE2=0.94D0 ELSE SAFE2=WORK(9) END IF C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- LFSAFE=2*KM*KM+KM IEDY=21 IEYH1=IEDY+N IEYH2=IEYH1+N IEDZ=IEYH2+N IESCAL=IEDZ+N IESCP=IESCAL+N IET=IESCP+N IETP=IET+KM*N IEFS=IETP+KM*N IEYS=IEFS+LFSAFE*NRDENS IEYSP=IEYS+KM*NRDENS IEHH=IEYSP+KM*NRDENS IEW=IEHH+KM IEA=IEW+KM IEFAC=IEA+KM C ------ TOTAL STORAGE REQUIREMENT ----------- IECO=IEFAC+2*KM ISTORE=IECO+(2*KM+6)*NRDENS-1 IF(ISTORE.GT.LWORK)THEN WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE ARRET=.TRUE. END IF C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- ICOM=21 IENJ=ICOM+NRDENS C --------- TOTAL REQUIREMENT --------------- ISTORE=IENJ+KM-1 IF(ISTORE.GT.LIWORK)THEN WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE ARRET=.TRUE. END IF C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 IF (ARRET) THEN IDID=-1 RETURN END IF C -------- CALL TO CORE INTEGRATOR ------------ NRD=MAX(1,NRDENS) NCOM=MAX(1,(2*KM+6)*NRDENS) CALL ODX2CO(N,FCN,X,Y,YP,XEND,HMAX,H,RTOL,ATOL,ITOL,KM, & SOLOUT,IOUT,IDID,NMAX,UROUND,WORK(IEDY),WORK(IEYH1), & WORK(IEYH2),WORK(IEDZ),WORK(IESCAL),WORK(IESCP),WORK(IEFS), & WORK(IEYS),WORK(IEYSP),WORK(IET),WORK(IETP),WORK(IEHH), & WORK(IEW),WORK(IEA),WORK(IECO),NCOM,IWORK(ICOM),IWORK(IENJ), & NSEQU,LFSAFE,SAFE1,SAFE2,SAFE3,FAC1,FAC2,FAC3,FAC4,IDERR, & WORK(IEFAC),MUDIF,NRD,RPAR,IPAR,NFCN,NSTEP,NACCPT,NREJCT) IWORK(17)=NFCN IWORK(18)=NSTEP IWORK(19)=NACCPT IWORK(20)=NREJCT C ----------- RETURN ----------- RETURN END C C C C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- C SUBROUTINE ODX2CO(N,FCN,X,Y,YP,XEND,HMAX,H,RTOL,ATOL,ITOL,KM, & SOLOUT,IOUT,IDID,NMAX,UROUND,DY,YH1,YH2,DZ,SCAL,SCALP,FSAFE, & YSAFE,YSAFEP,T,TP,HH,W,A,DENS,NCOM,ICOMP,NJ,NSEQU,LFSAFE, & SAFE1,SAFE2,SAFE3,FAC1,FAC2,FAC3,FAC4,IDERR,ERRFAC,MUDIF, & NRD,RPAR,IPAR,NFCN,NSTEP,NACCPT,NREJCT) C ---------------------------------------------------------- C CORE INTEGRATOR FOR ODEX2 C PARAMETERS SAME AS IN ODEX2 WITH WORKSPACE ADDED C ---------------------------------------------------------- C DECLARATIONS C ---------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL REJECT,LAST,ATOV EXTERNAL FCN DIMENSION Y(N),YP(N),DY(N),YH1(N),YH2(N),DZ(N),SCAL(N),SCALP(N) DIMENSION T(KM,N),TP(KM,N),NJ(KM),HH(KM),W(KM),A(KM) DIMENSION RTOL(*),ATOL(*),RPAR(*),IPAR(*) DIMENSION FSAFE(LFSAFE,NRD),YSAFE(KM,NRD),YSAFEP(KM,NRD) DIMENSION ERRFAC(2*KM),DENS(NCOM),ICOMP(NRD) COMMON /CONOD2/XOLDD,HHH,KMIT C --- DEFINE THE STEP SIZE SEQUENCE IF (NSEQU.EQ.1) THEN DO 1 I=1,KM 1 NJ(I)=2*I END IF IF (NSEQU.EQ.2) THEN NJ(1)=2 DO 2 I=2,KM 2 NJ(I)=4*I-4 END IF IF (NSEQU.EQ.3) THEN NJ(1)=2 NJ(2)=4 NJ(3)=6 DO 11 I=4,KM 11 NJ(I)=2*NJ(I-2) END IF IF (NSEQU.EQ.4) THEN DO 3 I=1,KM 3 NJ(I)=4*I-2 END IF C --- DEFINE THE A(I) FOR ORDER SELECTION A(1)=1.D0+NJ(1)/2 DO 4 I=2,KM 4 A(I)=A(I-1)+NJ(I)/2 C --- INITIAL SCALING IF (ITOL.EQ.0) THEN DO 6 I=1,N SCAL(I)=ATOL(1)+RTOL(1)*ABS(Y(I)) 6 SCALP(I)=ATOL(1)+RTOL(1)*ABS(YP(I)) ELSE DO 8 I=1,N SCAL(I)=ATOL(I)+RTOL(I)*ABS(Y(I)) 8 SCALP(I)=ATOL(I+N)+RTOL(I+N)*ABS(YP(I)) END IF C --- INITIAL PREPARATIONS POSNEG=SIGN(1.D0,XEND-X) K=MAX(2,MIN(KM-1,INT(-LOG10(RTOL(1))*0.6D0+1.5D0))) HMAX=ABS(HMAX) H=MAX(ABS(H),1.D-4) H=POSNEG*MIN(H,HMAX,ABS(XEND-X)) IF (IOUT.GE.1) THEN IF (IOUT.GE.2) THEN DO 9 MU=1,KM*2 ERRX=SQRT(MU/(MU+6.D0))*0.5D0 PROD=(1.5D0/(MU+6.D0))**3 DO 7 J=1,MU 7 PROD=PROD*ERRX/J 9 ERRFAC(MU)=PROD END IF NNRD=NRD IRTRN=0 XOLD=X CALL SOLOUT (NACCPT+1,XOLD,X,Y,YP,N,DENS,NCOM,ICOMP,NRD, & RPAR,IPAR,IRTRN) IF (IRTRN.LT.0) GOTO 120 END IF ERR=0.D0 ERROLD=1.D10 HOPTDE=POSNEG*HMAX W(1)=0.D0 REJECT=.FALSE. LAST=.FALSE. 10 ATOV=.FALSE. C --- IS XEND REACHED IN THE NEXT STEP? IF (0.1D0*ABS(XEND-X).LE.ABS(X)*UROUND)GOTO 110 H=POSNEG*MIN(ABS(H),ABS(XEND-X),HMAX,ABS(HOPTDE)) IF ((X+1.01D0*H-XEND)*POSNEG.GT.0.D0) THEN H=XEND-X LAST=.TRUE. END IF IF (NSTEP.EQ.0.OR.IOUT.NE.2) CALL FCN(N,X,Y,DZ,RPAR,IPAR) NFCN=NFCN+1 C --- THE FIRST AND LAST STEP IF (NSTEP.EQ.0.OR.LAST) THEN IPT=0 NSTEP=NSTEP+1 DO 20 J=1,K KC=J CALL STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W, 1 ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV, 2 SAFE3,REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE, 3 LFSAFE,IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN) IF (ATOV) GO TO 10 20 IF (J.GT.1.AND.ERR.LE.1.D0) GO TO 60 GO TO 55 END IF C --- BASIC INTEGRATION STEP 30 CONTINUE IPT=0 NSTEP=NSTEP+1 IF (NSTEP.GE.NMAX) GO TO 120 KC=K-1 DO 40 J=1,KC CALL STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W, 1 ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3, 2 REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE, 3 IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN) IF (ATOV) GO TO 10 40 CONTINUE C --- CONVERGENCE MONITOR IF (K.EQ.2.OR.REJECT) GO TO 50 IF (ERR.LE.1.D0) GO TO 60 IF (ERR.GT.((NJ(K+1)*NJ(K))/4.D0)**2) GO TO 100 50 CONTINUE CALL STOERM(K,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W, 1 ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3, 2 REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE, 3 IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN) IF (ATOV) GO TO 10 KC=K IF (ERR.LE.1.D0) GO TO 60 C --- HOPE FOR CONVERGENCE IN LINE K+1 55 CONTINUE IF (ERR.GT.(NJ(K+1)/2.D0)**2) GO TO 100 KC=K+1 CALL STOERM(KC,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W, 1 ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3, 2 REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE, 3 IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN) IF (ATOV) GO TO 10 IF (ERR.GT.1.D0) GO TO 100 C --- STEP IS ACCEPTED 60 XOLD=X X=X+H IF (IOUT.GE.2) THEN C --- KMIT = MU OF THE PAPER KMIT=MAX(1,2*KC-MUDIF+1) DO 478 I=1,N 478 YH1(I)=T(1,I) CALL FCN(N,X,YH1,YH2,RPAR,IPAR) H2=H*H DO 69 I=1,NRD DENS(I)=Y(ICOMP(I)) DENS(NRD+I)=H*YP(ICOMP(I)) DENS(2*NRD+I)=H2*DZ(ICOMP(I)) DENS(3*NRD+I)=YH1(ICOMP(I)) DENS(4*NRD+I)=H*TP(1,ICOMP(I)) DENS(5*NRD+I)=H2*YH2(ICOMP(I)) 69 CONTINUE XOLDD=XOLD HHH=H C --- COMPUTE SOLUTION AND FIRST DERIVATIVE AT MID-POINT ---- DO 473 J=2,KC DBLENJ=NJ(J) DO 473 L=J,2,-1 FACTOR=(DBLENJ/NJ(L-1))**2-1.D0 DO 473 I=1,NRD YSAFE(L-1,I)=YSAFE(L,I)+(YSAFE(L,I)-YSAFE(L-1,I))/FACTOR YSAFEP(L-1,I)=YSAFEP(L,I)+(YSAFEP(L,I)-YSAFEP(L-1,I))/FACTOR 473 CONTINUE DO 474 I=1,NRD DENS(6*NRD+I)=YSAFE(1,I) 474 DENS(7*NRD+I)=H*YSAFEP(1,I) C DO 180 KMI=2,KMIT C --- COMPUTE KMI-TH DERIVATIVE AT MID-POINT ---- KBEG=KMI/2 IF (KMI.EQ.2*KBEG) THEN IF (KMI.EQ.2) THEN DO 570 I=1,NRD 570 YSAFE(1,I)=(DZ(ICOMP(I))+FSAFE(1,I))/2.D0 KBEG=2 END IF DO 375 KK=KBEG,KC FACNJ=0.5D0*(NJ(KK)/2.D0)**(KMI-2) IPT=(KK-1)**2+KK+KMI/2-2 DO 371 I=1,NRD 371 YSAFE(KK,I)=(FSAFE(IPT,I)+FSAFE(IPT+1,I))*FACNJ 375 CONTINUE ELSE DO 575 KK=KBEG,KC FACNJ=(NJ(KK)/2.D0)**(KMI-2) IPT=(KK-1)**2+KK+KBEG-1 DO 571 I=1,NRD 571 YSAFE(KK,I)=FSAFE(IPT,I)*FACNJ 575 CONTINUE END IF C------ EXTRAPOLATION DO 373 J=KBEG+1,KC DBLENJ=NJ(J) DO 373 L=J,KBEG+1,-1 FACTOR=(DBLENJ/NJ(L-1))**2-1.D0 DO 373 I=1,NRD YSAFE(L-1,I)=YSAFE(L,I)+(YSAFE(L,I)-YSAFE(L-1,I))/FACTOR 373 CONTINUE KRN=(KMI+6)*NRD DO 374 I=1,NRD 374 DENS(KRN+I)=YSAFE(KBEG,I)*H2 IF (KMI.EQ.KMIT) GOTO 180 C ------ COMPUTE DIFFERENCES DO 66 KK=(KMI+1)/2,KC LBEG=(KK-1)**2+KMI-1 LEND=KK**2 IF (KMI.EQ.2) LBEG=LBEG+1 DO 64 L=LEND,LBEG,-1 DO 64 I=1,NRD 64 FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I) IF (KMI.EQ.2) THEN L=LBEG-1 DO 65 I=1,NRD 65 FSAFE(L,I)=FSAFE(L,I)-DZ(ICOMP(I)) END IF 66 CONTINUE 180 CONTINUE CALL INTPO2(NRD,DENS,KMIT) C --- ESTIMATION OF INTERPOLATION ERROR IF (IDERR.EQ.0.AND.KMIT.GE.1) THEN ERRINT=0.D0 DO 187 I=1,NRD 187 ERRINT=ERRINT+(DENS((KMIT+6)*NRD+I)/SCAL(ICOMP(I)))**2 ERRINT=SQRT(ERRINT/NRD)*ERRFAC(KMIT) HOPTDE=H/MAX((ERRINT)**(1.D0/(KMIT+6)),0.01D0) IF (ERRINT.GT.10.D0) THEN H=HOPTDE X=XOLD NREJCT=NREJCT+1 REJECT=.TRUE. GOTO 10 END IF END IF DO 189 I=1,N 189 DZ(I)=YH2(I) END IF DO 70 I=1,N Y(I)=T(1,I) 70 YP(I)=TP(1,I) NACCPT=NACCPT+1 IF (IOUT.GE.1) THEN CALL SOLOUT (NACCPT+1,XOLD,X,Y,YP,N,DENS,NCOM,ICOMP,NRD, & RPAR,IPAR,IRTRN) IF (IRTRN.LT.0) GOTO 120 END IF C --- COMPUTE OPTIMAL ORDER IF (KC.EQ.2) THEN KOPT=MIN(3,KM-1) IF (REJECT) KOPT=2 GO TO 80 END IF IF (KC.LE.K) THEN KOPT=KC IF (W(KC-1).LT.W(KC)*FAC3) KOPT=KC-1 IF (W(KC).LT.W(KC-1)*FAC4) KOPT=MIN(KC+1,KM-1) ELSE KOPT=KC-1 IF (KC.GT.3.AND.W(KC-2).LT.W(KC-1)*FAC3) KOPT=KC-2 IF (W(KC).LT.W(KOPT)*FAC4) KOPT=MIN(KC,KM-1) END IF C --- AFTER A REJECTED STEP 80 IF (REJECT) THEN K=MIN(KOPT,KC) H=POSNEG*MIN(ABS(H),ABS(HH(K))) REJECT=.FALSE. GO TO 10 END IF C --- COMPUTE STEPSIZE FOR NEXT STEP IF (KOPT.LE.KC) THEN H=HH(KOPT) ELSE IF (KC.LT.K.AND.W(KC).LT.W(KC-1)*FAC4) THEN H=HH(KC)*A(KOPT+1)/A(KC) ELSE H=HH(KC)*A(KOPT)/A(KC) END IF END IF K=KOPT H=POSNEG*ABS(H) IF (LAST) GOTO 110 GO TO 10 C --- STEP IS REJECTED 100 CONTINUE K=MIN(K,KC,KM-1) IF (K.GT.2.AND.W(K-1).LT.W(K)*FAC3) K=K-1 NREJCT=NREJCT+1 H=POSNEG*HH(K) REJECT=.TRUE. LAST=.FALSE. GO TO 30 C --- SOLUTION EXIT 110 CONTINUE IDID=1 RETURN C --- FAIL EXIT 120 WRITE (6,979) X,H 979 FORMAT(' EXIT OF ODEX2 AT X=',D14.7,' H=',D14.7) IDID=-1 RETURN END C SUBROUTINE STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ, 1 HH,W,ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV, 2 SAFE3,REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE, 3 LFSAFE,IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN) C --- THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE C --- EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATION C --- OF THE OPTIMAL STEPSIZE IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL FCN LOGICAL REJECT,ATOV DIMENSION Y(N),YP(N),DY(N),YH1(N),YH2(N),DZ(N),SCAL(N),SCALP(N) DIMENSION T(KM,N),TP(KM,N),NJ(KM),HH(KM),W(KM),A(KM) DIMENSION RTOL(*),ATOL(*),RPAR(*),IPAR(*),ICOMP(NRD) DIMENSION FSAFE(LFSAFE,NRD),YSAFE(KM,NRD),YSAFEP(KM,NRD) HJ=H/NJ(J) HJ2=HJ*2.D0 C --- EULER STARTING STEP DO 30 I=1,N YH1(I)=Y(I) 30 YH2(I)=YP(I)+HJ*DZ(I) C --- EXPLICIT MIDPOINT (STOERMER) RULE M=NJ(J)/2 DO 35 MM=1,M IF (IOUT.GE.2.AND.MM.EQ.(M+1)/2) THEN DO 31 I=1,NRD YSAFEP(J,I)=YH2(ICOMP(I)) 31 YSAFE(J,I) =YH1(ICOMP(I))+HJ*YH2(ICOMP(I)) END IF DO 33 I=1,N 33 YH1(I)=YH1(I)+HJ2*YH2(I) CALL FCN(N,X+HJ2*MM,YH1,DY,RPAR,IPAR) IF (IOUT.GE.2.AND.ABS(MM-(M+1)/2).LT.J) THEN IPT=IPT+1 DO 32 I=1,NRD 32 FSAFE(IPT,I)=DY(ICOMP(I)) END IF IF (MM.EQ.M) GOTO 37 DO 34 I=1,N 34 YH2(I)=YH2(I)+HJ2*DY(I) 35 CONTINUE C --- FINAL SMOOTHING STEP 37 CONTINUE DO 40 I=1,N T(J,I)=YH1(I) 40 TP(J,I)=YH2(I)+HJ*DY(I) NFCN=NFCN+M C --- POLYNOMIAL EXTRAPOLATION IF (J.EQ.1) RETURN DBLENJ=NJ(J) DO 60 L=J,2,-1 FAC=(DBLENJ/NJ(L-1))**2-1.D0 DO 60 I=1,N T(L-1,I)=T(L,I)+(T(L,I)-T(L-1,I))/FAC TP(L-1,I)=TP(L,I)+(TP(L,I)-TP(L-1,I))/FAC 60 CONTINUE ERR=0.D0 C --- SCALING IF (ITOL.EQ.0) THEN DO 65 I=1,N SCAL(I)=ATOL(1)+RTOL(1)*MAX(ABS(Y(I)),ABS(T(1,I))) 65 SCALP(I)=ATOL(1)+RTOL(1)*MAX(ABS(YP(I)),ABS(TP(1,I))) ELSE DO 66 I=1,N SCAL(I)=ATOL(I)+RTOL(I)*MAX(ABS(Y(I)),ABS(T(1,I))) 66 SCALP(I)=ATOL(I+N)+RTOL(I+N)*MAX(ABS(YP(I)),ABS(TP(1,I))) END IF DO 67 I=1,N 67 ERR=ERR+((T(1,I)-T(2,I))/SCAL(I))**2 & +((TP(1,I)-TP(2,I))/SCALP(I))**2 ERR=SQRT(ERR/(2*N)) IF (ERR*UROUND.GE.1.D0) GOTO 79 IF (J.GT.2.AND.ERR.GE.ERROLD) GOTO 79 ERROLD=MAX(4*ERR,1.D0) C --- COMPUTE OPTIMAL STEPSIZES EXPO=1.D0/(2*J-1) FACMIN=FAC1**EXPO FAC=MIN(FAC2/FACMIN,MAX(FACMIN,(ERR/SAFE1)**EXPO/SAFE2)) FAC=1.D0/FAC HH(J)=MIN(ABS(H)*FAC,HMAX) W(J)=A(J)/HH(J) RETURN 79 ATOV=.TRUE. H=H*SAFE3 REJECT=.TRUE. RETURN END C SUBROUTINE INTPO2(N,Y,IMIT) C --- COMPUTES THE COEFFICIENTS OF THE INTERPOLATION FORMULA IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N*(IMIT+7)),A(0:40) C --- BEGIN WITH HERMITE INTERPOLATION (ORDER 5) DO 100 I=1,N Y0=Y(I) Y1=Y(3*N+I) YP0=Y(N+I) YP1=Y(4*N+I) YPP0=Y(2*N+I)/2.D0 YPP1=Y(5*N+I)/2.D0 YDIFF=Y1-Y0 AH=YDIFF-YP0 BH=YP1-YDIFF CH=AH-YPP0 DH=BH-AH EH=YPP1-BH FH=DH-CH GH=EH-DH Y(N+I)=YDIFF Y(2*N+I)=-AH Y(3*N+I)=-DH Y(4*N+I)=GH Y(5*N+I)=FH IF (IMIT.LT.0) GOTO 100 C --- COMPUTE THE DERIVATIVES OF HERMITE AT MIDPOINT ABH=AH+BH GFH=GH+FH GMF=GH-FH PH0=0.5D0*(Y0+Y1+0.25D0*(-ABH+0.25D0*GFH)) PH1=YDIFF+0.25D0*(AH-BH+0.25D0*GMF) PH2=ABH-0.5D0*GFH PH3=6.D0*(BH-AH)-3.D0*GMF PH4=12.D0*GFH PH5=120.D0*GMF C --- COMPUTE THE FURTHER COEFFICIENTS IF (IMIT.LT.1) GOTO 20 A(1)=64.D0*(Y(7*N+I)-PH1) IF (IMIT.LT.3) GOTO 20 A(3)=64.D0*(Y(9*N+I)-PH3+A(1)*9.D0/8.D0) IF (IMIT.LT.5) GOTO 20 A(5)=64.D0*(Y(11*N+I)-PH5+A(3)*15.D0/4.D0-A(1)*90.D0) IF (IMIT.LT.7) GOTO 20 DO 10 IM=7,IMIT,2 FC1=IM*(IM-1)*3.D0/16.D0 FC2=FC1*(IM-2)*(IM-3)*4.D0 FC3=IM*(IM-1)*(IM-2)*(IM-3)*(IM-4)*(IM-5) 10 A(IM)=64.D0*(Y((IM+6)*N+I)+FC1*A(IM-2)-FC2*A(IM-4)+FC3*A(IM-6)) 20 CONTINUE A(0)=64.D0*(Y(6*N+I)-PH0) IF (IMIT.LT.2) GOTO 60 A(2)=64.D0*(Y(N*8+I)-PH2+A(0)*3.D0/8.D0) IF (IMIT.LT.4) GOTO 60 A(4)=64.D0*(Y(N*10+I)-PH4+A(2)*9.D0/4.D0-A(0)*18.D0) IF (IMIT.LT.6) GOTO 60 DO 30 IM=6,IMIT,2 FC1=IM*(IM-1)*3.D0/16.D0 FC2=FC1*(IM-2)*(IM-3)*4.D0 FC3=IM*(IM-1)*(IM-2)*(IM-3)*(IM-4)*(IM-5) 30 A(IM)=64.D0*(Y(N*(IM+6)+I)+A(IM-2)*FC1-A(IM-4)*FC2+A(IM-6)*FC3) 60 CONTINUE DO 70 IM=0,IMIT 70 Y(N*(IM+6)+I)=A(IM) 100 CONTINUE RETURN END C FUNCTION CONTX2(II,X,Y,NCON,ICOMP,N) C ---------------------------------------------------------- C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION C WITH THE OUTPUT-SUBROUTINE FOR ODEX2. IT PROVIDES AN C APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X. C ---------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(NCON),ICOMP(N) COMMON /CONOD2/XOLD,H,IMIT C ----- COMPUTE PLACE OF II-TH COMPONENT I=0 DO 5 J=1,N IF (ICOMP(J).EQ.II) I=J 5 CONTINUE IF (I.EQ.0) THEN WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',II RETURN END IF C ----- COMPUTE THE INTERPOLATED VALUE THETA=(X-XOLD)/H THETA1=1.D0-THETA PHTHET=Y(I)+THETA*(Y(N+I)+THETA1*(Y(2*N+I)+THETA*(Y(3*N+I)+ & THETA1*(Y(4*N+I)*THETA+Y(5*N+I)*THETA1)))) IF (IMIT.LT.0) THEN CONTX2=PHTHET RETURN END IF THETAH=THETA-0.5D0 CONTX2=Y(N*(IMIT+6)+I) DO 70 IM=IMIT,1,-1 70 CONTX2=Y(N*(IM+5)+I)+CONTX2*THETAH/IM CONTX2=PHTHET+(THETA*THETA1)**3*CONTX2 RETURN END