C ALGORITHM 700 , COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 4, DECEMBER, 1991, PP. 500-501. PROGRAM MAIN C C This sample driver calls ZCOUNT and SLEIGN to find the smallest C eigenvalue, EIG(1), of the semidefinite Sturm-Liouville problem C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (-pi,pi) C C where P(x) = 1, Q(x) = 4, and R(x) = 1 for x .le. -PI/2 or C x .gt. PI/2, and R(x) = 0 otherwise. The correct answer is C EIG(1) = -3, more generally, EIG(n) = n**2 - 4. C INTEGER I,IFLAG,INTAB,ISLFUN,JPAIRS,JSUM,NUMEIG,NWRITE REAL A,A1,A2,B,B1,B2,EIG,P0ATA,P0ATB,QFATA,QFATB,TOL REAL PAIRS(2),SLFUN(9) EXTERNAL PARAMS,SLEIGN,ZCOUNT C DATA NWRITE/6/ C CALL PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) C C Count number of zeros in (A,B) of eigenfunction associated C with smallest eigenvalue. C CALL ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) WRITE (NWRITE,10) ' zero count =',JSUM 10 FORMAT (A,I6) C C Calculate smallest eigenvalue and associated eigenfunction data. C NUMEIG = 1 + JSUM EIG = 0.0 TOL = 1.0E-3 ISLFUN = 0 CALL SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) WRITE (NWRITE,20) ' iflag,numeig,eig,tol = ',IFLAG,NUMEIG,EIG,TOL 20 FORMAT (A,2I6,F12.5,1PE12.2) WRITE (NWRITE,30) ' slfun(1-9) =',(SLFUN(I),I=1,9) 30 FORMAT (A,F12.5/(13X,3F12.5)) STOP END C SUBROUTINE PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) INTEGER INTAB,JPAIRS REAL A,A1,A2,B,B1,B2,ONE,PI,P0ATA,P0ATB,QFATA,QFATB REAL PAIRS(*) INTRINSIC ATAN COMMON PI ONE = 1.0 PI = 4.0*ATAN(ONE) A = -PI B = PI INTAB = 1 P0ATA = -1.0 QFATA = 1.0 P0ATB = -1.0 QFATB = 1.0 A1 = 1.0 A2 = 0.0 B1 = 1.0 B2 = 0.0 JPAIRS = 1 PAIRS(1) = -PI/2.0 PAIRS(2) = PI/2.0 RETURN END C FUNCTION P(X) REAL P,X P = 1.0 RETURN END C FUNCTION Q(X) REAL Q,X Q = 4.0 RETURN END C FUNCTION R(X) REAL PI,R,X COMMON PI R = 0.0 IF (X.LE.-PI/2.0 .OR. X.GT.PI/2.0) R = 1.0 RETURN END SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) INTEGER INTAB,NUMEIG,IFLAG,ISLFUN REAL A,B,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2,EIG,TOL REAL SLFUN(9) C ********** C C This subroutine is designed for the calculation of a specified C eigenvalue, EIG, of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C The problem may be either regular or singular. In the C regular case, boundary conditions of the form C C a1*y(a) + a2*p(a)*y'(a) = 0 C b1*y(b) + b2*p(b)*y'(b) = 0 C C are prescribed by specifying the numbers A1, A2, B1, and B2. C The index of the desired eigenvalue is specified in NUMEIG C and its requested accuracy in TOL. Initial data for the C associated eigenfunction are also computed along with values C at selected points, if desired, in array SLFUN. C C The SUBROUTINE statement is C C SUBROUTINE sleign(a,b,intab,p0ata,qfata,p0atb,qfatb,a1,a2,b1,b2, C numeig,eig,tol,iflag,islfun,slfun) C C where C C A and B are input variables defining the interval. If the C interval is finite, A must be less than B. (See INTAB below.) C C INTAB is an integer input variable specifying the nature of the C interval. It can have four values. C C INTAB = 1 - A and B are finite. C INTAB = 2 - A is finite and B is infinite (+). C INTAB = 3 - A is infinite (-) and B is finite. C INTAB = 4 - A is infinite (-) and B is infinite (+). C C If either A or B is infinite, it is classified singular and C its value is ignored. C C P0ATA, QFATA, P0ATB, and QFATB are input variables set to C 1.0 or -1.0 as the following properties of P, Q, and R at C the interval endpoints are true or false, respectively. C C P0ATA - p(a) is zero. (If true, A is singular.) C QFATA - q(a) and r(a) are finite. (If false, A is singular.) C P0ATB - p(b) is zero. (If true, B is singular.) C QFATB - q(b) and r(b) are finite. (If false, B is singular.) C C A1 and A2 are input variables set to prescribe the boundary C condition at A. Their values are ignored if A is singular. C C B1 and B2 are input variables set to prescribe the boundary C condition at B. Their values are ignored if B is singular. C C NUMEIG is an integer variable. On input, it should be set to C the index of the desired eigenvalue (increasing sequence). C On output, it is unchanged unless the problem (apparently) C lacks NUMEIG eigenvalues, in which case it is reset to the C index of the largest eigenvalue. C C EIG is a variable set on input to 0.0 or to an initial guess of C the eigenvalue. If EIG is set to 0.0, SLEIGN will generate C the initial guess. On output, EIG holds the calculated C eigenvalue if IFLAG (see below) signals success. C C TOL is a variable set on input to the desired accuracy of the C eigenvalue. On output, TOL is reset to the accuracy estimated C to have been achieved if IFLAG (see below) signals success. C This accuracy estimate is absolute if EIG is less than one C in magnitude, and relative otherwise. In addition, prefixing C TOL with a negative sign, removed after interrogation, serves C as a flag to request trace output from the calculation. C C IFLAG is an integer output variable set as follows. C C IFLAG = 1 - successful problem solution. C IFLAG = 2 - improper input parameters. C IFLAG = 3 - NUMEIG exceeds actual number of eigenvalues. C IFLAG = 4 - some uncertainty about accuracy estimate TOL. C IFLAG = 5 - convergence too slow, best results returned. C IFLAG = 6 - failure, integrator could not complete. C C ISLFUN is an integer input variable set to the number of C selected eigenfunction values desired. If no values are C desired, set ISLFUN zero or negative. C C SLFUN is an array of length at least 9. On output, the first 9 C locations contain the integration interval and initial data C that completely determine the eigenfunction. C C SLFUN(1) - point where two pieces of eigenfunction Y match. C SLFUN(2) - left endpoint XAA of the (truncated) interval. C SLFUN(3) - value of THETA at XAA. (Y = RHO*sin(THETA)) C SLFUN(4) - value of F at XAA. (RHO = exp(F)) C SLFUN(5) - right endpoint XBB of the (truncated) interval. C SLFUN(6) - value of THETA at XBB. C SLFUN(7) - value of F at XBB. C SLFUN(8) - final value of integration accuracy parameter EPS. C SLFUN(9) - the constant Z in the polar form transformation. C C F(XAA) and F(XBB) are chosen so that the eigenfunction is C continuous in the interval (XAA,XBB) and has weighted (by R) C L2-norm of 1.0 on the interval. If ISLFUN is positive, then C on input the further ISLFUN locations of SLFUN specify the C points, in ascending order, where the eigenfunction values C are desired and on output contain the values themselves. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... alfbet,dxdt,epslon,esteig,estpac,f,gerk C C This version dated 5/18/89. C Paul B. Bailey, Sandia Laboratories, Albuquerque C Burton S. Garbow, Argonne National Laboratory C C ********** C .. Scalars in Common .. INTEGER INTSAV LOGICAL EIGF REAL ASAV,BSAV,C1,C2,EIGSAV,Z C .. C .. Local Scalars .. INTEGER I,IA,IB,IMAX,IMIN,IOUT,IP,J,JFLAG,KFLAG,LFLAG,MF,ML, 1 NEIG,NMID LOGICAL AOK,BOK,BRACKT,CHNGAB,CHNGM,CONVRG,FYNYT,FYNYT1, 1 LIMIT,LOGIC,NEWTON,ONEDIG,PRIN,THEGT0,THELT0 REAL AA,AAA,ALFA,ASL,ASR,BB,BBB,BETA, 1 C,CHNG,CL,CR,DAV,DE,DEDW,DEN,DERIVL,DERIVR,DIST, 2 DPSIL,DPSIPL,DPSIPR,DPSIR,DT,DTHDA,DTHDAA,DTHDB, 3 DTHDBB,DTHDE,DTHDEA,DTHDEB,DTHETA,DTHOLD,E,EEE, 4 EIGLO,EIGLT,EIGRT,EIGUP,EL,ELIM,EMAX,EMIN,EOLD,EPS,EPSMIN, 5 ER1,ER2,ESTERR,FLO,FMAX,FUP,GMAX,GUESS,H,ONE,PI,PIN, 6 PSIL,PSIPL,PSIPR,PSIR,PX,QAV,QX,RATIO,RAV,RAY,RX, 7 SL,SQL,SQR,SR,T,T1,T2,T3,TAU,TEMP,THRESH,TMAX,TMID,TMP, 8 U,UL,UR,V,WL,X,X50,XAA,XBB,XBC,XMID,XSAV,ZAV C .. C .. Local Arrays .. INTEGER IWORK(5) REAL DS(98),ERL(3),ERR(3),PS(99),QS(99),RS(99), 1 WORK(27),YL(3),YR(3) C .. C .. External Functions .. REAL EPSLON,P,Q,R EXTERNAL EPSLON,P,Q,R C .. C .. External Subroutines .. EXTERNAL ALFBET,DXDT,ESTEIG,ESTPAC,F,GERK C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,COS,EXP,INT,LOG,MAX,MIN,SIGN,SIN,TAN C .. C .. Common blocks .. COMMON /DATADT/ASAV,BSAV,C1,C2,INTSAV COMMON /DATAF/EIGSAV,EIGF COMMON /ZEE/Z C .. C Set constants EPSMIN, the computer unit roundoff error, and PI. C (Variable ONE set to 1.0 eases precision conversion.) C ONE = 1.0 EPSMIN = EPSLON(ONE) PI = 4.0*ATAN(ONE) C C Set output device number. C IOUT = 6 C C Check input parameters for errors. If errors, return IFLAG=2. C LOGIC = TOL.NE.0.0 .AND. 1.LE.INTAB .AND. INTAB.LE.4 .AND. 1 P0ATA*QFATA*P0ATB*QFATB*NUMEIG.NE.0.0 IF (INTAB.EQ.1) LOGIC = LOGIC .AND. A.LT.B IF (.NOT.LOGIC) THEN IFLAG = 2 RETURN END IF C C Set PRIN = .true. to trigger trace printout of successive steps. C PRIN = .FALSE. IF (TOL.LT.0.0) PRIN = .TRUE. C C Set EPS to the (initial) integration accuracy. C EPS = 0.001 C C AOK (BOK) signals, if true, that endpoint A (B) is not singular. C AOK = INTAB.LT.3 .AND. P0ATA.LT.0.0 .AND. QFATA.GT.0.0 BOK = (INTAB.EQ.1 .OR. INTAB.EQ.3) .AND. 1 P0ATB.LT.0.0 .AND. QFATB.GT.0.0 NEIG = ABS(NUMEIG) - 1 C C Initial C1 and C2, used in the mapping between X and T intervals. C C1 = 1.0 C2 = 0.0 C DO (SAVE-INPUT-DATA) ASAV = A BSAV = B INTSAV = INTAB TAU = ABS(TOL) C END (SAVE-INPUT-DATA) C C Evaluate P, Q, R to obtain preliminary information about the C differential equation. C C DO (SAMPLE-COEFFICIENTS) THRESH = 1.0E+17 10 CONTINUE CALL DXDT(EPSMIN,TEMP,X50) PX = P(X50) QX = Q(X50) RX = R(X50) PS(50) = PX QS(50) = QX/PX RS(50) = RX/PX C C DAV,QAV,RAV are used in special case estimation when NUMEIG = 1,2. C EMIN = min(-Q/R), achieved at X for index value IMIN. C EMAX = max(-Q/R), achieved at X for index value IMAX. C MF and ML are the least and greatest index values, respectively. C DAV = 0.0 QAV = 0.0 RAV = 0.0 XSAV = X50 EMIN = THRESH EMAX = -THRESH IF (RX.NE.0.0) THEN EMIN = -QX/RX EMAX = EMIN IMIN = 50 IMAX = 50 END IF H = 0.9/40.0 DO 20 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I) = XSAV - X DAV = DAV + DS(I) QAV = QAV + DS(I)*(0.5*(QS(I)+QS(I+1))-QAV)/DAV RAV = RAV + DS(I)*(0.5*(RS(I)+RS(I+1))-RAV)/DAV XSAV = X C C Try to avoid overflow by stopping when functions are large near A. C FYNYT = (ABS(RX)+ABS(QX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF MF = I IF (.NOT.FYNYT) GO TO 30 20 CONTINUE 30 CONTINUE AAA = T IF (AOK) AAA = -1.0 XSAV = X50 DO 40 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I-1) = X - XSAV DAV = DAV + DS(I-1) QAV = QAV + DS(I-1)*(0.5*(QS(I-1)+QS(I))-QAV)/DAV RAV = RAV + DS(I-1)*(0.5*(RS(I-1)+RS(I))-RAV)/DAV XSAV = X C C Try to avoid overflow by stopping when functions are large near B. C FYNYT1 = (ABS(QX)+ABS(RX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF ML = I - 1 IF (.NOT.FYNYT1) GO TO 50 40 CONTINUE 50 CONTINUE BBB = T IF (BOK) BBB = 1.0 LOGIC = C1.EQ.1.0 .AND. (.NOT.FYNYT .OR. .NOT.FYNYT1) C C Modify (T,X) transformation corresponding to truncated interval. C IF (LOGIC) THEN C1 = 0.5*(BBB-AAA) C2 = 0.5*(AAA+BBB) GO TO 10 END IF C C Estimate upper bound ELIM for EIG such that boundary conditions C can be satisfied. C ELIM = EMAX + 1.0 IF (INTAB.EQ.3 .OR. (P0ATA.GT.0.0 .AND. QFATA.LT.0.0)) THEN IF (-QS(MF)/RS(MF).LE.ELIM) THEN ELIM = -QS(MF)/RS(MF) IMAX = MF END IF END IF IF (INTAB.EQ.2 .OR. (P0ATB.GT.0.0 .AND. QFATB.LT.0.0)) THEN IF (-QS(ML)/RS(ML).LE.ELIM) THEN ELIM = -QS(ML)/RS(ML) IMAX = ML END IF END IF IF (INTAB.EQ.4) THEN ELIM = MIN(ELIM,-QS(MF)/RS(MF),-QS(ML)/RS(ML)) IF (-QS(MF)/RS(MF).EQ.ELIM) IMAX = MF IF (-QS(ML)/RS(ML).EQ.ELIM) IMAX = ML END IF ELIM = ELIM - EPSMIN IF (ELIM.EQ.0.0) ELIM = -EPSMIN LIMIT = ELIM.LE.EMAX C END (SAMPLE-COEFFICIENTS) PIN = (NEIG+1)*PI IF (EIG.EQ.0.0) THEN C DO (ESTIMATE-EIG) CALL ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) C END (ESTIMATE-EIG) END IF GUESS = EIG C DO (SET-INITIAL-INTERVAL-AND-MATCHPOINT) IF (GUESS.NE.0.0) THEN C C Reduce overly large guess for EIG to upper bound if necessary. C IF (LIMIT .AND. EIG.GT.ELIM) EIG = ELIM EEE = EIG C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,TEMP,U,TMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) END IF C C Choose initial interval as large as possible that avoids overflow. C IA and IB are boundary indices for nonnegativity of EIG*R+Q. C AA = -1.0 IF (.NOT.AOK) THEN AA = H*(IA-50) IF (IA.LT.10) AA = -1.0 + 0.1/4.0**(10-IA) END IF BB = 1.0 IF (.NOT.BOK) THEN BB = H*(IB-50) IF (IB.GT.90) BB = 1.0 - 0.1/4.0**(IB-90) END IF AA = AA + 0.6*(AAA-AA) BB = BB + 0.6*(BBB-BB) C C Determine boundary values ALFA and BETA for theta at A and B. C Z = 1.0 CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA C C Take boundary conditions into account in estimation of EIG. C PIN = PIN + BETA - ALFA - PI IF (GUESS.EQ.0.0) EEE = EL + DEDW*(PIN-WL) C C Subroutine ESTPAC must be called again because PIN has changed. C C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,TEMP,U,ZAV) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C Choose the constant Z. C IF (U.GT.0.0) Z = ZAV/U C C Reset boundary values ALFA and BETA. C CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' alfa=',ALFA,' beta=',BETA C C Special case formula for estimation of EIG when NUMEIG = 1,2. C IF (U.EQ.0.0 .AND. NEIG.LE.1 .AND. (BETA+NEIG*PI).LT.ALFA) THEN XBC = MAX(-1.0/TAN(ALFA),1.0/TAN(BETA)) EEE = -(XBC*XBC-QAV)/RAV DEDW = XBC*(1.0+XBC*XBC)/RAV END IF C C Choose initial matching point TMID. C TMID = H*(IP-50) IF (TMID.LT.-0.8) TMID = -0.4 IF (TMID.GT.0.8) TMID = 0.4 IF (PRIN) WRITE(IOUT,'(A,E15.7,A,F11.8,A,E15.7)') 1 ' estim=',EEE,' tmid=',TMID,' z=',Z IF (PRIN) WRITE(IOUT,'(A,F11.8,A,F11.8,A,F11.8,A,F11.8)') 1 ' aaa=',AAA,' aa=',AA,' bb=',BB,' bbb=',BBB C END (SET-INITIAL-INTERVAL-AND-MATCHPOINT) C C End preliminary work, begin main task of computing EIG. C C Logical variables have the following meanings if true. C AOK - endpoint A is not singular. C BOK - endpoint B is not singular. C CHNGM - matching point TMID should be changed. C CHNGAB - one or both endpoints should be moved farther out. C BRACKT - EIG has been bracketed. C CONVRG - convergence test for EIG has been successfully passed. C NEWTON - Newton iteration may be employed. C THELT0 - lower bound for EIG has been found. C THEGT0 - upper bound for EIG has been found. C LIMIT - upper bound exists with boundary conditions satisfied. C ONEDIG - most significant digit can be expected to be correct. C EIGF - derivative argument is in original coordinate system. C EIG = EEE EIGF = .FALSE. CHNGM = .FALSE. CHNGAB = .TRUE. 60 CONTINUE IF (CHNGAB) THEN C DO (INITIAL-IZE) CHNGAB = .FALSE. BRACKT = .FALSE. CONVRG = .FALSE. THELT0 = .FALSE. THEGT0 = .FALSE. EIGLO = 0.0 EIGLT = 0.0 EIGRT = 0.0 EIGUP = 0.0 DTHOLD = 0.0 C END (INITIAL-IZE) END IF C C Recompute boundary conditions at singular endpoint(s). C C DO (RESET-BOUNDARY-CONDITIONS) DERIVL = 0.0 IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,ALFA,KFLAG,DERIVL) DERIVR = 0.0 IF (.NOT.BOK) THEN CALL ALFBET(B,INTAB,BB,B1,B2,EIG,P0ATB,QFATB,.FALSE., 1 BETA,JFLAG,DERIVR) BETA = PI - BETA END IF IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' alfa=',ALFA,' beta=',BETA C END (RESET-BOUNDARY-CONDITIONS) 70 CONTINUE IF (EIG.NE.GUESS .AND. .NOT.BRACKT) THEN C C If initial guess was supplied, check that boundary conditions C can be satisfied at singular endpoints. If not, try for C slightly lower EIG consistent with boundary conditions. C 80 CONTINUE IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,TMP,KFLAG,TEMP) IF (.NOT.BOK) CALL ALFBET(B,INTAB,BB,B1,B2,EIG, 1 P0ATB,QFATB,.FALSE.,TMP,JFLAG,TEMP) IF (KFLAG*JFLAG.NE.1) THEN IF (THEGT0) EIG = 0.6*EIG + 0.4*EIGUP IF (THELT0) EIG = 0.6*EIG + 0.4*EIGLO IF (THELT0 .AND. EIGLO.LT.ELIM) EIGUP = ELIM GO TO 80 END IF END IF C DO (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) 90 CONTINUE IF (PRIN) WRITE(IOUT,'(/A,E22.14,A,E10.3,A,E10.3)') 1 ' guess=',EIG,' eps=',EPS,' tmid=',TMID C DO (INTEGRATE-FOR-DTHETA) C DO (SET-INITIAL-CONDITIONS) DTHDEA = DERIVL DTHDAA = 0.0 IF (.NOT.AOK) THEN CALL DXDT(AA,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDAA = -(COS(ALFA)**2/PX + 1 C*SIN(ALFA)**2)*DT C C Two special cases for DTHDAA. C IF (C.GE.0.0 .AND. P0ATA.LT.0.0 .AND. 1 QFATA.LT.0.0) DTHDAA = DTHDAA + 1 ALFA*DT/(X-A) IF (C.GE.0.0 .AND. P0ATA.GT.0.0 .AND. 1 QFATA.GT.0.0) DTHDAA = DTHDAA + 2 (ALFA-0.5*PI)*DT/(X-A) END IF DTHDEB = -DERIVR DTHDBB = 0.0 IF (.NOT.BOK) THEN CALL DXDT(BB,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDBB = -(COS(BETA)**2/PX + 1 C*SIN(BETA)**2)*DT C C Two special cases for DTHDBB. C IF (C.GE.0.0 .AND. P0ATB.LT.0.0 .AND. 1 QFATB.LT.0.0) DTHDBB = DTHDBB + 2 (PI-BETA)*DT/(B-X) IF (C.GE.0.0 .AND. P0ATB.GT.0.0 .AND. 1 QFATB.GT.0.0) DTHDBB = DTHDBB + 2 (0.5*PI-BETA)*DT/(B-X) END IF TMAX = TMID GMAX = ABS(DTHDEA) EIGSAV = EIG C END (SET-INITIAL-CONDITIONS) C T C YL = (theta,d(theta)/d(eig),d(theta)/da) C T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 C C Use integrator in one-step mode towards change to different TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 100 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.GT.-0.8 .AND. 1 ABS(YL(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YL(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 100 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDA = DTHDAA*EXP(-2.0*YL(3)) C T C YR = (theta,d(theta)/d(eig),d(theta)/db) C T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 C C Use integrator in one-step mode towards change to different TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 110 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.LT.0.8 .AND. 1 ABS(YR(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YR(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 110 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDB = DTHDBB*EXP(-2.0*YR(3)) C C DTHETA measures theta difference from left and right integrations. C DTHETA = YL(1) - YR(1) DTHDE = YL(2) - YR(2) ER1 = ERL(1) - ERR(1) ER2 = ERL(2) - ERR(2) TMID = TMAX C END (INTEGRATE-FOR-DTHETA) C C Define ONEDIG to try to be sure of one correct digit in DTHETA. C Redo integrations with tighter tolerance until ONEDIG is true. C ONEDIG = ABS(ER1).LE.0.1*ABS(DTHETA) .AND. 1 ABS(ER2).LE.0.1*ABS(DTHDE) NEWTON = ABS(DTHETA).LT.0.06 IF (NEWTON) THEN C DO (COMPUTE-CONVRG) C C Measure convergence after adding separate contributions to error. C T1 = ABS(DTHETA)+50.0*ABS(ER1) T2 = (1.0+AA)*ABS(DTHDA) T3 = (1.0-BB)*ABS(DTHDB) ESTERR = (T1+T2+T3)/ABS(DTHDE)/MAX(ONE,ABS(EIG)) CONVRG = ESTERR.LE.TAU IF (PRIN) WRITE(IOUT,'(A,L2)') 1 ' converge=',CONVRG IF (PRIN .AND. .NOT.CONVRG) 1 WRITE(IOUT,'(A,E15.7)') 2 ' estim. acc.=',ESTERR C END (COMPUTE-CONVRG) END IF IF (.NOT.ONEDIG .OR. 1 ABS(ER1).GT.0.01*ABS(DTHETA)) THEN C C Reduce local error criterion, but return IFLAG=5 if too small. C EPS = 0.05*EPS IF (EPS.LE.EPSMIN) THEN IFLAG = 5 RETURN END IF END IF IF (.NOT.(ONEDIG .OR. CONVRG)) GO TO 90 IF (ABS(DTHETA).LT.0.1) CHNGM = .FALSE. IF (PRIN) WRITE(IOUT,'(A,E15.7,A,E15.7)') 1 ' dtheta=',DTHETA,' dthde=',DTHDE IF (PRIN) WRITE(IOUT,'(/A,E15.7,A,E15.7)') 1 ' thetal=',YL(1),' thetar=',YR(1) C END (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) C DO (SET-BRACKET-AND-FORM-NEWTON-ITERATES) C C EIG is bracketed when both THEGT0=.true. and THELT0=.true. C IF (DTHETA*DTHDE.GT.0.0) THEN IF (.NOT.THEGT0 .OR. EIG.LE.EIGUP) THEN THEGT0 = .TRUE. EIGUP = EIG FUP = DTHETA EIGRT = EIG - DTHETA/DTHDE END IF ELSE IF (.NOT.THELT0 .OR. EIG.GE.EIGLO) THEN THELT0 = .TRUE. EIGLO = EIG FLO = DTHETA EIGLT = EIG - DTHETA/DTHDE END IF END IF BRACKT = THELT0 .AND. THEGT0 IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' eigrt=',EIGRT,' eigup=',EIGUP IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' eiglt=',EIGLT,' eiglo=',EIGLO C END (SET-BRACKET-AND-FORM-NEWTON-ITERATES) IF (NEWTON) THEN CHNGM = .TRUE. IF (CONVRG) THEN CHNG = DTHDA*(-1.0-AA) - DTHDB*(1.0-BB) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) ELSE CHNGAB = T1.LT.0.5*(T2+T3) C C Move endpoint(s) out or take Newton step, according to CHNGAB. C IF (CHNGAB) THEN C DO (MOVE-ENDPOINTS) IF (T2.GT.T1 .AND. AA.GT.AAA) 1 AA = AA + 0.8*(-1.0-AA) IF (T3.GE.T1 .AND. BB.LT.BBB) 1 BB = BB + 0.8*(1.0-BB) AA = MAX(AA,AAA) BB = MIN(BB,BBB) IF ((AAA-AA).EQ.(BBB-BB)) THEN C C Cannot move endpoint(s) again. Store estimates and return IFLAG=5. C CHNG = (DTHDA-DTHDB)*(AAA-AA) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF EEE = EIG IF (PRIN) WRITE(IOUT,'(A,2E15.8)') 1 ' new endpoints ',AA,BB C END (MOVE-ENDPOINTS) ELSE EIG = EIG - DTHETA/DTHDE END IF END IF ELSE IF (BRACKT) THEN C C Obtain next estimate of EIG by bisection or linear interpolation. C FMAX = MAX(-FLO,FUP) EOLD = EIG RATIO = DTHETA/DTHOLD EIG = 0.5*(EIGLO+EIGUP) IF (FMAX.LE.1.5) THEN U = -FLO/(FUP-FLO) DIST = EIGUP - EIGLO EIG = EIGLO + U*DIST V = MIN(EIGLT,EIGRT) IF (EIG.LE.V) EIG = 0.5*(EIG+V) V = MAX(EIGLT,EIGRT) IF (EIG.GE.V) EIG = 0.5*(EIG+V) DE = EIG - EOLD CHNGAB = RATIO.GE.0.4 .AND. .NOT.(AOK.AND.BOK) IF (ABS(DE).LT.EPSMIN) THEN TOL = ABS(DE)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF END IF CHNGM = .NOT.CHNGM .AND. RATIO.GE.0.4 ELSE C DO (TRY-FOR-BRACKET) C C Take twice the Newton step in trying for a bracket. C IF (EIG.EQ.EEE) THEN IF (GUESS.NE.0.0) DEDW = 1.0/DTHDE CHNG = -(DEDW+1.0/DTHDE)*DTHETA IF (ABS(CHNG).GT.0.1*ABS(EIG)) 1 CHNG = -0.1*SIGN(EIG,DTHETA) ELSE CHNG = -2.0*DTHETA/DTHDE END IF LOGIC = EIG.NE.EEE .AND. DTHOLD.LT.0.0 .AND. 1 LIMIT .AND. CHNG.GT.(ELIM-EIG) IF (LOGIC) THEN CHNG = 0.99*(ELIM-EIG+EPSMIN) IF (CHNG.LT.EPSMIN) THEN C C If change is too small, EIG is presumed not to exist (IFLAG=3). C NUMEIG = NEIG - INT(-DTHETA/PI) IFLAG = 3 RETURN END IF C C Limit change in EIG to a factor of 10. C ELSE IF (ABS(CHNG).GT.(1.0+10.0*ABS(EIG))) THEN CHNG = SIGN(1.0+10.0*ABS(EIG),CHNG) ELSE IF (ABS(EIG).GE.1.0 .AND. 1 ABS(CHNG).LT.0.1*ABS(EIG)) THEN CHNG = 0.1*SIGN(EIG,CHNG) END IF EOLD = EIG EIG = EIG + CHNG C END (TRY-FOR-BRACKET) END IF END IF DTHOLD = DTHETA IF (.NOT.(CONVRG .OR. CHNGAB .OR. NEWTON)) GO TO 70 IF (.NOT.CONVRG) GO TO 60 IFLAG = 1 IF (PRIN) WRITE(IOUT,'(A,I7,A,E22.14,A,E10.3)') 1 ' numeig=',NUMEIG,' eig=',EIG,' tol=',TOL C DO (COMPUTE-EIGENFUNCTION-DATA) C C Convert from T to X values, fill 7 of first 9 locations of SLFUN. C CALL DXDT(TMID,TEMP,XMID) CALL DXDT(AA,TEMP,XAA) CALL DXDT(BB,TEMP,XBB) SLFUN(1) = XMID SLFUN(2) = XAA SLFUN(3) = ALFA SLFUN(5) = XBB SLFUN(6) = BETA + PI*NEIG SLFUN(8) = EPS SLFUN(9) = Z C C Compute SLFUN(4), SLFUN(7) towards normalizing the eigenfunction. C EIGSAV = EIG Z = -Z T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 LFLAG = 1 120 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 120 T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 LFLAG = 1 130 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 130 Z = -Z SL = SIN(YL(1)) SR = SIN(YR(1)) CL = COS(YL(1)) CR = COS(YR(1)) UL = (YL(2)-DTHDEA*EXP(-2.0*YL(3)))*Z UR = (YR(2)-DTHDEB*EXP(-2.0*YR(3)))*Z ASL = ABS(SL) ASR = ABS(SR) DEN = 0.5*LOG(UL*ASR*ASR-UR*ASL*ASL) SLFUN(4) = LOG(ASR) - YL(3) - DEN SLFUN(7) = LOG(ASL) - YR(3) - DEN C END (COMPUTE-EIGENFUNCTION-DATA) C DO (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) C C Perform final check on EIG. Return IFLAG=4 if not accurate enough. C E = ASR*EXP(-DEN) PSIL = E*SL PSIPL = E*CL SQL = E*E*UL DPSIL = PSIL*ERL(3) + PSIPL*ERL(1) DPSIPL = PSIL*ERL(1) + PSIPL*ERL(3) PSIPL = PSIPL*Z E = ASL*EXP(-DEN) PSIR = E*SR PSIPR = E*CR SQR = E*E*UR DPSIR = PSIR*ERR(3) + PSIPR*ERR(1) DPSIPR = PSIR*ERR(1) + PSIPR*ERR(3) PSIPR = PSIPR*Z RAY = EIG + (PSIL*PSIPL-PSIR*PSIPR)/(SQL-SQR) IF (PRIN) THEN WRITE(IOUT,'(A,E22.14)') ' ray=',RAY WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' psil=',PSIL,' psir=',PSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' psipl=',PSIPL,' psipr=',PSIPR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' sql=',SQL,' sqr=',SQR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' dpsil=',DPSIL,' dpsir=',DPSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' dpsipl=',DPSIPL,' dpsipr=',DPSIPR END IF C END (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) IF (ABS(RAY-EIG).GT.TAU*MAX(ONE,ABS(EIG))) IFLAG = 4 IF (ISLFUN.GT.0) THEN C C Calculate selected eigenfunction values by integration. C C DO (GENERATE-EIGENFUNCTION-VALUES) EIGF = .TRUE. NMID = 0 DO 140 I=1,ISLFUN IF (SLFUN(9+I).LE.SLFUN(1)) NMID = I 140 CONTINUE IF (NMID.GT.0) THEN X = SLFUN(2) YL(1) = SLFUN(3) YL(2) = 0.0 YL(3) = SLFUN(4) LFLAG = 1 DO 160 J=1,NMID 150 CONTINUE CALL GERK(F,3,YL,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 150 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YL(3))*SIN(YL(1)) 160 CONTINUE END IF IF (NMID.LT.ISLFUN) THEN X = SLFUN(5) YR(1) = SLFUN(6) YR(2) = 0.0 YR(3) = SLFUN(7) LFLAG = 1 DO 180 J=ISLFUN,NMID+1,-1 170 CONTINUE CALL GERK(F,3,YR,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 170 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YR(3))*SIN(YR(1)) 180 CONTINUE END IF C END (GENERATE-EIGENFUNCTION-VALUES) END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE ALFBET(XEND,INTAB,TT,COEF1,COEF2,EIG,P0,QF,OK, 1 VALUE,IFLAG,DERIV) INTEGER INTAB,IFLAG LOGICAL OK REAL XEND,TT,COEF1,COEF2,EIG,P0,QF,VALUE,DERIV C ********** C C This subroutine computes a boundary value for a specified endpoint C of the interval for a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is called C from SLEIGN. Both regular and singular endpoints are treated. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt,extrap C C ********** C .. Scalars in Common .. REAL Z C .. C .. Local Scalars .. LOGICAL LOGIC REAL C,CD,D,HH,ONE,PI,PX,QX,RX,T,TEMP,X C .. C .. External Functions .. REAL P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT,EXTRAP C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SIGN,SQRT C .. C .. Common blocks .. COMMON /ZEE/Z C .. C Set machine dependent constant. C C PI (variable ONE set to 1.0 eases precision conversion). ONE = 1.0 PI = 4.0*ATAN(ONE) C IFLAG = 1 DERIV = 0.0 IF (OK) THEN VALUE = 0.5*PI IF (COEF1.NE.0.0) VALUE = ATAN(-Z*COEF2/COEF1) LOGIC = (TT.LT.0.0 .AND. VALUE.LT.0.0) .OR. 1 (TT.GT.0.0 .AND. VALUE.LE.0.0) IF (LOGIC) VALUE = VALUE + PI ELSE LOGIC = (INTAB.EQ.2 .AND. TT.GT.0.0) .OR. 1 (INTAB.EQ.3 .AND. TT.LT.0.0) .OR. 2 INTAB.EQ.4 .OR. (P0.GT.0.0 .AND. QF.LT.0.0) IF (LOGIC) THEN T = SIGN(ONE,TT) CALL EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) ELSE CALL DXDT(TT,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = 2.0*(EIG*RX+QX) IF (C.LT.0.0) THEN VALUE = 0.0 IF (P0.GT.0.0) VALUE = 0.5*PI ELSE HH = ABS(XEND-X) D = 2.0*HH/PX CD = C*D*HH IF (P0.GT.0.0) THEN VALUE = C*HH IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) VALUE = VALUE + 0.5*PI ELSE VALUE = D IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) END IF END IF END IF END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE DXDT(T,DT,X) REAL T,DT,X C ********** C C This subroutine transforms coordinates from T on (-1,1) to C X on (A,B) in the solution of a Sturm-Liouville problem. C It is called from subroutines SLEIGN, ALFBET, F, and EXTRAP. C C ********** C .. Scalars in Common .. INTEGER INTAB REAL A,B,C1,C2 C .. C .. Local Scalars .. REAL U C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /DATADT/A,B,C1,C2,INTAB C .. U = C1*T + C2 GO TO (10,20,30,40), INTAB 10 CONTINUE DT = C1*0.5*(B-A) X = 0.5*((B+A)+(B-A)*U) RETURN 20 CONTINUE DT = C1*2.0/(1.0-U)**2 X = A + (1.0+U)/(1.0-U) RETURN 30 CONTINUE DT = C1*2.0/(1.0+U)**2 X = B - (1.0-U)/(1.0+U) RETURN 40 CONTINUE DT = C1/(1.0-ABS(U))**2 X = U/(1.0-ABS(U)) RETURN END C ---------------------------------------------------------------------- SUBROUTINE ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) INTEGER MF,ML,IMAX,IMIN,IA,IB LOGICAL LIMIT REAL ELIM,EMAX,EMIN,PIN,EEE,EIG,EL,WL,DEDW REAL QS(ML),RS(ML),DS(ML),PS(ML) C ********** C C This subroutine generates an initial guess for a specified C eigenvalue of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is C called from SLEIGN when no initial guess is provided by the user. C C The method used is to approximately solve the equation for EIG C C Integral (sqrt((eig*r+q)/p)) = numeig*pi C C where the integral is taken over those X in (A,B) for which C C (eig*r+q)/p .gt. 0 C C and NUMEIG is the index of the desired eigenvalue (PIN=NUMEIG*pi). C C Subprograms called C C sleign-supplied ... estpac C C ********** C .. Scalars in Common .. INTEGER INTAB REAL A,B,C1,C2 C .. C .. Local Scalars .. INTEGER IE,IP LOGICAL LOGIC REAL BALLPK,EU,FNEW,FOLD,SUM,TEMP,U,WU C .. C .. External Subroutines .. EXTERNAL ESTPAC C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. C .. Common blocks .. COMMON /DATADT/A,B,C1,C2,INTAB C .. EEE = MIN(ELIM,EMAX) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C Choose bounds for EIG and associate function (integral) values. C EL = EMIN WL = 0.0 EU = EEE WU = SUM IF (LIMIT .AND. WU.LT.PIN) THEN EIG = ELIM ELSE IF (U.EQ.0.0) THEN EL = EMAX EEE = EMAX + 1.0 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM END IF 10 CONTINUE IF (WU.LE.PIN) THEN C C Increase trial value if integral is still too small. C EL = EU WL = WU EEE = EU + ((PIN-WU+3.0)/U)**2 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM GO TO 10 END IF C C EIG is bracketed. Now try to reduce the size of the bracket C by searching among the saved values of -QS()/RS(). C 20 CONTINUE IF (ABS(IMAX-IMIN).GE.2 .AND. EU.LE.EMAX) THEN IE = (IMAX+IMIN)/2 IF (RS(IE).NE.0.0) THEN EEE = -QS(IE)/RS(IE) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.GT.PIN) THEN IMAX = IE WU = SUM EU = EEE ELSE IMIN = IE WL = SUM EL = EEE END IF GO TO 20 END IF END IF C C Improve approximation for EIG using bisection or secant method. C Substitute 'ballpark' estimate if approximation grows too large. C DEDW = (EU-EL)/(WU-WL) FOLD = 0.0 IF (INTAB.EQ.1) BALLPK = (PIN/(A-B))**2 LOGIC = .TRUE. 30 CONTINUE IF (LOGIC) THEN LOGIC = (WL.LT.(PIN-1.0) .OR. WU.GT.(PIN+1.0)) EEE = EL + DEDW*(PIN-WL) FNEW = MIN(PIN-WL,WU-PIN) IF (FNEW.GT.0.4*FOLD .OR. FNEW.LE.1.0) 1 EEE = 0.5*(EL+EU) IF (INTAB.EQ.1 .AND. ABS(EEE).GT.1.0E+3*BALLPK) THEN EIG = BALLPK RETURN ELSE IF (INTAB.NE.1 .AND. ABS(EEE).GT.1.0E+6) THEN EIG = 1.0 RETURN ELSE FOLD = FNEW C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.LT.PIN) THEN EL = EEE WL = SUM ELSE EU = EEE WU = SUM END IF DEDW = (EU-EL)/(WU-WL) GO TO 30 END IF END IF END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,ZAV) INTEGER MF,ML,IA,IB,IP REAL EEE,PIN,SUM,U,ZAV REAL QS(ML),RS(ML),DS(ML),PS(ML) C ********** C C This subroutine estimates the change in 'phase angle' in the C eigenvalue determination of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is C called from SLEIGN, and also from ESTEIG when no initial guess C is provided by the user. C C The subroutine approximates the (trapezoidal rule) integral of C C sqrt((eig*r+q)/p) C C where the integral is taken over those X in (A,B) for which C C (eig*r+q)/p .gt. 0 C C ********** C .. Local Scalars .. INTEGER J REAL PSUM,RT,RTSAV,V,W,WSAV C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN,SQRT C .. IA = MF IB = 80 IP = MF C C SUM accumulates the integral approximation. U measures the total C length of subintervals where (EIG*R+Q)/P .gt. 0.0. ZAV is the C average value of sqrt((EIG*R+Q)*P) over those subintervals. C SUM = 0.0 U = 0.0 ZAV = 0.0 WSAV = EEE*RS(MF) + QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 10 J=MF+1,ML W = EEE*RS(J) + QS(J) IF (W.GT.0.0) THEN IF (J.GE.80) IB = J U = U + DS(J-1) RT = SQRT(W) ELSE RT = 0.0 IF (U.EQ.0.0 .AND. RTSAV.EQ.0.0 .AND. IA.LE.19) IA = IA + 1 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(J-1)*V IF (PSUM.LT.(PIN-SUM)) IP = J SUM = SUM + PSUM IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(J)+PS(J-1)) 10 CONTINUE SUM = 0.5*SUM ZAV = 0.25*ZAV RETURN END C ---------------------------------------------------------------------- SUBROUTINE EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) INTEGER IFLAG REAL T,TT,EIG,VALUE,DERIV C ********** C C This subroutine is called from ALFBET in determining boundary C values at a singular endpoint of the interval for a C Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C C EXTRAP, which in turn calls INTPOL, extrapolates the function C C arctan(1.0/sqrt(-p*(eig*r+q))) C C from its values for T within (-1,1) to an endpoint. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt,intpol C C ********** C .. Scalars in Common .. REAL Z C .. C .. Local Scalars .. INTEGER KGOOD REAL ANS,CTN,ERROR,PROD,PX,QX,RX,T1,TEMP,X C .. C .. Local Arrays .. REAL FN1(5),XN(5) C .. C .. External Functions .. REAL P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT,INTPOL C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SQRT,TAN C .. C .. Common blocks .. COMMON /ZEE/Z C .. IFLAG = 1 KGOOD = 0 T1 = TT 10 CONTINUE CALL DXDT(T1,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z PROD = -PX*(EIG*RX+QX) IF (PROD.LE.0.0) THEN T1 = 0.5*(T1+T) IF ((1.0+(T1-T)**2).GT.1.0) GO TO 10 IFLAG = 5 RETURN ELSE KGOOD = KGOOD + 1 XN(KGOOD) = T1 FN1(KGOOD) = ATAN(1.0/SQRT(PROD)) T1 = 0.5*(T+T1) IF (KGOOD.LT.5) GO TO 10 END IF T1 = 0.01 CALL INTPOL(5,XN,FN1,T,T1,3,ANS,ERROR) VALUE = ABS(ANS) CTN = 1.0/TAN(VALUE) DERIV = 0.5*PX*RX/CTN/(1.0+CTN**2) TT = XN(1) RETURN END C ---------------------------------------------------------------------- SUBROUTINE F(T,Y,YP) REAL T REAL Y(2),YP(3) C ********** C C This subroutine evaluates the derivative functions for use with C integrator GERK in solving a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt C C ********** C .. Scalars in Common .. LOGICAL EIGF REAL EIG,Z C .. C .. Local Scalars .. REAL C,C2,DT,QX,RX,S,S2,V,W,X,XP,ZP C .. C .. External Functions .. REAL P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,SIN C .. C .. Common blocks .. COMMON /DATAF/EIG,EIGF COMMON /ZEE/Z C .. DT = 1.0 X = T IF (.NOT.EIGF) CALL DXDT(T,DT,X) ZP = ABS(Z) XP = ZP/P(X) QX = Q(X)/ZP RX = R(X)/ZP V = EIG*RX + QX S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = DT*(XP*C2+V*S2) W = (XP-V)*S*C IF (Z.LT.0.0) RX = ABS(RX) YP(2) = DT*(-2.0*W*Y(2)+RX*S2) YP(3) = DT*W RETURN END C ---------------------------------------------------------------------- SUBROUTINE G(X,Y,YP) REAL X REAL Y(1),YP(1) C ********** C C This subroutine evaluates the derivative function for use with C integrator GERK in solving a differential equation in the form C C (p(x)*y'(x))' + q(x)*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P and Q. C C Subprograms called C C user-supplied ..... p,q C C ********** C .. Scalars in Common .. REAL Z C .. C .. Local Scalars .. REAL C,C2,QX,S,S2,XP C .. C .. External Functions .. REAL P,Q EXTERNAL P,Q C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /ZEE/Z C .. XP = Z/P(X) QX = Q(X)/Z S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = XP*C2+QX*S2 RETURN END C ---------------------------------------------------------------------- SUBROUTINE INTPOL(N,XN,FN,X,ABSERR,MAXDEG,ANS,ERROR) INTEGER N,MAXDEG REAL X,ABSERR,ANS,ERROR REAL XN(N),FN(N) C ********** C C This subroutine forms an interpolating polynomial for data pairs. C It is called from EXTRAP in solving a Sturm-Liouville problem. C C ********** C .. Local Scalars .. INTEGER I,I1,II,IJ,IK,IKM1,J,K,L,LIMIT REAL PROD C .. C .. Local Arrays .. INTEGER INDEX(10) REAL V(10,10) C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. L = MIN(MAXDEG,N-2) + 2 LIMIT = MIN(L,N-1) DO 10 I = 1,N V(I,1) = ABS(XN(I)-X) INDEX(I) = I 10 CONTINUE DO 30 I=1,LIMIT DO 20 J=I+1,N II = INDEX(I) IJ = INDEX(J) IF (V(II,1).GT.V(IJ,1)) THEN INDEX(I) = IJ INDEX(J) = II END IF 20 CONTINUE 30 CONTINUE PROD = 1.0 I1 = INDEX(1) ANS = FN(I1) V(1,1) = FN(I1) DO 50 K=2,L IK = INDEX(K) V(K,1) = FN(IK) DO 40 I=1,K-1 II = INDEX(I) V(K,I+1) = (V(I,I)-V(K,I))/(XN(II)-XN(IK)) 40 CONTINUE IKM1 = INDEX(K-1) PROD = (X-XN(IKM1))*PROD ERROR = PROD*V(K,K) IF(ABS(ERROR).LE.ABSERR) RETURN ANS = ANS + ERROR 50 CONTINUE ANS = ANS - ERROR RETURN END C ---------------------------------------------------------------------- SUBROUTINE ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) INTEGER JPAIRS,JSUM REAL A,B,A1,A2,B1,B2 REAL PAIRS(2*JPAIRS) C ********** C C This subroutine counts the zeros, over specified subintervals, of C the solutions of a second order differential equation in the form C C (p(x)*y'(x))' + q(x)*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P and Q. This count in C turn corresponds to the number of zeros, in the interior of (A,B), C of the first eigenfunction of the related Sturm-Liouville problem C whose (semidefinite) weight function vanishes identically in the C subintervals. The problem is restricted to be regular. C C The applicable initial condition depends upon three cases. C C Case 1 -- On a subinterval with left endpoint A, C A1*Y(A) + A2*Y'(A)*P(A) = 0. C C Case 2 -- On a subinterval with right endpoint B, C B1*Y(B) + B2*Y'(B)*P(B) = 0. C C Case 3 -- On a subinterval with neither A nor B as endpoint, C Y(XAA) = 0, where XAA is the left endpoint. C C The SUBROUTINE statement is C C SUBROUTINE zcount(a,b,a1,a2,b1,b2,jpairs,pairs,jsum) C C where C C A and B are input variables defining the full interval. C A must be less than B. C C A1 and A2 are input variables set to prescribe the initial C condition at A (Case 1). C C B1 and B2 are input variables set to prescribe the initial C condition at B (Case 2). C C JPAIRS is an integer input variable set to the number of C specified subintervals of (A,B). C C PAIRS is an input array of length 2*JPAIRS whose successive C ordered element pairs specify the subintervals. C C JSUM is an integer output variable set to the total zero count. C C Subprograms called C C sleign-supplied ... epslon,g,gerk C C user-supplied ..... p,q C C This version dated 5/18/89. C Burton S. Garbow, Argonne National Laboratory C C ********** C .. Scalars in Common .. REAL Z C .. C .. Local Scalars .. INTEGER I,J,LFLAG,MF,ML REAL EPS,EPSMIN,H,ONE,PI,PSUM,PX,QX,RT,RTSAV, 1 T,U,V,W,WSAV,X,X50,XAA,XBB,XSAV,ZAV C .. C .. Local Arrays .. INTEGER IWORK(5) REAL DS(98),GERROR(1),PS(99),QS(99),WORK(11),Y(1) C .. C .. External Functions .. REAL EPSLON,P,Q EXTERNAL EPSLON,P,Q C .. C .. External Subroutines .. EXTERNAL G,GERK C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,INT,SIGN,SQRT C .. C .. Common blocks .. COMMON /ZEE/Z C .. C Set constants EPSMIN, the computer unit roundoff error, and PI. C (Variable ONE set to 1.0 eases precision conversion.) C ONE = 1.0 EPSMIN = EPSLON(ONE) PI = 4.0*ATAN(ONE) C C Set relative and absolute error tolerances for GERK. C EPS = SQRT(EPSMIN) C JSUM = 0 DO 70 J = 1,JPAIRS XAA = PAIRS(2*J-1) XBB = PAIRS(2*J) C DO (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) C C Evaluate P, Q to obtain preliminary information about the C differential equation. C C DO (SAMPLE-COEFFICIENTS) X50 = 0.5*((XBB+XAA)+(XBB-XAA)*EPSMIN) PX = P(X50) QX = Q(X50) PS(50) = PX QS(50) = QX/PX C C MF and ML are the least and greatest index values, respectively. C XSAV = X50 H = 0.9/40.0 DO 10 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I) = XSAV - X XSAV = X MF = I 10 CONTINUE XSAV = X50 DO 20 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I-1) = X - XSAV XSAV = X ML = I - 1 20 CONTINUE C END (SAMPLE-COEFFICIENTS) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) C C U measures the total length of subintervals where Q/P .gt. 0.0. C ZAV is the average value of sqrt(Q*P) over those subintervals. C U = 0.0 ZAV = 0.0 WSAV = QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 30 I=MF+1,ML W = QS(I) IF (W.GT.0.0) THEN U = U + DS(I-1) RT = SQRT(W) ELSE RT = 0.0 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. 1 W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(I-1)*V IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(I)+PS(I-1)) 30 CONTINUE ZAV = 0.25*ZAV C END (ESTIMATE-PHASE-ANGLE-CHANGE) Z = 1.0 IF (U.GT.0.0) Z = ZAV/U C END (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) LFLAG = 1 IF (XAA.EQ.A) THEN C C Case 1 ---------- C Y(1) = PI/2.0 IF (A1.NE.0.0) Y(1) = ATAN(-Z*A2/A1) IF (Y(1).LT.0.0) Y(1) = Y(1) + PI 40 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 40 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) ELSE IF (XBB.EQ.B) THEN C C Case 2 ---------- C Y(1) = PI/2.0 IF (B1.NE.0.0) Y(1) = ATAN(-Z*B2/B1) IF (Y(1).GT.0.0) Y(1) = Y(1) - PI 50 CONTINUE CALL GERK(G,1,Y,XBB,XAA,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 50 JSUM = JSUM - INT((Y(1)-ABS(EPS))/PI) ELSE C C Case 3 ---------- C Y(1) = 0.0 60 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 60 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) END IF 70 CONTINUE RETURN END SUBROUTINE GERK(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK, IWORK) C C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C C WRITTEN BY H.A.WATTS AND L.F.SHAMPINE C SANDIA LABORATORIES C C GERK IS DESIGNED TO SOLVE SYSTEMS OF DIFFERENTIAL EQUATIONS C WHEN IT IS IMPORTANT TO HAVE A READILY AVAILABLE GLOBAL ERROR C ESTIMATE. PARALLEL INTEGRATION IS PERFORMED TO YIELD TWO C SOLUTIONS ON DIFFERENT MESH SPACINGS AND GLOBAL EXTRAPOLATION C IS APPLIED TO PROVIDE AN ESTIMATE OF THE GLOBAL ERROR IN THE C MORE ACCURATE SOLUTION. C C FOR IBM SYSTEM 360 AND 370 AND OTHER MACHINES OF SIMILAR C ARITHMETIC CHARACTERISTICS, THIS CODE SHOULD BE CONVERTED TO C DOUBLE PRECISION. C C******************************************************************* C ABSTRACT C******************************************************************* C C SUBROUTINE GERK INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN)) C WHERE THE Y(I) ARE GIVEN AT T . C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TOUT C BUT IT CAN BE USED AS A ONE-STEP INTEGRATOR TO ADVANCE THE C SOLUTION A SINGLE STEP IN THE DIRECTION OF TOUT. ON RETURN, AN C ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T IS PROVIDED C AND THE PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE C INTEGRATION. THE USER HAS ONLY TO CALL GERK AGAIN (AND PERHAPS C DEFINE A NEW VALUE FOR TOUT). ACTUALLY, GERK IS MERELY AN C INTERFACING ROUTINE WHICH ALLOCATES VIRTUAL STORAGE IN THE C ARRAYS WORK, IWORK AND CALLS SUBROUTINE GERKS FOR THE SOLUTION. C GERKS IN TURN CALLS SUBROUTINE FEHL WHICH COMPUTES AN APPROX- C IMATE SOLUTION OVER ONE STEP. C C GERK USES THE RUNGE-KUTTA-FEHLBERG (4,5) METHOD DESCRIBED C IN THE REFERENCE C E.FEHLBERG , LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITH C STEPSIZE CONTROL , NASA TR R-315 C C C THE PARAMETERS REPRESENT- C F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES C YP(I)=DY(I)/DT C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- SOLUTION VECTOR AT T C T -- INDEPENDENT VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR C LOCAL ERROR TEST. AT EACH STEP THE CODE REQUIRES THAT C ABS(LOCAL ERROR) .LE. RELERR*ABS(Y) + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION C VECTORS. C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION C GERROR(*) -- VECTOR WHICH ESTIMATES THE GLOBAL ERROR AT T. C THAT IS, GERROR(I) APPROXIMATES Y(I)-TRUE C SOLUTION(I). C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO GERK WHICH C IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED C AT LEAST 3+8*NEQN C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL C TO GERK WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST C BE DIMENSIONED AT LEAST 5 C C C******************************************************************* C FIRST CALL TO GERK C******************************************************************* C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE C ARRAYS IN THE CALL LIST - Y(NEQN), WORK(3+8*NEQN), IWORK(5), C DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE F(T,Y,YP) C AND INITIALIZE THE FOLLOWING PARAMETERS- C C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED. (NEQN .GE. 1) C Y(*) -- VECTOR OF INITIAL CONDITIONS C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED. C T=TOUT IS ALLOWED ON THE FIRST CALL ONLY,IN WHICH CASE C GERK RETURNS WITH IFLAG=2 IF CONTINUATION IS POSSIBLE. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C WHICH MUST BE NON-NEGATIVE BUT MAY BE CONSTANTS. WE CAN C USUALLY EXPECT THE GLOBAL ERRORS TO BE SOMEWHAT SMALLER C THAN THE REQUESTED LOCAL ERROR TOLERANCES. TO AVOID C LIMITING PRECISION DIFFICULTIES THE CODE ALWAYS USES C THE LARGER OF RELERR AND AN INTERNAL RELATIVE ERROR C PARAMETER WHICH IS MACHINE DEPENDENT. C IFLAG -- +1,-1 INDICATOR TO INITIALIZE THE CODE FOR EACH NEW C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG= C -1 ONLY WHEN ONE-STEP INTEGRATOR CONTROL IS ESSENTIAL. C IN THIS CASE, GERK ATTEMPTS TO ADVANCE THE SOLUTION A C SINGLE STEP IN THE DIRECTION OF TOUT EACH TIME IT IS C CALLED. SINCE THIS MODE OF OPERATION RESULTS IN EXTRA C COMPUTING OVERHEAD, IT SHOULD BE AVOIDED UNLESS NEEDED. C C C******************************************************************* C OUTPUT FROM GERK C******************************************************************* C C Y(*) -- SOLUTION AT T C T -- LAST POINT REACHED IN INTEGRATION. C IFLAG = 2 -- INTEGRATION REACHED TOUT. INDICATES SUCCESSFUL C RETURN AND IS THE NORMAL MODE FOR CONTINUING C INTEGRATION. C =-2 -- A SINGLE SUCCESSFUL STEP IN THE DIRECTION OF C TOUT HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING C INTEGRATION ONE STEP AT A TIME. C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN C 9000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS C IS APPROXIMATELY 500 STEPS. C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION C VANISHED MAKING A PURE RELATIVE ERROR TEST C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE. C USING THE ONE-STEP INTEGRATION MODE FOR ONE STEP C IS A GOOD WAY TO PROCEED. C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE C ATTEMPTED. C = 6 -- GERK IS BEING USED INEFFICIENTLY IN SOLVING C THIS PROBLEM. TOO MUCH OUTPUT IS RESTRICTING THE C NATURAL STEPSIZE CHOICE. USE THE ONE-STEP C INTEGRATOR MODE. C = 7 -- INVALID INPUT PARAMETERS C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS C SATISFIED - NEQN .LE. 0 C T=TOUT AND IFLAG .NE. +1 OR -1 C RELERR OR ABSERR .LT. 0. C IFLAG .EQ. 0 OR .LT. -2 OR .GT. 7 C GERROR(*) -- ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO C INTEREST TO THE USER BUT NECESSARY FOR SUBSEQUENT C CALLS. WORK(1),...,WORK(NEQN) CONTAIN THE FIRST C DERIVATIVES OF THE SOLUTION VECTOR Y AT T. C WORK(NEQN+1) CONTAINS THE STEPSIZE H TO BE C ATTEMPTED ON THE NEXT STEP. IWORK(1) CONTAINS C THE DERIVATIVE EVALUATION COUNTER. C C C******************************************************************* C SUBSEQUENT CALLS TO GERK C******************************************************************* C C SUBROUTINE GERK RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT, THE USER NEED C ONLY DEFINE A NEW TOUT AND CALL GERK AGAIN. IN THE ONE-STEP C INTEGRATOR MODE (IFLAG=-2) THE USER MUST KEEP IN MIND THAT EACH C STEP TAKEN IS IN THE DIRECTION OF THE CURRENT TOUT. UPON C REACHING TOUT (INDICATED BY CHANGING IFLAG TO 2), THE USER MUST C THEN DEFINE A NEW TOUT AND RESET IFLAG TO -2 TO CONTINUE IN THE C ONE-STEP INTEGRATOR MODE. C C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS C TO CONTINUE (IFLAG=3 CASE), HE JUST CALLS GERK AGAIN. THE C FUNCTION COUNTER IS THEN RESET TO 0 AND ANOTHER 9000 FUNCTION C EVALUATIONS ARE ALLOWED. C C HOWEVER, IN THE CASE IFLAG=4, THE USER MUST FIRST ALTER THE C ERROR CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE C INTEGRATION CAN PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED. C C ALSO, IN THE CASE IFLAG=5, IT IS NECESSARY FOR THE USER TO C RESET IFLAG TO 2 (OR -2 WHEN THE ONE-STEP INTEGRATION MODE IS C BEING USED) AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH C BEFORE THE INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, C EXECUTION WILL BE TERMINATED. THE OCCURRENCE OF IFLAG=5 C INDICATES A TROUBLE SPOT (SOLUTION IS CHANGING RAPIDLY, C SINGULARITY MAY BE PRESENT) AND IT OFTEN IS INADVISABLE TO C CONTINUE. C C IF IFLAG=6 IS ENCOUNTERED, THE USER SHOULD USE THE ONE-STEP C INTEGRATION MODE WITH THE STEPSIZE DETERMINED BY THE CODE. IF C THE USER INSISTS UPON CONTINUING THE INTEGRATION WITH GERK IN C THE INTERVAL MODE, HE MUST RESET IFLAG TO 2 BEFORE CALLING GERK C AGAIN. OTHERWISE,EXECUTION WILL BE TERMINATED. C C IF IFLAG=7 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS C THE INVALID INPUT PARAMETERS ARE CORRECTED. C C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN C INFORMATION REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, C WORK AND IWORK SHOULD NOT BE ALTERED. C C******************************************************************* C C .. SCALAR ARGUMENTS .. INTEGER IFLAG,NEQN REAL ABSERR,RELERR,T,TOUT C .. C .. ARRAY ARGUMENTS .. INTEGER IWORK(5) REAL GERROR(NEQN),WORK(3+8*NEQN),Y(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K1,K1M,K2,K3,K4,K5,K6,K7,K8 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GERKS C .. C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY K1M = NEQN + 1 K1 = K1M + 1 K2 = K1 + NEQN K3 = K2 + NEQN K4 = K3 + NEQN K5 = K4 + NEQN K6 = K5 + NEQN K7 = K6 + NEQN K8 = K7 + NEQN C ******************************************************************* C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER, C HE MUST USE GERKS DIRECTLY. C ******************************************************************* CALL GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK(1), WORK(K1M), WORK(K1), WORK(K2), WORK(K3), * WORK(K4), WORK(K5), WORK(K6), WORK(K7), WORK(K8), * WORK(K8+1), IWORK(1), IWORK(2), IWORK(3), IWORK(4), IWORK(5)) RETURN END SUBROUTINE GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, YP, H, F1, F2, F3, F4, F5, YG, YGP, SAVRE, SAVAE, * NFE, KOP, INIT, JFLAG, KFLAG) C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C ******************************************************************* C GERKS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS AS DESCRIBED IN THE COMMENTS FOR GERK. THE ARRAYS C YP,F1,F2,F3,F4,F5,YG AND YGP (OF DIMENSION AT LEAST NEQN) AND C THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE C USED INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO C ELIMINATE LOCAL RETENTION OF VARIABLES BETWEEN CALLS. C ACCORDINGLY, THEY SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE C INTEREST ARE C YP - DERIVATIVE OF SOLUTION VECTOR AT T C H - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP C NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION C EVALUATIONS. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER IFLAG,INIT,JFLAG,KFLAG,KOP,NEQN,NFE REAL ABSERR,H,RELERR,SAVAE,SAVRE,T,TOUT C .. C .. ARRAY ARGUMENTS .. REAL F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 GERROR(NEQN),Y(NEQN),YG(NEQN),YGP(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K,MAXNFE,MFLAG LOGICAL HFAILD,OUTPUT REAL A,AE,DT,EE,EEOET,ESTTOL,ET,HH,HMIN,ONE,REMIN,RER, 1 S,SCALE,TOL,TOLN,TS,U,U26,YPK C .. C .. EXTERNAL FUNCTIONS .. REAL EPSLON EXTERNAL EPSLON C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FEHL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MAX,MIN,SIGN C .. C ******************************************************************* C REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE C INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL C GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING C PRECISION ON COMPUTERS WITH LONG WORDLENGTHS. DATA REMIN /3.E-11/ C ******************************************************************* C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE. C AS SET,THIS CORRESPONDS TO ABOUT 500 STEPS. DATA MAXNFE /9000/ C ******************************************************************* C U - THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE C VALUE REPRESENTABLE IN THE MACHINE SUCH THAT 1.+ U .GT. 1. C (VARIABLE ONE SET TO 1.0 EASES PRECISION CONVERSION.) C ONE = 1.0 U = EPSLON(ONE) C ******************************************************************* C CHECK INPUT PARAMETERS IF (NEQN.LT.1) GO TO 10 IF ((RELERR.LT.0.) .OR. (ABSERR.LT.0.)) GO TO 10 MFLAG = ABS(IFLAG) IF ((MFLAG.GE.1) .AND. (MFLAG.LE.7)) GO TO 20 C INVALID INPUT 10 IFLAG = 7 RETURN C IS THIS THE FIRST CALL 20 IF (MFLAG.EQ.1) GO TO 70 C CHECK CONTINUATION POSSIBILITIES IF (T.EQ.TOUT) GO TO 10 IF (MFLAG.NE.2) GO TO 30 C IFLAG = +2 OR -2 IF (INIT.EQ.0) GO TO 60 IF (KFLAG.EQ.3) GO TO 50 IF ((KFLAG.EQ.4) .AND. (ABSERR.EQ.0.)) GO TO 40 IF ((KFLAG.EQ.5) .AND. (RELERR.LE.SAVRE) .AND. * (ABSERR.LE.SAVAE)) GO TO 40 GO TO 70 C IFLAG = 3,4,5,6, OR 7 30 IF (IFLAG.EQ.3) GO TO 50 IF ((IFLAG.EQ.4) .AND. (ABSERR.GT.0.)) GO TO 60 C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO C THE INSTRUCTIONS PERTAINING TO IFLAG=4,5,6 OR 7 40 STOP C ******************************************************************* C RESET FUNCTION EVALUATION COUNTER 50 NFE = 0 IF (MFLAG.EQ.2) GO TO 70 C RESET FLAG VALUE FROM PREVIOUS CALL 60 IFLAG = JFLAG C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT C INPUT CHECKING 70 JFLAG = IFLAG KFLAG = 0 C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS SAVRE = RELERR SAVAE = ABSERR C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS C 32U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING C FROM IMPOSSIBLE ACCURACY REQUESTS RER = MAX(RELERR,32.*U+REMIN) U26 = 26.*U DT = TOUT - T IF (MFLAG.EQ.1) GO TO 80 IF (INIT.EQ.0) GO TO 90 GO TO 110 C ******************************************************************* C INITIALIZATION -- C SET INITIALIZATION COMPLETION INDICATOR,INIT C SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP C EVALUATE INITIAL DERIVATIVES C COPY INITIAL VALUES AND DERIVATIVES FOR THE C PARALLEL SOLUTION C SET COUNTER FOR FUNCTION EVALUATIONS,NFE C ESTIMATE STARTING STEPSIZE 80 INIT = 0 KOP = 0 A = T CALL F(A, Y, YP) NFE = 1 IF (T.NE.TOUT) GO TO 90 IFLAG = 2 RETURN 90 INIT = 1 H = ABS(DT) TOLN = 0. DO 100 K=1,NEQN YG(K) = Y(K) YGP(K) = YP(K) TOL = RER*ABS(Y(K)) + ABSERR IF (TOL.LE.0.) GO TO 100 TOLN = TOL YPK = ABS(YP(K)) IF (YPK*H**5.GT.TOL) H = (TOL/YPK)**0.2 100 CONTINUE IF (TOLN.LE.0.) H = 0. H = MAX(H,U26*MAX(ABS(T),ABS(DT))) C ******************************************************************* C SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT 110 H = SIGN(H,DT) C TEST TO SEE IF GERK IS BEING SEVERELY IMPACTED BY TOO MANY C OUTPUT POINTS IF (ABS(H).GT.2.*ABS(DT)) KOP = KOP + 1 IF (KOP.NE.100) GO TO 120 KOP = 0 IFLAG = 6 RETURN 120 IF (ABS(DT).GT.U26*ABS(T)) GO TO 140 C IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN DO 130 K=1,NEQN YG(K) = YG(K) + DT*YGP(K) Y(K) = Y(K) + DT*YP(K) 130 CONTINUE A = TOUT CALL F(A, YG, YGP) CALL F(A, Y, YP) NFE = NFE + 2 GO TO 230 C INITIALIZE OUTPUT POINT INDICATOR 140 OUTPUT = .FALSE. C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION, C SCALE THE ERROR TOLERANCES SCALE = 2./RER AE = SCALE*ABSERR C ******************************************************************* C ******************************************************************* C STEP BY STEP INTEGRATION 150 HFAILD = .FALSE. C SET SMALLEST ALLOWABLE STEPSIZE HMIN = U26*ABS(T) C ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT. C LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE C AND THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE. DT = TOUT - T IF (ABS(DT).GE.2.*ABS(H)) GO TO 170 IF (ABS(DT).GT.ABS(H)) GO TO 160 C THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE C OUTPUT POINT OUTPUT = .TRUE. H = DT GO TO 170 160 H = 0.5*DT C ******************************************************************* C CORE INTEGRATOR FOR TAKING A SINGLE STEP C ******************************************************************* C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW C IN COMPUTING THE ERROR TOLERANCE FUNCTION ET. C TO AVOID PROBLEMS WITH ZERO CROSSINGS, RELATIVE ERROR IS C MEASURED USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION C AT THE BEGINNING AND END OF A STEP. C THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF C SIGNIFICANCE. C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T. C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES. C TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE C STEPSIZE IT ESTIMATES WILL SUCCEED. C AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE C FOR THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE C EFFICIENT ON PROBLEMS HAVING DISCONTINUITIES AND MORE C EFFECTIVE IN GENERAL SINCE LOCAL EXTRAPOLATION IS BEING USED C AND THE ERROR ESTIMATE MAY BE UNRELIABLE OR UNACCEPTABLE WHEN C A STEP FAILS. C ******************************************************************* C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS. C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H 170 IF (NFE.LE.MAXNFE) GO TO 180 C TOO MUCH WORK IFLAG = 3 KFLAG = 3 RETURN C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H 180 CALL FEHL(F, NEQN, YG, T, H, YGP, F1, F2, F3, F4, F5, F1) NFE = NFE + 5 C COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR C ESTIMATES AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE C ERROR IS MEASURED WITH RESPECT TO THE AVERAGE MAGNITUDES OF THE C OF THE SOLUTION AT THE BEGINNING AND END OF THE STEP. EEOET = 0. DO 200 K=1,NEQN ET = ABS(YG(K)) + ABS(F1(K)) + AE IF (ET.GT.0.) GO TO 190 C INAPPROPRIATE ERROR TOLERANCE IFLAG = 4 KFLAG = 4 RETURN 190 EE = ABS((-2090.*YGP(K)+(21970.*F3(K)-15048.*F4(K))) * +(22528.*F2(K)-27360.*F5(K))) EEOET = MAX(EEOET,EE/ET) 200 CONTINUE ESTTOL = ABS(H)*EEOET*SCALE/752400. IF (ESTTOL.LE.1.) GO TO 210 C UNSUCCESSFUL STEP C REDUCE THE STEPSIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR OF 1/10 HFAILD = .TRUE. OUTPUT = .FALSE. S = 0.1 IF (ESTTOL.LT.59049.) S = 0.9/ESTTOL**0.2 H = S*H IF (ABS(H).GT.HMIN) GO TO 170 C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE IFLAG = 5 KFLAG = 5 RETURN C SUCCESSFUL STEP C STORE ONE-STEP SOLUTION YG AT T+H C AND EVALUATE DERIVATIVES THERE 210 TS = T T = T + H DO 220 K=1,NEQN YG(K) = F1(K) 220 CONTINUE A = T CALL F(A, YG, YGP) NFE = NFE + 1 C NOW ADVANCE THE Y SOLUTION OVER TWO STEPS OF C LENGTH H/2 AND EVALUATE DERIVATIVES THERE HH = 0.5*H CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) TS = TS + HH A = TS CALL F(A, Y, YP) CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) A = T CALL F(A, Y, YP) NFE = NFE + 12 C CHOOSE NEXT STEPSIZE C THE INCREASE IS LIMITED TO A FACTOR OF 5 C IF STEP FAILURE HAS JUST OCCURRED, NEXT C STEPSIZE IS NOT ALLOWED TO INCREASE S = 5. IF (ESTTOL.GT.1.889568E-4) S = 0.9/ESTTOL**0.2 IF (HFAILD) S = MIN(S,ONE) H = SIGN(MAX(S*ABS(H),HMIN),H) C ******************************************************************* C END OF CORE INTEGRATOR C ******************************************************************* C SHOULD WE TAKE ANOTHER STEP IF (OUTPUT) GO TO 230 IF (IFLAG.GT.0) GO TO 150 C ******************************************************************* C ******************************************************************* C INTEGRATION SUCCESSFULLY COMPLETED C ONE-STEP MODE IFLAG = -2 GO TO 240 C INTERVAL MODE 230 T = TOUT IFLAG = 2 240 DO 250 K=1,NEQN GERROR(K) = (YG(K)-Y(K))/31. 250 CONTINUE RETURN END SUBROUTINE FEHL(F, NEQN, Y, T, H, YP, F1, F2, F3, F4, F5, S) C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD C ******************************************************************* C FEHL INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT=F(T,Y(1),---,Y(NEQN)) C WHERE THE INITIAL VALUES Y(I) AND THE INITIAL DERIVATIVES C YP(I) ARE SPECIFIED AT THE STARTING POINT T. FEHL ADVANCES C THE SOLUTION OVER THE FIXED STEP H AND RETURNS C THE FIFTH ORDER (SIXTH ORDER ACCURATE LOCALLY) SOLUTION C APPROXIMATION AT T+H IN ARRAY S(I). C F1,---,F5 ARE ARRAYS OF DIMENSION NEQN WHICH ARE NEEDED C FOR INTERNAL STORAGE. C THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF SIGNIFICANCE. C FEHL SHOULD BE CALLED WITH AN H NOT SMALLER THAN 13 UNITS OF C ROUNDOFF IN T SO THAT THE VARIOUS INDEPENDENT ARGUMENTS CAN BE C DISTINGUISHED. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER NEQN REAL H,T C .. C .. ARRAY ARGUMENTS .. REAL F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 S(NEQN),Y(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K REAL CH C .. CH = 0.25*H DO 10 K=1,NEQN F5(K) = Y(K) + CH*YP(K) 10 CONTINUE CALL F(T+0.25*H, F5, F1) CH = 0.09375*H DO 20 K=1,NEQN F5(K) = Y(K) + CH*(YP(K)+3.*F1(K)) 20 CONTINUE CALL F(T+0.375*H, F5, F2) CH = H/2197. DO 30 K=1,NEQN F5(K) = Y(K) + CH*(1932.*YP(K)+(7296.*F2(K)-7200.*F1(K))) 30 CONTINUE CALL F(T+12./13.*H, F5, F3) CH = H/4104. DO 40 K=1,NEQN F5(K) = Y(K) + CH*((8341.*YP(K)-845.*F3(K))+(29440.*F2(K) * -32832.*F1(K))) 40 CONTINUE CALL F(T+H, F5, F4) CH = H/20520. DO 50 K=1,NEQN F1(K) = Y(K) + CH*((-6080.*YP(K)+(9295.*F3(K)-5643.*F4(K))) * +(41040.*F1(K)-28352.*F2(K))) 50 CONTINUE CALL F(T+0.5*H, F1, F5) C COMPUTE APPROXIMATE SOLUTION AT T+H CH = H/7618050. DO 60 K=1,NEQN S(K) = Y(K) + CH*((902880.*YP(K)+(3855735.*F3(K)-1371249.* * F4(K)))+(3953664.*F2(K)+277020.*F5(K))) 60 CONTINUE RETURN END REAL FUNCTION EPSLON (X) REAL X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C REAL A,B,C,EPS,FOUR,THREE C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C FOUR = 4.0 THREE = 3.0 A = FOUR/THREE 10 B = A - 1.0 C = B + B + B EPS = ABS(C-1.0) IF (EPS .EQ. 0.0) GO TO 10 EPSLON = EPS*ABS(X) RETURN END PROGRAM MAIN C C This sample driver calls ZCOUNT and SLEIGN to find the smallest C eigenvalue, EIG(1), of the semidefinite Sturm-Liouville problem C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (-pi,pi) C C where P(x) = 1, Q(x) = 4, and R(x) = 1 for x .le. -PI/2 or C x .gt. PI/2, and R(x) = 0 otherwise. The correct answer is C EIG(1) = -3, more generally, EIG(n) = n**2 - 4. C INTEGER I,IFLAG,INTAB,ISLFUN,JPAIRS,JSUM,NUMEIG,NWRITE DOUBLE PRECISION A,A1,A2,B,B1,B2,EIG,P0ATA,P0ATB,QFATA,QFATB,TOL DOUBLE PRECISION PAIRS(2),SLFUN(9) EXTERNAL PARAMS,SLEIGN,ZCOUNT C DATA NWRITE/6/ C CALL PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) C C Count number of zeros in (A,B) of eigenfunction associated C with smallest eigenvalue. C CALL ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) WRITE (NWRITE,10) ' zero count =',JSUM 10 FORMAT (A,I6) C C Calculate smallest eigenvalue and associated eigenfunction data. C NUMEIG = 1 + JSUM EIG = 0.0 TOL = 1.0D-6 ISLFUN = 0 CALL SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) WRITE (NWRITE,20) ' iflag,numeig,eig,tol = ',IFLAG,NUMEIG,EIG,TOL 20 FORMAT (A,2I6,F12.5,1PE12.2) WRITE (NWRITE,30) ' slfun(1-9) =',(SLFUN(I),I=1,9) 30 FORMAT (A,F12.5/(13X,3F12.5)) STOP END C SUBROUTINE PARAMS(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 JPAIRS,PAIRS) INTEGER INTAB,JPAIRS DOUBLE PRECISION A,A1,A2,B,B1,B2,ONE,PI,P0ATA,P0ATB,QFATA,QFATB DOUBLE PRECISION PAIRS(*) INTRINSIC ATAN COMMON PI ONE = 1.0 PI = 4.0*ATAN(ONE) A = -PI B = PI INTAB = 1 P0ATA = -1.0 QFATA = 1.0 P0ATB = -1.0 QFATB = 1.0 A1 = 1.0 A2 = 0.0 B1 = 1.0 B2 = 0.0 JPAIRS = 1 PAIRS(1) = -PI/2.0 PAIRS(2) = PI/2.0 RETURN END C FUNCTION P(X) DOUBLE PRECISION P,X P = 1.0 RETURN END C FUNCTION Q(X) DOUBLE PRECISION Q,X Q = 4.0 RETURN END C FUNCTION R(X) DOUBLE PRECISION PI,R,X COMMON PI R = 0.0 IF (X.LE.-PI/2.0 .OR. X.GT.PI/2.0) R = 1.0 RETURN END SUBROUTINE SLEIGN(A,B,INTAB,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2, 1 NUMEIG,EIG,TOL,IFLAG,ISLFUN,SLFUN) INTEGER INTAB,NUMEIG,IFLAG,ISLFUN DOUBLE PRECISION A,B,P0ATA,QFATA,P0ATB,QFATB,A1,A2,B1,B2,EIG,TOL DOUBLE PRECISION SLFUN(9) C ********** C C This subroutine is designed for the calculation of a specified C eigenvalue, EIG, of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C The problem may be either regular or singular. In the C regular case, boundary conditions of the form C C a1*y(a) + a2*p(a)*y'(a) = 0 C b1*y(b) + b2*p(b)*y'(b) = 0 C C are prescribed by specifying the numbers A1, A2, B1, and B2. C The index of the desired eigenvalue is specified in NUMEIG C and its requested accuracy in TOL. Initial data for the C associated eigenfunction are also computed along with values C at selected points, if desired, in array SLFUN. C C The SUBROUTINE statement is C C SUBROUTINE sleign(a,b,intab,p0ata,qfata,p0atb,qfatb,a1,a2,b1,b2, C numeig,eig,tol,iflag,islfun,slfun) C C where C C A and B are input variables defining the interval. If the C interval is finite, A must be less than B. (See INTAB below.) C C INTAB is an integer input variable specifying the nature of the C interval. It can have four values. C C INTAB = 1 - A and B are finite. C INTAB = 2 - A is finite and B is infinite (+). C INTAB = 3 - A is infinite (-) and B is finite. C INTAB = 4 - A is infinite (-) and B is infinite (+). C C If either A or B is infinite, it is classified singular and C its value is ignored. C C P0ATA, QFATA, P0ATB, and QFATB are input variables set to C 1.0 or -1.0 as the following properties of P, Q, and R at C the interval endpoints are true or false, respectively. C C P0ATA - p(a) is zero. (If true, A is singular.) C QFATA - q(a) and r(a) are finite. (If false, A is singular.) C P0ATB - p(b) is zero. (If true, B is singular.) C QFATB - q(b) and r(b) are finite. (If false, B is singular.) C C A1 and A2 are input variables set to prescribe the boundary C condition at A. Their values are ignored if A is singular. C C B1 and B2 are input variables set to prescribe the boundary C condition at B. Their values are ignored if B is singular. C C NUMEIG is an integer variable. On input, it should be set to C the index of the desired eigenvalue (increasing sequence). C On output, it is unchanged unless the problem (apparently) C lacks NUMEIG eigenvalues, in which case it is reset to the C index of the largest eigenvalue. C C EIG is a variable set on input to 0.0 or to an initial guess of C the eigenvalue. If EIG is set to 0.0, SLEIGN will generate C the initial guess. On output, EIG holds the calculated C eigenvalue if IFLAG (see below) signals success. C C TOL is a variable set on input to the desired accuracy of the C eigenvalue. On output, TOL is reset to the accuracy estimated C to have been achieved if IFLAG (see below) signals success. C This accuracy estimate is absolute if EIG is less than one C in magnitude, and relative otherwise. In addition, prefixing C TOL with a negative sign, removed after interrogation, serves C as a flag to request trace output from the calculation. C C IFLAG is an integer output variable set as follows. C C IFLAG = 1 - successful problem solution. C IFLAG = 2 - improper input parameters. C IFLAG = 3 - NUMEIG exceeds actual number of eigenvalues. C IFLAG = 4 - some uncertainty about accuracy estimate TOL. C IFLAG = 5 - convergence too slow, best results returned. C IFLAG = 6 - failure, integrator could not complete. C C ISLFUN is an integer input variable set to the number of C selected eigenfunction values desired. If no values are C desired, set ISLFUN zero or negative. C C SLFUN is an array of length at least 9. On output, the first 9 C locations contain the integration interval and initial data C that completely determine the eigenfunction. C C SLFUN(1) - point where two pieces of eigenfunction Y match. C SLFUN(2) - left endpoint XAA of the (truncated) interval. C SLFUN(3) - value of THETA at XAA. (Y = RHO*sin(THETA)) C SLFUN(4) - value of F at XAA. (RHO = exp(F)) C SLFUN(5) - right endpoint XBB of the (truncated) interval. C SLFUN(6) - value of THETA at XBB. C SLFUN(7) - value of F at XBB. C SLFUN(8) - final value of integration accuracy parameter EPS. C SLFUN(9) - the constant Z in the polar form transformation. C C F(XAA) and F(XBB) are chosen so that the eigenfunction is C continuous in the interval (XAA,XBB) and has weighted (by R) C L2-norm of 1.0 on the interval. If ISLFUN is positive, then C on input the further ISLFUN locations of SLFUN specify the C points, in ascending order, where the eigenfunction values C are desired and on output contain the values themselves. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... alfbet,dxdt,epslon,esteig,estpac,f,gerk C C This version dated 5/18/89. C Paul B. Bailey, Sandia Laboratories, Albuquerque C Burton S. Garbow, Argonne National Laboratory C C ********** C .. Scalars in Common .. INTEGER INTSAV LOGICAL EIGF DOUBLE PRECISION ASAV,BSAV,C1,C2,EIGSAV,Z C .. C .. Local Scalars .. INTEGER I,IA,IB,IMAX,IMIN,IOUT,IP,J,JFLAG,KFLAG,LFLAG,MF,ML, 1 NEIG,NMID LOGICAL AOK,BOK,BRACKT,CHNGAB,CHNGM,CONVRG,FYNYT,FYNYT1, 1 LIMIT,LOGIC,NEWTON,ONEDIG,PRIN,THEGT0,THELT0 DOUBLE PRECISION AA,AAA,ALFA,ASL,ASR,BB,BBB,BETA, 1 C,CHNG,CL,CR,DAV,DE,DEDW,DEN,DERIVL,DERIVR,DIST, 2 DPSIL,DPSIPL,DPSIPR,DPSIR,DT,DTHDA,DTHDAA,DTHDB, 3 DTHDBB,DTHDE,DTHDEA,DTHDEB,DTHETA,DTHOLD,E,EEE, 4 EIGLO,EIGLT,EIGRT,EIGUP,EL,ELIM,EMAX,EMIN,EOLD,EPS,EPSMIN, 5 ER1,ER2,ESTERR,FLO,FMAX,FUP,GMAX,GUESS,H,ONE,PI,PIN, 6 PSIL,PSIPL,PSIPR,PSIR,PX,QAV,QX,RATIO,RAV,RAY,RX, 7 SL,SQL,SQR,SR,T,T1,T2,T3,TAU,TEMP,THRESH,TMAX,TMID,TMP, 8 U,UL,UR,V,WL,X,X50,XAA,XBB,XBC,XMID,XSAV,ZAV C .. C .. Local Arrays .. INTEGER IWORK(5) DOUBLE PRECISION DS(98),ERL(3),ERR(3),PS(99),QS(99),RS(99), 1 WORK(27),YL(3),YR(3) C .. C .. External Functions .. DOUBLE PRECISION EPSLON,P,Q,R EXTERNAL EPSLON,P,Q,R C .. C .. External Subroutines .. EXTERNAL ALFBET,DXDT,ESTEIG,ESTPAC,F,GERK C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,COS,EXP,INT,LOG,MAX,MIN,SIGN,SIN,TAN C .. C .. Common blocks .. COMMON /DATADT/ASAV,BSAV,C1,C2,INTSAV COMMON /DATAF/EIGSAV,EIGF COMMON /ZEE/Z C .. C Set constants EPSMIN, the computer unit roundoff error, and PI. C (Variable ONE set to 1.0 eases precision conversion.) C ONE = 1.0 EPSMIN = EPSLON(ONE) PI = 4.0*ATAN(ONE) C C Set output device number. C IOUT = 6 C C Check input parameters for errors. If errors, return IFLAG=2. C LOGIC = TOL.NE.0.0 .AND. 1.LE.INTAB .AND. INTAB.LE.4 .AND. 1 P0ATA*QFATA*P0ATB*QFATB*NUMEIG.NE.0.0 IF (INTAB.EQ.1) LOGIC = LOGIC .AND. A.LT.B IF (.NOT.LOGIC) THEN IFLAG = 2 RETURN END IF C C Set PRIN = .true. to trigger trace printout of successive steps. C PRIN = .FALSE. IF (TOL.LT.0.0) PRIN = .TRUE. C C Set EPS to the (initial) integration accuracy. C EPS = 0.001 C C AOK (BOK) signals, if true, that endpoint A (B) is not singular. C AOK = INTAB.LT.3 .AND. P0ATA.LT.0.0 .AND. QFATA.GT.0.0 BOK = (INTAB.EQ.1 .OR. INTAB.EQ.3) .AND. 1 P0ATB.LT.0.0 .AND. QFATB.GT.0.0 NEIG = ABS(NUMEIG) - 1 C C Initial C1 and C2, used in the mapping between X and T intervals. C C1 = 1.0 C2 = 0.0 C DO (SAVE-INPUT-DATA) ASAV = A BSAV = B INTSAV = INTAB TAU = ABS(TOL) C END (SAVE-INPUT-DATA) C C Evaluate P, Q, R to obtain preliminary information about the C differential equation. C C DO (SAMPLE-COEFFICIENTS) THRESH = 1.0E+17 10 CONTINUE CALL DXDT(EPSMIN,TEMP,X50) PX = P(X50) QX = Q(X50) RX = R(X50) PS(50) = PX QS(50) = QX/PX RS(50) = RX/PX C C DAV,QAV,RAV are used in special case estimation when NUMEIG = 1,2. C EMIN = min(-Q/R), achieved at X for index value IMIN. C EMAX = max(-Q/R), achieved at X for index value IMAX. C MF and ML are the least and greatest index values, respectively. C DAV = 0.0 QAV = 0.0 RAV = 0.0 XSAV = X50 EMIN = THRESH EMAX = -THRESH IF (RX.NE.0.0) THEN EMIN = -QX/RX EMAX = EMIN IMIN = 50 IMAX = 50 END IF H = 0.9/40.0 DO 20 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I) = XSAV - X DAV = DAV + DS(I) QAV = QAV + DS(I)*(0.5*(QS(I)+QS(I+1))-QAV)/DAV RAV = RAV + DS(I)*(0.5*(RS(I)+RS(I+1))-RAV)/DAV XSAV = X C C Try to avoid overflow by stopping when functions are large near A. C FYNYT = (ABS(RX)+ABS(QX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF MF = I IF (.NOT.FYNYT) GO TO 30 20 CONTINUE 30 CONTINUE AAA = T IF (AOK) AAA = -1.0 XSAV = X50 DO 40 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) CALL DXDT(T,TEMP,X) PX = P(X) QX = Q(X) RX = R(X) PS(I) = PX QS(I) = QX/PX RS(I) = RX/PX DS(I-1) = X - XSAV DAV = DAV + DS(I-1) QAV = QAV + DS(I-1)*(0.5*(QS(I-1)+QS(I))-QAV)/DAV RAV = RAV + DS(I-1)*(0.5*(RS(I-1)+RS(I))-RAV)/DAV XSAV = X C C Try to avoid overflow by stopping when functions are large near B. C FYNYT1 = (ABS(QX)+ABS(RX)+1.0/PX).LE.THRESH IF (RX.NE.0.0) THEN IF (-QX/RX.LT.EMIN) THEN EMIN = -QX/RX IMIN = I END IF IF (-QX/RX.GT.EMAX) THEN EMAX = -QX/RX IMAX = I END IF END IF ML = I - 1 IF (.NOT.FYNYT1) GO TO 50 40 CONTINUE 50 CONTINUE BBB = T IF (BOK) BBB = 1.0 LOGIC = C1.EQ.1.0 .AND. (.NOT.FYNYT .OR. .NOT.FYNYT1) C C Modify (T,X) transformation corresponding to truncated interval. C IF (LOGIC) THEN C1 = 0.5*(BBB-AAA) C2 = 0.5*(AAA+BBB) GO TO 10 END IF C C Estimate upper bound ELIM for EIG such that boundary conditions C can be satisfied. C ELIM = EMAX + 1.0 IF (INTAB.EQ.3 .OR. (P0ATA.GT.0.0 .AND. QFATA.LT.0.0)) THEN IF (-QS(MF)/RS(MF).LE.ELIM) THEN ELIM = -QS(MF)/RS(MF) IMAX = MF END IF END IF IF (INTAB.EQ.2 .OR. (P0ATB.GT.0.0 .AND. QFATB.LT.0.0)) THEN IF (-QS(ML)/RS(ML).LE.ELIM) THEN ELIM = -QS(ML)/RS(ML) IMAX = ML END IF END IF IF (INTAB.EQ.4) THEN ELIM = MIN(ELIM,-QS(MF)/RS(MF),-QS(ML)/RS(ML)) IF (-QS(MF)/RS(MF).EQ.ELIM) IMAX = MF IF (-QS(ML)/RS(ML).EQ.ELIM) IMAX = ML END IF ELIM = ELIM - EPSMIN IF (ELIM.EQ.0.0) ELIM = -EPSMIN LIMIT = ELIM.LE.EMAX C END (SAMPLE-COEFFICIENTS) PIN = (NEIG+1)*PI IF (EIG.EQ.0.0) THEN C DO (ESTIMATE-EIG) CALL ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) C END (ESTIMATE-EIG) END IF GUESS = EIG C DO (SET-INITIAL-INTERVAL-AND-MATCHPOINT) IF (GUESS.NE.0.0) THEN C C Reduce overly large guess for EIG to upper bound if necessary. C IF (LIMIT .AND. EIG.GT.ELIM) EIG = ELIM EEE = EIG C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,TEMP,U,TMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) END IF C C Choose initial interval as large as possible that avoids overflow. C IA and IB are boundary indices for nonnegativity of EIG*R+Q. C AA = -1.0 IF (.NOT.AOK) THEN AA = H*(IA-50) IF (IA.LT.10) AA = -1.0 + 0.1/4.0**(10-IA) END IF BB = 1.0 IF (.NOT.BOK) THEN BB = H*(IB-50) IF (IB.GT.90) BB = 1.0 - 0.1/4.0**(IB-90) END IF AA = AA + 0.6*(AAA-AA) BB = BB + 0.6*(BBB-BB) C C Determine boundary values ALFA and BETA for theta at A and B. C Z = 1.0 CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA C C Take boundary conditions into account in estimation of EIG. C PIN = PIN + BETA - ALFA - PI IF (GUESS.EQ.0.0) EEE = EL + DEDW*(PIN-WL) C C Subroutine ESTPAC must be called again because PIN has changed. C C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,TEMP,U,ZAV) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C Choose the constant Z. C IF (U.GT.0.0) Z = ZAV/U C C Reset boundary values ALFA and BETA. C CALL ALFBET(A,INTAB,AA,A1,A2,EEE,P0ATA,QFATA,AOK, 1 ALFA,KFLAG,DERIVL) CALL ALFBET(B,INTAB,BB,B1,B2,EEE,P0ATB,QFATB,BOK, 1 BETA,JFLAG,DERIVR) IF (.NOT.BOK) BETA = PI - BETA IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' alfa=',ALFA,' beta=',BETA C C Special case formula for estimation of EIG when NUMEIG = 1,2. C IF (U.EQ.0.0 .AND. NEIG.LE.1 .AND. (BETA+NEIG*PI).LT.ALFA) THEN XBC = MAX(-1.0/TAN(ALFA),1.0/TAN(BETA)) EEE = -(XBC*XBC-QAV)/RAV DEDW = XBC*(1.0+XBC*XBC)/RAV END IF C C Choose initial matching point TMID. C TMID = H*(IP-50) IF (TMID.LT.-0.8) TMID = -0.4 IF (TMID.GT.0.8) TMID = 0.4 IF (PRIN) WRITE(IOUT,'(A,E15.7,A,F11.8,A,E15.7)') 1 ' estim=',EEE,' tmid=',TMID,' z=',Z IF (PRIN) WRITE(IOUT,'(A,F11.8,A,F11.8,A,F11.8,A,F11.8)') 1 ' aaa=',AAA,' aa=',AA,' bb=',BB,' bbb=',BBB C END (SET-INITIAL-INTERVAL-AND-MATCHPOINT) C C End preliminary work, begin main task of computing EIG. C C Logical variables have the following meanings if true. C AOK - endpoint A is not singular. C BOK - endpoint B is not singular. C CHNGM - matching point TMID should be changed. C CHNGAB - one or both endpoints should be moved farther out. C BRACKT - EIG has been bracketed. C CONVRG - convergence test for EIG has been successfully passed. C NEWTON - Newton iteration may be employed. C THELT0 - lower bound for EIG has been found. C THEGT0 - upper bound for EIG has been found. C LIMIT - upper bound exists with boundary conditions satisfied. C ONEDIG - most significant digit can be expected to be correct. C EIGF - derivative argument is in original coordinate system. C EIG = EEE EIGF = .FALSE. CHNGM = .FALSE. CHNGAB = .TRUE. 60 CONTINUE IF (CHNGAB) THEN C DO (INITIAL-IZE) CHNGAB = .FALSE. BRACKT = .FALSE. CONVRG = .FALSE. THELT0 = .FALSE. THEGT0 = .FALSE. EIGLO = 0.0 EIGLT = 0.0 EIGRT = 0.0 EIGUP = 0.0 DTHOLD = 0.0 C END (INITIAL-IZE) END IF C C Recompute boundary conditions at singular endpoint(s). C C DO (RESET-BOUNDARY-CONDITIONS) DERIVL = 0.0 IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,ALFA,KFLAG,DERIVL) DERIVR = 0.0 IF (.NOT.BOK) THEN CALL ALFBET(B,INTAB,BB,B1,B2,EIG,P0ATB,QFATB,.FALSE., 1 BETA,JFLAG,DERIVR) BETA = PI - BETA END IF IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' alfa=',ALFA,' beta=',BETA C END (RESET-BOUNDARY-CONDITIONS) 70 CONTINUE IF (EIG.NE.GUESS .AND. .NOT.BRACKT) THEN C C If initial guess was supplied, check that boundary conditions C can be satisfied at singular endpoints. If not, try for C slightly lower EIG consistent with boundary conditions. C 80 CONTINUE IF (.NOT.AOK) CALL ALFBET(A,INTAB,AA,A1,A2,EIG, 1 P0ATA,QFATA,.FALSE.,TMP,KFLAG,TEMP) IF (.NOT.BOK) CALL ALFBET(B,INTAB,BB,B1,B2,EIG, 1 P0ATB,QFATB,.FALSE.,TMP,JFLAG,TEMP) IF (KFLAG*JFLAG.NE.1) THEN IF (THEGT0) EIG = 0.6*EIG + 0.4*EIGUP IF (THELT0) EIG = 0.6*EIG + 0.4*EIGLO IF (THELT0 .AND. EIGLO.LT.ELIM) EIGUP = ELIM GO TO 80 END IF END IF C DO (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) 90 CONTINUE IF (PRIN) WRITE(IOUT,'(/A,E22.14,A,E10.3,A,E10.3)') 1 ' guess=',EIG,' eps=',EPS,' tmid=',TMID C DO (INTEGRATE-FOR-DTHETA) C DO (SET-INITIAL-CONDITIONS) DTHDEA = DERIVL DTHDAA = 0.0 IF (.NOT.AOK) THEN CALL DXDT(AA,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDAA = -(COS(ALFA)**2/PX + 1 C*SIN(ALFA)**2)*DT C C Two special cases for DTHDAA. C IF (C.GE.0.0 .AND. P0ATA.LT.0.0 .AND. 1 QFATA.LT.0.0) DTHDAA = DTHDAA + 1 ALFA*DT/(X-A) IF (C.GE.0.0 .AND. P0ATA.GT.0.0 .AND. 1 QFATA.GT.0.0) DTHDAA = DTHDAA + 2 (ALFA-0.5*PI)*DT/(X-A) END IF DTHDEB = -DERIVR DTHDBB = 0.0 IF (.NOT.BOK) THEN CALL DXDT(BB,DT,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = EIG*RX + QX DTHDBB = -(COS(BETA)**2/PX + 1 C*SIN(BETA)**2)*DT C C Two special cases for DTHDBB. C IF (C.GE.0.0 .AND. P0ATB.LT.0.0 .AND. 1 QFATB.LT.0.0) DTHDBB = DTHDBB + 2 (PI-BETA)*DT/(B-X) IF (C.GE.0.0 .AND. P0ATB.GT.0.0 .AND. 1 QFATB.GT.0.0) DTHDBB = DTHDBB + 2 (0.5*PI-BETA)*DT/(B-X) END IF TMAX = TMID GMAX = ABS(DTHDEA) EIGSAV = EIG C END (SET-INITIAL-CONDITIONS) C T C YL = (theta,d(theta)/d(eig),d(theta)/da) C T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 C C Use integrator in one-step mode towards change to different TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 100 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.GT.-0.8 .AND. 1 ABS(YL(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YL(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 100 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDA = DTHDAA*EXP(-2.0*YL(3)) C T C YR = (theta,d(theta)/d(eig),d(theta)/db) C T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 C C Use integrator in one-step mode towards change to different TMID. C LFLAG = 1 IF (CHNGM) LFLAG = -1 110 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR, 1 WORK,IWORK) IF (LFLAG.EQ.-2 .AND. T.LT.0.8 .AND. 1 ABS(YR(2)).GT.GMAX) THEN TMAX = T GMAX = ABS(YR(2)) END IF IF (LFLAG.EQ.3 .OR. LFLAG.EQ.-2) GO TO 110 IF (LFLAG.GT.3) THEN IFLAG = 6 RETURN END IF DTHDB = DTHDBB*EXP(-2.0*YR(3)) C C DTHETA measures theta difference from left and right integrations. C DTHETA = YL(1) - YR(1) DTHDE = YL(2) - YR(2) ER1 = ERL(1) - ERR(1) ER2 = ERL(2) - ERR(2) TMID = TMAX C END (INTEGRATE-FOR-DTHETA) C C Define ONEDIG to try to be sure of one correct digit in DTHETA. C Redo integrations with tighter tolerance until ONEDIG is true. C ONEDIG = ABS(ER1).LE.0.1*ABS(DTHETA) .AND. 1 ABS(ER2).LE.0.1*ABS(DTHDE) NEWTON = ABS(DTHETA).LT.0.06 IF (NEWTON) THEN C DO (COMPUTE-CONVRG) C C Measure convergence after adding separate contributions to error. C T1 = ABS(DTHETA)+50.0*ABS(ER1) T2 = (1.0+AA)*ABS(DTHDA) T3 = (1.0-BB)*ABS(DTHDB) ESTERR = (T1+T2+T3)/ABS(DTHDE)/MAX(ONE,ABS(EIG)) CONVRG = ESTERR.LE.TAU IF (PRIN) WRITE(IOUT,'(A,L2)') 1 ' converge=',CONVRG IF (PRIN .AND. .NOT.CONVRG) 1 WRITE(IOUT,'(A,E15.7)') 2 ' estim. acc.=',ESTERR C END (COMPUTE-CONVRG) END IF IF (.NOT.ONEDIG .OR. 1 ABS(ER1).GT.0.01*ABS(DTHETA)) THEN C C Reduce local error criterion, but return IFLAG=5 if too small. C EPS = 0.05*EPS IF (EPS.LE.EPSMIN) THEN IFLAG = 5 RETURN END IF END IF IF (.NOT.(ONEDIG .OR. CONVRG)) GO TO 90 IF (ABS(DTHETA).LT.0.1) CHNGM = .FALSE. IF (PRIN) WRITE(IOUT,'(A,E15.7,A,E15.7)') 1 ' dtheta=',DTHETA,' dthde=',DTHDE IF (PRIN) WRITE(IOUT,'(/A,E15.7,A,E15.7)') 1 ' thetal=',YL(1),' thetar=',YR(1) C END (OBTAIN-DTHETA-WITH-ONE-CORRECT-DIGIT) C DO (SET-BRACKET-AND-FORM-NEWTON-ITERATES) C C EIG is bracketed when both THEGT0=.true. and THELT0=.true. C IF (DTHETA*DTHDE.GT.0.0) THEN IF (.NOT.THEGT0 .OR. EIG.LE.EIGUP) THEN THEGT0 = .TRUE. EIGUP = EIG FUP = DTHETA EIGRT = EIG - DTHETA/DTHDE END IF ELSE IF (.NOT.THELT0 .OR. EIG.GE.EIGLO) THEN THELT0 = .TRUE. EIGLO = EIG FLO = DTHETA EIGLT = EIG - DTHETA/DTHDE END IF END IF BRACKT = THELT0 .AND. THEGT0 IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' eigrt=',EIGRT,' eigup=',EIGUP IF (PRIN) WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' eiglt=',EIGLT,' eiglo=',EIGLO C END (SET-BRACKET-AND-FORM-NEWTON-ITERATES) IF (NEWTON) THEN CHNGM = .TRUE. IF (CONVRG) THEN CHNG = DTHDA*(-1.0-AA) - DTHDB*(1.0-BB) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) ELSE CHNGAB = T1.LT.0.5*(T2+T3) C C Move endpoint(s) out or take Newton step, according to CHNGAB. C IF (CHNGAB) THEN C DO (MOVE-ENDPOINTS) IF (T2.GT.T1 .AND. AA.GT.AAA) 1 AA = AA + 0.8*(-1.0-AA) IF (T3.GE.T1 .AND. BB.LT.BBB) 1 BB = BB + 0.8*(1.0-BB) AA = MAX(AA,AAA) BB = MIN(BB,BBB) IF ((AAA-AA).EQ.(BBB-BB)) THEN C C Cannot move endpoint(s) again. Store estimates and return IFLAG=5. C CHNG = (DTHDA-DTHDB)*(AAA-AA) TEMP = (DTHETA+CHNG)/DTHDE EIG = EIG - TEMP TOL = ABS(TEMP)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF EEE = EIG IF (PRIN) WRITE(IOUT,'(A,2E15.8)') 1 ' new endpoints ',AA,BB C END (MOVE-ENDPOINTS) ELSE EIG = EIG - DTHETA/DTHDE END IF END IF ELSE IF (BRACKT) THEN C C Obtain next estimate of EIG by bisection or linear interpolation. C FMAX = MAX(-FLO,FUP) EOLD = EIG RATIO = DTHETA/DTHOLD EIG = 0.5*(EIGLO+EIGUP) IF (FMAX.LE.1.5) THEN U = -FLO/(FUP-FLO) DIST = EIGUP - EIGLO EIG = EIGLO + U*DIST V = MIN(EIGLT,EIGRT) IF (EIG.LE.V) EIG = 0.5*(EIG+V) V = MAX(EIGLT,EIGRT) IF (EIG.GE.V) EIG = 0.5*(EIG+V) DE = EIG - EOLD CHNGAB = RATIO.GE.0.4 .AND. .NOT.(AOK.AND.BOK) IF (ABS(DE).LT.EPSMIN) THEN TOL = ABS(DE)/MAX(ONE,ABS(EIG)) IFLAG = 5 RETURN END IF END IF CHNGM = .NOT.CHNGM .AND. RATIO.GE.0.4 ELSE C DO (TRY-FOR-BRACKET) C C Take twice the Newton step in trying for a bracket. C IF (EIG.EQ.EEE) THEN IF (GUESS.NE.0.0) DEDW = 1.0/DTHDE CHNG = -(DEDW+1.0/DTHDE)*DTHETA IF (ABS(CHNG).GT.0.1*ABS(EIG)) 1 CHNG = -0.1*SIGN(EIG,DTHETA) ELSE CHNG = -2.0*DTHETA/DTHDE END IF LOGIC = EIG.NE.EEE .AND. DTHOLD.LT.0.0 .AND. 1 LIMIT .AND. CHNG.GT.(ELIM-EIG) IF (LOGIC) THEN CHNG = 0.99*(ELIM-EIG+EPSMIN) IF (CHNG.LT.EPSMIN) THEN C C If change is too small, EIG is presumed not to exist (IFLAG=3). C NUMEIG = NEIG - INT(-DTHETA/PI) IFLAG = 3 RETURN END IF C C Limit change in EIG to a factor of 10. C ELSE IF (ABS(CHNG).GT.(1.0+10.0*ABS(EIG))) THEN CHNG = SIGN(1.0+10.0*ABS(EIG),CHNG) ELSE IF (ABS(EIG).GE.1.0 .AND. 1 ABS(CHNG).LT.0.1*ABS(EIG)) THEN CHNG = 0.1*SIGN(EIG,CHNG) END IF EOLD = EIG EIG = EIG + CHNG C END (TRY-FOR-BRACKET) END IF END IF DTHOLD = DTHETA IF (.NOT.(CONVRG .OR. CHNGAB .OR. NEWTON)) GO TO 70 IF (.NOT.CONVRG) GO TO 60 IFLAG = 1 IF (PRIN) WRITE(IOUT,'(A,I7,A,E22.14,A,E10.3)') 1 ' numeig=',NUMEIG,' eig=',EIG,' tol=',TOL C DO (COMPUTE-EIGENFUNCTION-DATA) C C Convert from T to X values, fill 7 of first 9 locations of SLFUN. C CALL DXDT(TMID,TEMP,XMID) CALL DXDT(AA,TEMP,XAA) CALL DXDT(BB,TEMP,XBB) SLFUN(1) = XMID SLFUN(2) = XAA SLFUN(3) = ALFA SLFUN(5) = XBB SLFUN(6) = BETA + PI*NEIG SLFUN(8) = EPS SLFUN(9) = Z C C Compute SLFUN(4), SLFUN(7) towards normalizing the eigenfunction. C EIGSAV = EIG Z = -Z T = AA YL(1) = ALFA YL(2) = DTHDEA YL(3) = 0.0 LFLAG = 1 120 CONTINUE CALL GERK(F,3,YL,T,TMID,EPS,EPS,LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 120 T = BB YR(1) = BETA + PI*NEIG YR(2) = DTHDEB YR(3) = 0.0 LFLAG = 1 130 CONTINUE CALL GERK(F,3,YR,T,TMID,EPS,EPS,LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 130 Z = -Z SL = SIN(YL(1)) SR = SIN(YR(1)) CL = COS(YL(1)) CR = COS(YR(1)) UL = (YL(2)-DTHDEA*EXP(-2.0*YL(3)))*Z UR = (YR(2)-DTHDEB*EXP(-2.0*YR(3)))*Z ASL = ABS(SL) ASR = ABS(SR) DEN = 0.5*LOG(UL*ASR*ASR-UR*ASL*ASL) SLFUN(4) = LOG(ASR) - YL(3) - DEN SLFUN(7) = LOG(ASL) - YR(3) - DEN C END (COMPUTE-EIGENFUNCTION-DATA) C DO (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) C C Perform final check on EIG. Return IFLAG=4 if not accurate enough. C E = ASR*EXP(-DEN) PSIL = E*SL PSIPL = E*CL SQL = E*E*UL DPSIL = PSIL*ERL(3) + PSIPL*ERL(1) DPSIPL = PSIL*ERL(1) + PSIPL*ERL(3) PSIPL = PSIPL*Z E = ASL*EXP(-DEN) PSIR = E*SR PSIPR = E*CR SQR = E*E*UR DPSIR = PSIR*ERR(3) + PSIPR*ERR(1) DPSIPR = PSIR*ERR(1) + PSIPR*ERR(3) PSIPR = PSIPR*Z RAY = EIG + (PSIL*PSIPL-PSIR*PSIPR)/(SQL-SQR) IF (PRIN) THEN WRITE(IOUT,'(A,E22.14)') ' ray=',RAY WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' psil=',PSIL,' psir=',PSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' psipl=',PSIPL,' psipr=',PSIPR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' sql=',SQL,' sqr=',SQR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' dpsil=',DPSIL,' dpsir=',DPSIR WRITE(IOUT,'(A,E22.14,A,E22.14)') 1 ' dpsipl=',DPSIPL,' dpsipr=',DPSIPR END IF C END (CHECK-MATCHING-VALUES-OF-EIGENFUNCTION) IF (ABS(RAY-EIG).GT.TAU*MAX(ONE,ABS(EIG))) IFLAG = 4 IF (ISLFUN.GT.0) THEN C C Calculate selected eigenfunction values by integration. C C DO (GENERATE-EIGENFUNCTION-VALUES) EIGF = .TRUE. NMID = 0 DO 140 I=1,ISLFUN IF (SLFUN(9+I).LE.SLFUN(1)) NMID = I 140 CONTINUE IF (NMID.GT.0) THEN X = SLFUN(2) YL(1) = SLFUN(3) YL(2) = 0.0 YL(3) = SLFUN(4) LFLAG = 1 DO 160 J=1,NMID 150 CONTINUE CALL GERK(F,3,YL,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERL,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 150 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YL(3))*SIN(YL(1)) 160 CONTINUE END IF IF (NMID.LT.ISLFUN) THEN X = SLFUN(5) YR(1) = SLFUN(6) YR(2) = 0.0 YR(3) = SLFUN(7) LFLAG = 1 DO 180 J=ISLFUN,NMID+1,-1 170 CONTINUE CALL GERK(F,3,YR,X,SLFUN(J+9),SLFUN(8),SLFUN(8), 1 LFLAG,ERR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 170 IF (LFLAG.EQ.6) LFLAG = 2 SLFUN(J+9) = EXP(YR(3))*SIN(YR(1)) 180 CONTINUE END IF C END (GENERATE-EIGENFUNCTION-VALUES) END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE ALFBET(XEND,INTAB,TT,COEF1,COEF2,EIG,P0,QF,OK, 1 VALUE,IFLAG,DERIV) INTEGER INTAB,IFLAG LOGICAL OK DOUBLE PRECISION XEND,TT,COEF1,COEF2,EIG,P0,QF,VALUE,DERIV C ********** C C This subroutine computes a boundary value for a specified endpoint C of the interval for a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is called C from SLEIGN. Both regular and singular endpoints are treated. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt,extrap C C ********** C .. Scalars in Common .. DOUBLE PRECISION Z C .. C .. Local Scalars .. LOGICAL LOGIC DOUBLE PRECISION C,CD,D,HH,ONE,PI,PX,QX,RX,T,TEMP,X C .. C .. External Functions .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT,EXTRAP C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SIGN,SQRT C .. C .. Common blocks .. COMMON /ZEE/Z C .. C Set machine dependent constant. C C PI (variable ONE set to 1.0 eases precision conversion). ONE = 1.0 PI = 4.0*ATAN(ONE) C IFLAG = 1 DERIV = 0.0 IF (OK) THEN VALUE = 0.5*PI IF (COEF1.NE.0.0) VALUE = ATAN(-Z*COEF2/COEF1) LOGIC = (TT.LT.0.0 .AND. VALUE.LT.0.0) .OR. 1 (TT.GT.0.0 .AND. VALUE.LE.0.0) IF (LOGIC) VALUE = VALUE + PI ELSE LOGIC = (INTAB.EQ.2 .AND. TT.GT.0.0) .OR. 1 (INTAB.EQ.3 .AND. TT.LT.0.0) .OR. 2 INTAB.EQ.4 .OR. (P0.GT.0.0 .AND. QF.LT.0.0) IF (LOGIC) THEN T = SIGN(ONE,TT) CALL EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) ELSE CALL DXDT(TT,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z C = 2.0*(EIG*RX+QX) IF (C.LT.0.0) THEN VALUE = 0.0 IF (P0.GT.0.0) VALUE = 0.5*PI ELSE HH = ABS(XEND-X) D = 2.0*HH/PX CD = C*D*HH IF (P0.GT.0.0) THEN VALUE = C*HH IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) VALUE = VALUE + 0.5*PI ELSE VALUE = D IF (CD.LT.1.0) VALUE = VALUE/(1.0+SQRT(1.0-CD)) END IF END IF END IF END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE DXDT(T,DT,X) DOUBLE PRECISION T,DT,X C ********** C C This subroutine transforms coordinates from T on (-1,1) to C X on (A,B) in the solution of a Sturm-Liouville problem. C It is called from subroutines SLEIGN, ALFBET, F, and EXTRAP. C C ********** C .. Scalars in Common .. INTEGER INTAB DOUBLE PRECISION A,B,C1,C2 C .. C .. Local Scalars .. DOUBLE PRECISION U C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /DATADT/A,B,C1,C2,INTAB C .. U = C1*T + C2 GO TO (10,20,30,40), INTAB 10 CONTINUE DT = C1*0.5*(B-A) X = 0.5*((B+A)+(B-A)*U) RETURN 20 CONTINUE DT = C1*2.0/(1.0-U)**2 X = A + (1.0+U)/(1.0-U) RETURN 30 CONTINUE DT = C1*2.0/(1.0+U)**2 X = B - (1.0-U)/(1.0+U) RETURN 40 CONTINUE DT = C1/(1.0-ABS(U))**2 X = U/(1.0-ABS(U)) RETURN END C ---------------------------------------------------------------------- SUBROUTINE ESTEIG(MF,ML,LIMIT,ELIM,EMAX,EMIN,PIN,QS,RS,DS,PS, 1 IMAX,IMIN,EEE,EIG,IA,IB,EL,WL,DEDW) INTEGER MF,ML,IMAX,IMIN,IA,IB LOGICAL LIMIT DOUBLE PRECISION ELIM,EMAX,EMIN,PIN,EEE,EIG,EL,WL,DEDW DOUBLE PRECISION QS(ML),RS(ML),DS(ML),PS(ML) C ********** C C This subroutine generates an initial guess for a specified C eigenvalue of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is C called from SLEIGN when no initial guess is provided by the user. C C The method used is to approximately solve the equation for EIG C C Integral (sqrt((eig*r+q)/p)) = numeig*pi C C where the integral is taken over those X in (A,B) for which C C (eig*r+q)/p .gt. 0 C C and NUMEIG is the index of the desired eigenvalue (PIN=NUMEIG*pi). C C Subprograms called C C sleign-supplied ... estpac C C ********** C .. Scalars in Common .. INTEGER INTAB DOUBLE PRECISION A,B,C1,C2 C .. C .. Local Scalars .. INTEGER IE,IP LOGICAL LOGIC DOUBLE PRECISION BALLPK,EU,FNEW,FOLD,SUM,TEMP,U,WU C .. C .. External Subroutines .. EXTERNAL ESTPAC C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. C .. Common blocks .. COMMON /DATADT/A,B,C1,C2,INTAB C .. EEE = MIN(ELIM,EMAX) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) C C Choose bounds for EIG and associate function (integral) values. C EL = EMIN WL = 0.0 EU = EEE WU = SUM IF (LIMIT .AND. WU.LT.PIN) THEN EIG = ELIM ELSE IF (U.EQ.0.0) THEN EL = EMAX EEE = EMAX + 1.0 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM END IF 10 CONTINUE IF (WU.LE.PIN) THEN C C Increase trial value if integral is still too small. C EL = EU WL = WU EEE = EU + ((PIN-WU+3.0)/U)**2 C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) EU = EEE WU = SUM GO TO 10 END IF C C EIG is bracketed. Now try to reduce the size of the bracket C by searching among the saved values of -QS()/RS(). C 20 CONTINUE IF (ABS(IMAX-IMIN).GE.2 .AND. EU.LE.EMAX) THEN IE = (IMAX+IMIN)/2 IF (RS(IE).NE.0.0) THEN EEE = -QS(IE)/RS(IE) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.GT.PIN) THEN IMAX = IE WU = SUM EU = EEE ELSE IMIN = IE WL = SUM EL = EEE END IF GO TO 20 END IF END IF C C Improve approximation for EIG using bisection or secant method. C Substitute 'ballpark' estimate if approximation grows too large. C DEDW = (EU-EL)/(WU-WL) FOLD = 0.0 IF (INTAB.EQ.1) BALLPK = (PIN/(A-B))**2 LOGIC = .TRUE. 30 CONTINUE IF (LOGIC) THEN LOGIC = (WL.LT.(PIN-1.0) .OR. WU.GT.(PIN+1.0)) EEE = EL + DEDW*(PIN-WL) FNEW = MIN(PIN-WL,WU-PIN) IF (FNEW.GT.0.4*FOLD .OR. FNEW.LE.1.0) 1 EEE = 0.5*(EL+EU) IF (INTAB.EQ.1 .AND. ABS(EEE).GT.1.0E+3*BALLPK) THEN EIG = BALLPK RETURN ELSE IF (INTAB.NE.1 .AND. ABS(EEE).GT.1.0E+6) THEN EIG = 1.0 RETURN ELSE FOLD = FNEW C DO (ESTIMATE-PHASE-ANGLE-CHANGE) CALL ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS, 1 IA,IB,IP,SUM,U,TEMP) C END (ESTIMATE-PHASE-ANGLE-CHANGE) IF (SUM.LT.PIN) THEN EL = EEE WL = SUM ELSE EU = EEE WU = SUM END IF DEDW = (EU-EL)/(WU-WL) GO TO 30 END IF END IF END IF RETURN END C ---------------------------------------------------------------------- SUBROUTINE ESTPAC(MF,ML,EEE,PIN,QS,RS,DS,PS,IA,IB,IP,SUM,U,ZAV) INTEGER MF,ML,IA,IB,IP DOUBLE PRECISION EEE,PIN,SUM,U,ZAV DOUBLE PRECISION QS(ML),RS(ML),DS(ML),PS(ML) C ********** C C This subroutine estimates the change in 'phase angle' in the C eigenvalue determination of a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. It is C called from SLEIGN, and also from ESTEIG when no initial guess C is provided by the user. C C The subroutine approximates the (trapezoidal rule) integral of C C sqrt((eig*r+q)/p) C C where the integral is taken over those X in (A,B) for which C C (eig*r+q)/p .gt. 0 C C ********** C .. Local Scalars .. INTEGER J DOUBLE PRECISION PSUM,RT,RTSAV,V,W,WSAV C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN,SQRT C .. IA = MF IB = 80 IP = MF C C SUM accumulates the integral approximation. U measures the total C length of subintervals where (EIG*R+Q)/P .gt. 0.0. ZAV is the C average value of sqrt((EIG*R+Q)*P) over those subintervals. C SUM = 0.0 U = 0.0 ZAV = 0.0 WSAV = EEE*RS(MF) + QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 10 J=MF+1,ML W = EEE*RS(J) + QS(J) IF (W.GT.0.0) THEN IF (J.GE.80) IB = J U = U + DS(J-1) RT = SQRT(W) ELSE RT = 0.0 IF (U.EQ.0.0 .AND. RTSAV.EQ.0.0 .AND. IA.LE.19) IA = IA + 1 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(J-1)*V IF (PSUM.LT.(PIN-SUM)) IP = J SUM = SUM + PSUM IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(J)+PS(J-1)) 10 CONTINUE SUM = 0.5*SUM ZAV = 0.25*ZAV RETURN END C ---------------------------------------------------------------------- SUBROUTINE EXTRAP(T,TT,EIG,VALUE,DERIV,IFLAG) INTEGER IFLAG DOUBLE PRECISION T,TT,EIG,VALUE,DERIV C ********** C C This subroutine is called from ALFBET in determining boundary C values at a singular endpoint of the interval for a C Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C C EXTRAP, which in turn calls INTPOL, extrapolates the function C C arctan(1.0/sqrt(-p*(eig*r+q))) C C from its values for T within (-1,1) to an endpoint. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt,intpol C C ********** C .. Scalars in Common .. DOUBLE PRECISION Z C .. C .. Local Scalars .. INTEGER KGOOD DOUBLE PRECISION ANS,CTN,ERROR,PROD,PX,QX,RX,T1,TEMP,X C .. C .. Local Arrays .. DOUBLE PRECISION FN1(5),XN(5) C .. C .. External Functions .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT,INTPOL C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,SQRT,TAN C .. C .. Common blocks .. COMMON /ZEE/Z C .. IFLAG = 1 KGOOD = 0 T1 = TT 10 CONTINUE CALL DXDT(T1,TEMP,X) PX = P(X)/Z QX = Q(X)/Z RX = R(X)/Z PROD = -PX*(EIG*RX+QX) IF (PROD.LE.0.0) THEN T1 = 0.5*(T1+T) IF ((1.0+(T1-T)**2).GT.1.0) GO TO 10 IFLAG = 5 RETURN ELSE KGOOD = KGOOD + 1 XN(KGOOD) = T1 FN1(KGOOD) = ATAN(1.0/SQRT(PROD)) T1 = 0.5*(T+T1) IF (KGOOD.LT.5) GO TO 10 END IF T1 = 0.01 CALL INTPOL(5,XN,FN1,T,T1,3,ANS,ERROR) VALUE = ABS(ANS) CTN = 1.0/TAN(VALUE) DERIV = 0.5*PX*RX/CTN/(1.0+CTN**2) TT = XN(1) RETURN END C ---------------------------------------------------------------------- SUBROUTINE F(T,Y,YP) DOUBLE PRECISION T DOUBLE PRECISION Y(2),YP(3) C ********** C C This subroutine evaluates the derivative functions for use with C integrator GERK in solving a Sturm-Liouville problem in the form C C (p(x)*y'(x))' + (q(x) + eig*r(x))*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P, Q, and R. C C Subprograms called C C user-supplied ..... p,q,r C C sleign-supplied ... dxdt C C ********** C .. Scalars in Common .. LOGICAL EIGF DOUBLE PRECISION EIG,Z C .. C .. Local Scalars .. DOUBLE PRECISION C,C2,DT,QX,RX,S,S2,V,W,X,XP,ZP C .. C .. External Functions .. DOUBLE PRECISION P,Q,R EXTERNAL P,Q,R C .. C .. External Subroutines .. EXTERNAL DXDT C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,SIN C .. C .. Common blocks .. COMMON /DATAF/EIG,EIGF COMMON /ZEE/Z C .. DT = 1.0 X = T IF (.NOT.EIGF) CALL DXDT(T,DT,X) ZP = ABS(Z) XP = ZP/P(X) QX = Q(X)/ZP RX = R(X)/ZP V = EIG*RX + QX S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = DT*(XP*C2+V*S2) W = (XP-V)*S*C IF (Z.LT.0.0) RX = ABS(RX) YP(2) = DT*(-2.0*W*Y(2)+RX*S2) YP(3) = DT*W RETURN END C ---------------------------------------------------------------------- SUBROUTINE G(X,Y,YP) DOUBLE PRECISION X DOUBLE PRECISION Y(1),YP(1) C ********** C C This subroutine evaluates the derivative function for use with C integrator GERK in solving a differential equation in the form C C (p(x)*y'(x))' + q(x)*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P and Q. C C Subprograms called C C user-supplied ..... p,q C C ********** C .. Scalars in Common .. DOUBLE PRECISION Z C .. C .. Local Scalars .. DOUBLE PRECISION C,C2,QX,S,S2,XP C .. C .. External Functions .. DOUBLE PRECISION P,Q EXTERNAL P,Q C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN C .. C .. Common blocks .. COMMON /ZEE/Z C .. XP = Z/P(X) QX = Q(X)/Z S = SIN(Y(1)) C = COS(Y(1)) S2 = S*S C2 = C*C YP(1) = XP*C2+QX*S2 RETURN END C ---------------------------------------------------------------------- SUBROUTINE INTPOL(N,XN,FN,X,ABSERR,MAXDEG,ANS,ERROR) INTEGER N,MAXDEG DOUBLE PRECISION X,ABSERR,ANS,ERROR DOUBLE PRECISION XN(N),FN(N) C ********** C C This subroutine forms an interpolating polynomial for data pairs. C It is called from EXTRAP in solving a Sturm-Liouville problem. C C ********** C .. Local Scalars .. INTEGER I,I1,II,IJ,IK,IKM1,J,K,L,LIMIT DOUBLE PRECISION PROD C .. C .. Local Arrays .. INTEGER INDEX(10) DOUBLE PRECISION V(10,10) C .. C .. Intrinsic Functions .. INTRINSIC ABS,MIN C .. L = MIN(MAXDEG,N-2) + 2 LIMIT = MIN(L,N-1) DO 10 I = 1,N V(I,1) = ABS(XN(I)-X) INDEX(I) = I 10 CONTINUE DO 30 I=1,LIMIT DO 20 J=I+1,N II = INDEX(I) IJ = INDEX(J) IF (V(II,1).GT.V(IJ,1)) THEN INDEX(I) = IJ INDEX(J) = II END IF 20 CONTINUE 30 CONTINUE PROD = 1.0 I1 = INDEX(1) ANS = FN(I1) V(1,1) = FN(I1) DO 50 K=2,L IK = INDEX(K) V(K,1) = FN(IK) DO 40 I=1,K-1 II = INDEX(I) V(K,I+1) = (V(I,I)-V(K,I))/(XN(II)-XN(IK)) 40 CONTINUE IKM1 = INDEX(K-1) PROD = (X-XN(IKM1))*PROD ERROR = PROD*V(K,K) IF(ABS(ERROR).LE.ABSERR) RETURN ANS = ANS + ERROR 50 CONTINUE ANS = ANS - ERROR RETURN END C ---------------------------------------------------------------------- SUBROUTINE ZCOUNT(A,B,A1,A2,B1,B2,JPAIRS,PAIRS,JSUM) INTEGER JPAIRS,JSUM DOUBLE PRECISION A,B,A1,A2,B1,B2 DOUBLE PRECISION PAIRS(2*JPAIRS) C ********** C C This subroutine counts the zeros, over specified subintervals, of C the solutions of a second order differential equation in the form C C (p(x)*y'(x))' + q(x)*y(x) = 0 on (a,b) C C for user-supplied coefficient functions P and Q. This count in C turn corresponds to the number of zeros, in the interior of (A,B), C of the first eigenfunction of the related Sturm-Liouville problem C whose (semidefinite) weight function vanishes identically in the C subintervals. The problem is restricted to be regular. C C The applicable initial condition depends upon three cases. C C Case 1 -- On a subinterval with left endpoint A, C A1*Y(A) + A2*Y'(A)*P(A) = 0. C C Case 2 -- On a subinterval with right endpoint B, C B1*Y(B) + B2*Y'(B)*P(B) = 0. C C Case 3 -- On a subinterval with neither A nor B as endpoint, C Y(XAA) = 0, where XAA is the left endpoint. C C The SUBROUTINE statement is C C SUBROUTINE zcount(a,b,a1,a2,b1,b2,jpairs,pairs,jsum) C C where C C A and B are input variables defining the full interval. C A must be less than B. C C A1 and A2 are input variables set to prescribe the initial C condition at A (Case 1). C C B1 and B2 are input variables set to prescribe the initial C condition at B (Case 2). C C JPAIRS is an integer input variable set to the number of C specified subintervals of (A,B). C C PAIRS is an input array of length 2*JPAIRS whose successive C ordered element pairs specify the subintervals. C C JSUM is an integer output variable set to the total zero count. C C Subprograms called C C sleign-supplied ... epslon,g,gerk C C user-supplied ..... p,q C C This version dated 5/18/89. C Burton S. Garbow, Argonne National Laboratory C C ********** C .. Scalars in Common .. DOUBLE PRECISION Z C .. C .. Local Scalars .. INTEGER I,J,LFLAG,MF,ML DOUBLE PRECISION EPS,EPSMIN,H,ONE,PI,PSUM,PX,QX,RT,RTSAV, 1 T,U,V,W,WSAV,X,X50,XAA,XBB,XSAV,ZAV C .. C .. Local Arrays .. INTEGER IWORK(5) DOUBLE PRECISION DS(98),GERROR(1),PS(99),QS(99),WORK(11),Y(1) C .. C .. External Functions .. DOUBLE PRECISION EPSLON,P,Q EXTERNAL EPSLON,P,Q C .. C .. External Subroutines .. EXTERNAL G,GERK C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,INT,SIGN,SQRT C .. C .. Common blocks .. COMMON /ZEE/Z C .. C Set constants EPSMIN, the computer unit roundoff error, and PI. C (Variable ONE set to 1.0 eases precision conversion.) C ONE = 1.0 EPSMIN = EPSLON(ONE) PI = 4.0*ATAN(ONE) C C Set relative and absolute error tolerances for GERK. C EPS = SQRT(EPSMIN) C JSUM = 0 DO 70 J = 1,JPAIRS XAA = PAIRS(2*J-1) XBB = PAIRS(2*J) C DO (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) C C Evaluate P, Q to obtain preliminary information about the C differential equation. C C DO (SAMPLE-COEFFICIENTS) X50 = 0.5*((XBB+XAA)+(XBB-XAA)*EPSMIN) PX = P(X50) QX = Q(X50) PS(50) = PX QS(50) = QX/PX C C MF and ML are the least and greatest index values, respectively. C XSAV = X50 H = 0.9/40.0 DO 10 I=49,1,-1 IF (I.GE.10) T = H*(I-50) IF (I.LT.10) T = T - 0.75*(1.0+T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I) = XSAV - X XSAV = X MF = I 10 CONTINUE XSAV = X50 DO 20 I=51,99 IF (I.LE.90) T = H*(I-50) IF (I.GT.90) T = T + 0.75*(1.0-T) X = 0.5*((XBB+XAA)+(XBB-XAA)*T) PX = P(X) QX = Q(X) PS(I) = PX QS(I) = QX/PX DS(I-1) = X - XSAV XSAV = X ML = I - 1 20 CONTINUE C END (SAMPLE-COEFFICIENTS) C DO (ESTIMATE-PHASE-ANGLE-CHANGE) C C U measures the total length of subintervals where Q/P .gt. 0.0. C ZAV is the average value of sqrt(Q*P) over those subintervals. C U = 0.0 ZAV = 0.0 WSAV = QS(MF) IF (WSAV.GT.0.0) THEN RTSAV = SQRT(WSAV) ELSE RTSAV = 0.0 END IF DO 30 I=MF+1,ML W = QS(I) IF (W.GT.0.0) THEN U = U + DS(I-1) RT = SQRT(W) ELSE RT = 0.0 END IF IF (W.EQ.0.0 .OR. WSAV.EQ.0.0 .OR. 1 W.EQ.SIGN(W,WSAV)) THEN V = RT + RTSAV ELSE V = (W*RT+WSAV*RTSAV)/ABS(W-WSAV) END IF WSAV = W RTSAV = RT PSUM = DS(I-1)*V IF (U.GT.0.0) ZAV = ZAV + PSUM*(PS(I)+PS(I-1)) 30 CONTINUE ZAV = 0.25*ZAV C END (ESTIMATE-PHASE-ANGLE-CHANGE) Z = 1.0 IF (U.GT.0.0) Z = ZAV/U C END (CALCULATE-MODIFIED-PRUFER-TRANSFORMATION-CONSTANT) LFLAG = 1 IF (XAA.EQ.A) THEN C C Case 1 ---------- C Y(1) = PI/2.0 IF (A1.NE.0.0) Y(1) = ATAN(-Z*A2/A1) IF (Y(1).LT.0.0) Y(1) = Y(1) + PI 40 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 40 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) ELSE IF (XBB.EQ.B) THEN C C Case 2 ---------- C Y(1) = PI/2.0 IF (B1.NE.0.0) Y(1) = ATAN(-Z*B2/B1) IF (Y(1).GT.0.0) Y(1) = Y(1) - PI 50 CONTINUE CALL GERK(G,1,Y,XBB,XAA,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 50 JSUM = JSUM - INT((Y(1)-ABS(EPS))/PI) ELSE C C Case 3 ---------- C Y(1) = 0.0 60 CONTINUE CALL GERK(G,1,Y,XAA,XBB,EPS,EPS,LFLAG,GERROR,WORK,IWORK) IF (LFLAG.EQ.3) GO TO 60 JSUM = JSUM + INT((Y(1)+ABS(EPS))/PI) END IF 70 CONTINUE RETURN END SUBROUTINE GERK(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK, IWORK) C C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C C WRITTEN BY H.A.WATTS AND L.F.SHAMPINE C SANDIA LABORATORIES C C GERK IS DESIGNED TO SOLVE SYSTEMS OF DIFFERENTIAL EQUATIONS C WHEN IT IS IMPORTANT TO HAVE A READILY AVAILABLE GLOBAL ERROR C ESTIMATE. PARALLEL INTEGRATION IS PERFORMED TO YIELD TWO C SOLUTIONS ON DIFFERENT MESH SPACINGS AND GLOBAL EXTRAPOLATION C IS APPLIED TO PROVIDE AN ESTIMATE OF THE GLOBAL ERROR IN THE C MORE ACCURATE SOLUTION. C C FOR IBM SYSTEM 360 AND 370 AND OTHER MACHINES OF SIMILAR C ARITHMETIC CHARACTERISTICS, THIS CODE SHOULD BE CONVERTED TO C DOUBLE PRECISION. C C******************************************************************* C ABSTRACT C******************************************************************* C C SUBROUTINE GERK INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN)) C WHERE THE Y(I) ARE GIVEN AT T . C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TOUT C BUT IT CAN BE USED AS A ONE-STEP INTEGRATOR TO ADVANCE THE C SOLUTION A SINGLE STEP IN THE DIRECTION OF TOUT. ON RETURN, AN C ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T IS PROVIDED C AND THE PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE C INTEGRATION. THE USER HAS ONLY TO CALL GERK AGAIN (AND PERHAPS C DEFINE A NEW VALUE FOR TOUT). ACTUALLY, GERK IS MERELY AN C INTERFACING ROUTINE WHICH ALLOCATES VIRTUAL STORAGE IN THE C ARRAYS WORK, IWORK AND CALLS SUBROUTINE GERKS FOR THE SOLUTION. C GERKS IN TURN CALLS SUBROUTINE FEHL WHICH COMPUTES AN APPROX- C IMATE SOLUTION OVER ONE STEP. C C GERK USES THE RUNGE-KUTTA-FEHLBERG (4,5) METHOD DESCRIBED C IN THE REFERENCE C E.FEHLBERG , LOW-ORDER CLASSICAL RUNGE-KUTTA FORMULAS WITH C STEPSIZE CONTROL , NASA TR R-315 C C C THE PARAMETERS REPRESENT- C F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES C YP(I)=DY(I)/DT C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- SOLUTION VECTOR AT T C T -- INDEPENDENT VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR C LOCAL ERROR TEST. AT EACH STEP THE CODE REQUIRES THAT C ABS(LOCAL ERROR) .LE. RELERR*ABS(Y) + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION C VECTORS. C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION C GERROR(*) -- VECTOR WHICH ESTIMATES THE GLOBAL ERROR AT T. C THAT IS, GERROR(I) APPROXIMATES Y(I)-TRUE C SOLUTION(I). C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO GERK WHICH C IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED C AT LEAST 3+8*NEQN C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL C TO GERK WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST C BE DIMENSIONED AT LEAST 5 C C C******************************************************************* C FIRST CALL TO GERK C******************************************************************* C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE C ARRAYS IN THE CALL LIST - Y(NEQN), WORK(3+8*NEQN), IWORK(5), C DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE F(T,Y,YP) C AND INITIALIZE THE FOLLOWING PARAMETERS- C C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED. (NEQN .GE. 1) C Y(*) -- VECTOR OF INITIAL CONDITIONS C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE C TOUT -- OUTPUT POINT AT WHICH SOLUTION IS DESIRED. C T=TOUT IS ALLOWED ON THE FIRST CALL ONLY,IN WHICH CASE C GERK RETURNS WITH IFLAG=2 IF CONTINUATION IS POSSIBLE. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C WHICH MUST BE NON-NEGATIVE BUT MAY BE CONSTANTS. WE CAN C USUALLY EXPECT THE GLOBAL ERRORS TO BE SOMEWHAT SMALLER C THAN THE REQUESTED LOCAL ERROR TOLERANCES. TO AVOID C LIMITING PRECISION DIFFICULTIES THE CODE ALWAYS USES C THE LARGER OF RELERR AND AN INTERNAL RELATIVE ERROR C PARAMETER WHICH IS MACHINE DEPENDENT. C IFLAG -- +1,-1 INDICATOR TO INITIALIZE THE CODE FOR EACH NEW C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG= C -1 ONLY WHEN ONE-STEP INTEGRATOR CONTROL IS ESSENTIAL. C IN THIS CASE, GERK ATTEMPTS TO ADVANCE THE SOLUTION A C SINGLE STEP IN THE DIRECTION OF TOUT EACH TIME IT IS C CALLED. SINCE THIS MODE OF OPERATION RESULTS IN EXTRA C COMPUTING OVERHEAD, IT SHOULD BE AVOIDED UNLESS NEEDED. C C C******************************************************************* C OUTPUT FROM GERK C******************************************************************* C C Y(*) -- SOLUTION AT T C T -- LAST POINT REACHED IN INTEGRATION. C IFLAG = 2 -- INTEGRATION REACHED TOUT. INDICATES SUCCESSFUL C RETURN AND IS THE NORMAL MODE FOR CONTINUING C INTEGRATION. C =-2 -- A SINGLE SUCCESSFUL STEP IN THE DIRECTION OF C TOUT HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING C INTEGRATION ONE STEP AT A TIME. C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN C 9000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS C IS APPROXIMATELY 500 STEPS. C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION C VANISHED MAKING A PURE RELATIVE ERROR TEST C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE. C USING THE ONE-STEP INTEGRATION MODE FOR ONE STEP C IS A GOOD WAY TO PROCEED. C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE C ATTEMPTED. C = 6 -- GERK IS BEING USED INEFFICIENTLY IN SOLVING C THIS PROBLEM. TOO MUCH OUTPUT IS RESTRICTING THE C NATURAL STEPSIZE CHOICE. USE THE ONE-STEP C INTEGRATOR MODE. C = 7 -- INVALID INPUT PARAMETERS C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS C SATISFIED - NEQN .LE. 0 C T=TOUT AND IFLAG .NE. +1 OR -1 C RELERR OR ABSERR .LT. 0. C IFLAG .EQ. 0 OR .LT. -2 OR .GT. 7 C GERROR(*) -- ESTIMATE OF THE GLOBAL ERROR IN THE SOLUTION AT T C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO C INTEREST TO THE USER BUT NECESSARY FOR SUBSEQUENT C CALLS. WORK(1),...,WORK(NEQN) CONTAIN THE FIRST C DERIVATIVES OF THE SOLUTION VECTOR Y AT T. C WORK(NEQN+1) CONTAINS THE STEPSIZE H TO BE C ATTEMPTED ON THE NEXT STEP. IWORK(1) CONTAINS C THE DERIVATIVE EVALUATION COUNTER. C C C******************************************************************* C SUBSEQUENT CALLS TO GERK C******************************************************************* C C SUBROUTINE GERK RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT, THE USER NEED C ONLY DEFINE A NEW TOUT AND CALL GERK AGAIN. IN THE ONE-STEP C INTEGRATOR MODE (IFLAG=-2) THE USER MUST KEEP IN MIND THAT EACH C STEP TAKEN IS IN THE DIRECTION OF THE CURRENT TOUT. UPON C REACHING TOUT (INDICATED BY CHANGING IFLAG TO 2), THE USER MUST C THEN DEFINE A NEW TOUT AND RESET IFLAG TO -2 TO CONTINUE IN THE C ONE-STEP INTEGRATOR MODE. C C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS C TO CONTINUE (IFLAG=3 CASE), HE JUST CALLS GERK AGAIN. THE C FUNCTION COUNTER IS THEN RESET TO 0 AND ANOTHER 9000 FUNCTION C EVALUATIONS ARE ALLOWED. C C HOWEVER, IN THE CASE IFLAG=4, THE USER MUST FIRST ALTER THE C ERROR CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE C INTEGRATION CAN PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED. C C ALSO, IN THE CASE IFLAG=5, IT IS NECESSARY FOR THE USER TO C RESET IFLAG TO 2 (OR -2 WHEN THE ONE-STEP INTEGRATION MODE IS C BEING USED) AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH C BEFORE THE INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, C EXECUTION WILL BE TERMINATED. THE OCCURRENCE OF IFLAG=5 C INDICATES A TROUBLE SPOT (SOLUTION IS CHANGING RAPIDLY, C SINGULARITY MAY BE PRESENT) AND IT OFTEN IS INADVISABLE TO C CONTINUE. C C IF IFLAG=6 IS ENCOUNTERED, THE USER SHOULD USE THE ONE-STEP C INTEGRATION MODE WITH THE STEPSIZE DETERMINED BY THE CODE. IF C THE USER INSISTS UPON CONTINUING THE INTEGRATION WITH GERK IN C THE INTERVAL MODE, HE MUST RESET IFLAG TO 2 BEFORE CALLING GERK C AGAIN. OTHERWISE,EXECUTION WILL BE TERMINATED. C C IF IFLAG=7 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS C THE INVALID INPUT PARAMETERS ARE CORRECTED. C C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN C INFORMATION REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, C WORK AND IWORK SHOULD NOT BE ALTERED. C C******************************************************************* C C .. SCALAR ARGUMENTS .. INTEGER IFLAG,NEQN DOUBLE PRECISION ABSERR,RELERR,T,TOUT C .. C .. ARRAY ARGUMENTS .. INTEGER IWORK(5) DOUBLE PRECISION GERROR(NEQN),WORK(3+8*NEQN),Y(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K1,K1M,K2,K3,K4,K5,K6,K7,K8 C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL GERKS C .. C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY K1M = NEQN + 1 K1 = K1M + 1 K2 = K1 + NEQN K3 = K2 + NEQN K4 = K3 + NEQN K5 = K4 + NEQN K6 = K5 + NEQN K7 = K6 + NEQN K8 = K7 + NEQN C ******************************************************************* C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER, C HE MUST USE GERKS DIRECTLY. C ******************************************************************* CALL GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, WORK(1), WORK(K1M), WORK(K1), WORK(K2), WORK(K3), * WORK(K4), WORK(K5), WORK(K6), WORK(K7), WORK(K8), * WORK(K8+1), IWORK(1), IWORK(2), IWORK(3), IWORK(4), IWORK(5)) RETURN END SUBROUTINE GERKS(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, * GERROR, YP, H, F1, F2, F3, F4, F5, YG, YGP, SAVRE, SAVAE, * NFE, KOP, INIT, JFLAG, KFLAG) C FEHLBERG FOURTH(FIFTH) ORDER RUNGE-KUTTA METHOD WITH C GLOBAL ERROR ASSESSMENT C ******************************************************************* C GERKS INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS AS DESCRIBED IN THE COMMENTS FOR GERK. THE ARRAYS C YP,F1,F2,F3,F4,F5,YG AND YGP (OF DIMENSION AT LEAST NEQN) AND C THE VARIABLES H,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,AND KFLAG ARE C USED INTERNALLY BY THE CODE AND APPEAR IN THE CALL LIST TO C ELIMINATE LOCAL RETENTION OF VARIABLES BETWEEN CALLS. C ACCORDINGLY, THEY SHOULD NOT BE ALTERED. ITEMS OF POSSIBLE C INTEREST ARE C YP - DERIVATIVE OF SOLUTION VECTOR AT T C H - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP C NFE- COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION C EVALUATIONS. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER IFLAG,INIT,JFLAG,KFLAG,KOP,NEQN,NFE DOUBLE PRECISION ABSERR,H,RELERR,SAVAE,SAVRE,T,TOUT C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 GERROR(NEQN),Y(NEQN),YG(NEQN),YGP(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K,MAXNFE,MFLAG LOGICAL HFAILD,OUTPUT DOUBLE PRECISION A,AE,DT,EE,EEOET,ESTTOL,ET,HH,HMIN,ONE,REMIN,RER, 1 S,SCALE,TOL,TOLN,TS,U,U26,YPK C .. C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION EPSLON EXTERNAL EPSLON C .. C .. EXTERNAL SUBROUTINES .. EXTERNAL FEHL C .. C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS,MAX,MIN,SIGN C .. C ******************************************************************* C REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE C INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL C GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING C PRECISION ON COMPUTERS WITH LONG WORDLENGTHS. DATA REMIN /3.E-11/ C ******************************************************************* C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE. C AS SET,THIS CORRESPONDS TO ABOUT 500 STEPS. DATA MAXNFE /9000/ C ******************************************************************* C U - THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE C VALUE REPRESENTABLE IN THE MACHINE SUCH THAT 1.+ U .GT. 1. C (VARIABLE ONE SET TO 1.0 EASES PRECISION CONVERSION.) C ONE = 1.0 U = EPSLON(ONE) C ******************************************************************* C CHECK INPUT PARAMETERS IF (NEQN.LT.1) GO TO 10 IF ((RELERR.LT.0.) .OR. (ABSERR.LT.0.)) GO TO 10 MFLAG = ABS(IFLAG) IF ((MFLAG.GE.1) .AND. (MFLAG.LE.7)) GO TO 20 C INVALID INPUT 10 IFLAG = 7 RETURN C IS THIS THE FIRST CALL 20 IF (MFLAG.EQ.1) GO TO 70 C CHECK CONTINUATION POSSIBILITIES IF (T.EQ.TOUT) GO TO 10 IF (MFLAG.NE.2) GO TO 30 C IFLAG = +2 OR -2 IF (INIT.EQ.0) GO TO 60 IF (KFLAG.EQ.3) GO TO 50 IF ((KFLAG.EQ.4) .AND. (ABSERR.EQ.0.)) GO TO 40 IF ((KFLAG.EQ.5) .AND. (RELERR.LE.SAVRE) .AND. * (ABSERR.LE.SAVAE)) GO TO 40 GO TO 70 C IFLAG = 3,4,5,6, OR 7 30 IF (IFLAG.EQ.3) GO TO 50 IF ((IFLAG.EQ.4) .AND. (ABSERR.GT.0.)) GO TO 60 C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO C THE INSTRUCTIONS PERTAINING TO IFLAG=4,5,6 OR 7 40 STOP C ******************************************************************* C RESET FUNCTION EVALUATION COUNTER 50 NFE = 0 IF (MFLAG.EQ.2) GO TO 70 C RESET FLAG VALUE FROM PREVIOUS CALL 60 IFLAG = JFLAG C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT C INPUT CHECKING 70 JFLAG = IFLAG KFLAG = 0 C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS SAVRE = RELERR SAVAE = ABSERR C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS C 32U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING C FROM IMPOSSIBLE ACCURACY REQUESTS RER = MAX(RELERR,32.*U+REMIN) U26 = 26.*U DT = TOUT - T IF (MFLAG.EQ.1) GO TO 80 IF (INIT.EQ.0) GO TO 90 GO TO 110 C ******************************************************************* C INITIALIZATION -- C SET INITIALIZATION COMPLETION INDICATOR,INIT C SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP C EVALUATE INITIAL DERIVATIVES C COPY INITIAL VALUES AND DERIVATIVES FOR THE C PARALLEL SOLUTION C SET COUNTER FOR FUNCTION EVALUATIONS,NFE C ESTIMATE STARTING STEPSIZE 80 INIT = 0 KOP = 0 A = T CALL F(A, Y, YP) NFE = 1 IF (T.NE.TOUT) GO TO 90 IFLAG = 2 RETURN 90 INIT = 1 H = ABS(DT) TOLN = 0. DO 100 K=1,NEQN YG(K) = Y(K) YGP(K) = YP(K) TOL = RER*ABS(Y(K)) + ABSERR IF (TOL.LE.0.) GO TO 100 TOLN = TOL YPK = ABS(YP(K)) IF (YPK*H**5.GT.TOL) H = (TOL/YPK)**0.2 100 CONTINUE IF (TOLN.LE.0.) H = 0. H = MAX(H,U26*MAX(ABS(T),ABS(DT))) C ******************************************************************* C SET STEPSIZE FOR INTEGRATION IN THE DIRECTION FROM T TO TOUT 110 H = SIGN(H,DT) C TEST TO SEE IF GERK IS BEING SEVERELY IMPACTED BY TOO MANY C OUTPUT POINTS IF (ABS(H).GT.2.*ABS(DT)) KOP = KOP + 1 IF (KOP.NE.100) GO TO 120 KOP = 0 IFLAG = 6 RETURN 120 IF (ABS(DT).GT.U26*ABS(T)) GO TO 140 C IF TOO CLOSE TO OUTPUT POINT,EXTRAPOLATE AND RETURN DO 130 K=1,NEQN YG(K) = YG(K) + DT*YGP(K) Y(K) = Y(K) + DT*YP(K) 130 CONTINUE A = TOUT CALL F(A, YG, YGP) CALL F(A, Y, YP) NFE = NFE + 2 GO TO 230 C INITIALIZE OUTPUT POINT INDICATOR 140 OUTPUT = .FALSE. C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION, C SCALE THE ERROR TOLERANCES SCALE = 2./RER AE = SCALE*ABSERR C ******************************************************************* C ******************************************************************* C STEP BY STEP INTEGRATION 150 HFAILD = .FALSE. C SET SMALLEST ALLOWABLE STEPSIZE HMIN = U26*ABS(T) C ADJUST STEPSIZE IF NECESSARY TO HIT THE OUTPUT POINT. C LOOK AHEAD TWO STEPS TO AVOID DRASTIC CHANGES IN THE STEPSIZE C AND THUS LESSEN THE IMPACT OF OUTPUT POINTS ON THE CODE. DT = TOUT - T IF (ABS(DT).GE.2.*ABS(H)) GO TO 170 IF (ABS(DT).GT.ABS(H)) GO TO 160 C THE NEXT SUCCESSFUL STEP WILL COMPLETE THE INTEGRATION TO THE C OUTPUT POINT OUTPUT = .TRUE. H = DT GO TO 170 160 H = 0.5*DT C ******************************************************************* C CORE INTEGRATOR FOR TAKING A SINGLE STEP C ******************************************************************* C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW C IN COMPUTING THE ERROR TOLERANCE FUNCTION ET. C TO AVOID PROBLEMS WITH ZERO CROSSINGS, RELATIVE ERROR IS C MEASURED USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION C AT THE BEGINNING AND END OF A STEP. C THE ERROR ESTIMATE FORMULA HAS BEEN GROUPED TO CONTROL LOSS OF C SIGNIFICANCE. C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T. C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES. C TO PREVENT UNNECESSARY FAILURES, THE CODE USES 9/10 THE C STEPSIZE IT ESTIMATES WILL SUCCEED. C AFTER A STEP FAILURE, THE STEPSIZE IS NOT ALLOWED TO INCREASE C FOR THE NEXT ATTEMPTED STEP. THIS MAKES THE CODE MORE C EFFICIENT ON PROBLEMS HAVING DISCONTINUITIES AND MORE C EFFECTIVE IN GENERAL SINCE LOCAL EXTRAPOLATION IS BEING USED C AND THE ERROR ESTIMATE MAY BE UNRELIABLE OR UNACCEPTABLE WHEN C A STEP FAILS. C ******************************************************************* C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS. C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H 170 IF (NFE.LE.MAXNFE) GO TO 180 C TOO MUCH WORK IFLAG = 3 KFLAG = 3 RETURN C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H 180 CALL FEHL(F, NEQN, YG, T, H, YGP, F1, F2, F3, F4, F5, F1) NFE = NFE + 5 C COMPUTE AND TEST ALLOWABLE TOLERANCES VERSUS LOCAL ERROR C ESTIMATES AND REMOVE SCALING OF TOLERANCES. NOTE THAT RELATIVE C ERROR IS MEASURED WITH RESPECT TO THE AVERAGE MAGNITUDES OF THE C OF THE SOLUTION AT THE BEGINNING AND END OF THE STEP. EEOET = 0. DO 200 K=1,NEQN ET = ABS(YG(K)) + ABS(F1(K)) + AE IF (ET.GT.0.) GO TO 190 C INAPPROPRIATE ERROR TOLERANCE IFLAG = 4 KFLAG = 4 RETURN 190 EE = ABS((-2090.*YGP(K)+(21970.*F3(K)-15048.*F4(K))) * +(22528.*F2(K)-27360.*F5(K))) EEOET = MAX(EEOET,EE/ET) 200 CONTINUE ESTTOL = ABS(H)*EEOET*SCALE/752400. IF (ESTTOL.LE.1.) GO TO 210 C UNSUCCESSFUL STEP C REDUCE THE STEPSIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR OF 1/10 HFAILD = .TRUE. OUTPUT = .FALSE. S = 0.1 IF (ESTTOL.LT.59049.) S = 0.9/ESTTOL**0.2 H = S*H IF (ABS(H).GT.HMIN) GO TO 170 C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE IFLAG = 5 KFLAG = 5 RETURN C SUCCESSFUL STEP C STORE ONE-STEP SOLUTION YG AT T+H C AND EVALUATE DERIVATIVES THERE 210 TS = T T = T + H DO 220 K=1,NEQN YG(K) = F1(K) 220 CONTINUE A = T CALL F(A, YG, YGP) NFE = NFE + 1 C NOW ADVANCE THE Y SOLUTION OVER TWO STEPS OF C LENGTH H/2 AND EVALUATE DERIVATIVES THERE HH = 0.5*H CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) TS = TS + HH A = TS CALL F(A, Y, YP) CALL FEHL(F, NEQN, Y, TS, HH, YP, F1, F2, F3, F4, F5, Y) A = T CALL F(A, Y, YP) NFE = NFE + 12 C CHOOSE NEXT STEPSIZE C THE INCREASE IS LIMITED TO A FACTOR OF 5 C IF STEP FAILURE HAS JUST OCCURRED, NEXT C STEPSIZE IS NOT ALLOWED TO INCREASE S = 5. IF (ESTTOL.GT.1.889568E-4) S = 0.9/ESTTOL**0.2 IF (HFAILD) S = MIN(S,ONE) H = SIGN(MAX(S*ABS(H),HMIN),H) C ******************************************************************* C END OF CORE INTEGRATOR C ******************************************************************* C SHOULD WE TAKE ANOTHER STEP IF (OUTPUT) GO TO 230 IF (IFLAG.GT.0) GO TO 150 C ******************************************************************* C ******************************************************************* C INTEGRATION SUCCESSFULLY COMPLETED C ONE-STEP MODE IFLAG = -2 GO TO 240 C INTERVAL MODE 230 T = TOUT IFLAG = 2 240 DO 250 K=1,NEQN GERROR(K) = (YG(K)-Y(K))/31. 250 CONTINUE RETURN END SUBROUTINE FEHL(F, NEQN, Y, T, H, YP, F1, F2, F3, F4, F5, S) C FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD C ******************************************************************* C FEHL INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT=F(T,Y(1),---,Y(NEQN)) C WHERE THE INITIAL VALUES Y(I) AND THE INITIAL DERIVATIVES C YP(I) ARE SPECIFIED AT THE STARTING POINT T. FEHL ADVANCES C THE SOLUTION OVER THE FIXED STEP H AND RETURNS C THE FIFTH ORDER (SIXTH ORDER ACCURATE LOCALLY) SOLUTION C APPROXIMATION AT T+H IN ARRAY S(I). C F1,---,F5 ARE ARRAYS OF DIMENSION NEQN WHICH ARE NEEDED C FOR INTERNAL STORAGE. C THE FORMULAS HAVE BEEN GROUPED TO CONTROL LOSS OF SIGNIFICANCE. C FEHL SHOULD BE CALLED WITH AN H NOT SMALLER THAN 13 UNITS OF C ROUNDOFF IN T SO THAT THE VARIOUS INDEPENDENT ARGUMENTS CAN BE C DISTINGUISHED. C ******************************************************************* C .. SCALAR ARGUMENTS .. INTEGER NEQN DOUBLE PRECISION H,T C .. C .. ARRAY ARGUMENTS .. DOUBLE PRECISION F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN), 1 S(NEQN),Y(NEQN),YP(NEQN) C .. C .. SUBROUTINE ARGUMENTS .. EXTERNAL F C .. C .. LOCAL SCALARS .. INTEGER K DOUBLE PRECISION CH C .. CH = 0.25*H DO 10 K=1,NEQN F5(K) = Y(K) + CH*YP(K) 10 CONTINUE CALL F(T+0.25*H, F5, F1) CH = 0.09375*H DO 20 K=1,NEQN F5(K) = Y(K) + CH*(YP(K)+3.*F1(K)) 20 CONTINUE CALL F(T+0.375*H, F5, F2) CH = H/2197. DO 30 K=1,NEQN F5(K) = Y(K) + CH*(1932.*YP(K)+(7296.*F2(K)-7200.*F1(K))) 30 CONTINUE CALL F(T+12./13.*H, F5, F3) CH = H/4104. DO 40 K=1,NEQN F5(K) = Y(K) + CH*((8341.*YP(K)-845.*F3(K))+(29440.*F2(K) * -32832.*F1(K))) 40 CONTINUE CALL F(T+H, F5, F4) CH = H/20520. DO 50 K=1,NEQN F1(K) = Y(K) + CH*((-6080.*YP(K)+(9295.*F3(K)-5643.*F4(K))) * +(41040.*F1(K)-28352.*F2(K))) 50 CONTINUE CALL F(T+0.5*H, F1, F5) C COMPUTE APPROXIMATE SOLUTION AT T+H CH = H/7618050. DO 60 K=1,NEQN S(K) = Y(K) + CH*((902880.*YP(K)+(3855735.*F3(K)-1371249.* * F4(K)))+(3953664.*F2(K)+277020.*F5(K))) 60 CONTINUE RETURN END DOUBLE PRECISION FUNCTION EPSLON (X) DOUBLE PRECISION X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C DOUBLE PRECISION A,B,C,EPS,FOUR,THREE C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C FOUR = 4.0 THREE = 3.0 A = FOUR/THREE 10 B = A - 1.0 C = B + B + B EPS = ABS(C-1.0) IF (EPS .EQ. 0.0) GO TO 10 EPSLON = EPS*ABS(X) RETURN END