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
.