C Routines for 2nd order non-selfadjoint problems SUBROUTINE SLMISS(XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS,F,CONST, & IFLAG) IMPLICIT NONE INTEGER NMESH,IMATCH,IFLAG DOUBLE PRECISION XMESH(0:NMESH),CONST DOUBLE COMPLEX ELAM,F EXTERNAL SL2COF,SL2BCS C This routine computes the miss-distance function for the second order C equation C C -(p(x)y')' + q_1(x) y' + q_0(x) y = 0 C C Here any of the coefficients is (in principle) allowed to depend on the C eigenparameter lambda = ELAM. C C The equation is cast in the quasi-Hamiltonian form C C J z' = S z C C where J = /0 -1\ and S = /-q_0(x) -q_1(x)/p(x)\ C \1 0/ \ 0 1/p(x) / C C and z = (y,py')^T. C C For this problem the Riccati-like variables are y and py', so the C Riccati-like system is just the same as the system for z. C INTEGER IEND DOUBLE COMPLEX ZL(2),ZR(2),Y,PDY DOUBLE PRECISION CONSTL, CONSTR C Get boundary conditions at the left hand endpoint: IEND = 0 CALL SL2BCS(IEND,XMESH(0),ELAM,Y,PDY) ZL(1) = Y ZL(2) = PDY C WRITE (7,FMT=*) 'Left BCS Y, PDY:',Y,PDY C Integrate from left endpoint to centre: CALL EXPINT(ZL,CONSTL,XMESH,NMESH,0,IMATCH,SL2COF,ELAM,IFLAG) IF (IFLAG.NE.0) RETURN C WRITE (7,FMT=*) 'At matchpoint, ZL, CONST:',ZL(1),ZL(2),CONSTL C Get boundary conditions at the right hand endpoint: IEND = 1 CALL SL2BCS(IEND,XMESH(NMESH),ELAM,Y,PDY) ZR(1) = Y ZR(2) = PDY C WRITE (7,FMT=*) 'Right BCS Y, PDY:',Y,PDY CALL EXPINT(ZR,CONSTR,XMESH,NMESH,NMESH,IMATCH,SL2COF,ELAM,IFLAG) IF (IFLAG.NE.0) RETURN C WRITE (7,FMT=*) 'At matchpoint, ZR, CONST:',ZR(1),ZR(2),CONSTL C Miss distance is the determinant of the matrix C C / yl yr \ C \ pyl' pyr'/ C C F = ZL(1)*ZR(2)-ZL(2)*ZR(1) C WRITE (7,FMT=*) 'ELAM,F,CONST:',ELAM,F,CONSTL+CONSTR CONST = CONSTL+CONSTR RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPINT(Z,CONST,XMESH,NMESH,ISTART,IEND,SL2COF,ELAM, & IFAIL) C C This routine integrates the linear system for the 2nd order C Riccati variables across a mesh in each interval of which the C coefficients in the 2nd order ODE are approximated using the C Magnus method (and are therefore treated as constant). C C Z: on entry, the array of Riccati variables at XMESH(ISTART) C on exit, the Riccati variables at XMESH(IEND) scaled by C exp(-CONST) to avoid overflow. C C CONST: on entry, CONST = 0. On exit, CONST specifies the number C such that exp(CONST).Z recovers the Riccati variables. N.B. C computation of exp(CONST) may cause overflow. C C XMESH: input. The mesh, on each interval of which the coefficients C in the ODE are assumed to be constant. C C SL2COF: specifies the coefficients in the ODE, which is C C -(P y')' + Q1 y' + Q0 y = 0. C C The dependence on the eigenparameter ELAM is hidden in the C coefficients. C C IFAIL: usual NAG-type error flag. C IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST INTEGER IEND,IFAIL,ISTART,NMESH C .. C .. Array Arguments .. DOUBLE COMPLEX Z(2) DOUBLE PRECISION XMESH(0:NMESH) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE COMPLEX CONE,CZERO DOUBLE PRECISION CNLOC,ZLMAX INTEGER I,IDI,II C .. C .. Local Arrays .. DOUBLE COMPLEX F(2,2),ZLOC(2) C .. C .. External Subroutines .. EXTERNAL EXPSOL,F06SAF,F06TFF C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,LOG,MAX,MIN C .. CONE = CMPLX(1.D0,0.D0) CZERO = CMPLX(0.D0,0.D0) IF (MIN(ISTART,IEND).LT.0 .OR. MAX(ISTART,IEND).GT.NMESH) THEN C Cannot integrate over this range IFAIL = 1 RETURN END IF CONST = 0.D0 IDI = 1 IF (IEND.LT.ISTART) IDI = -1 IF (ISTART.EQ.IEND) RETURN C WRITE (7,FMT=*) '--- ELAM = ',ELAM DO 30 I = ISTART + IDI,IEND,IDI CNLOC = 0.D0 IFAIL = 0 CALL EXPSOL(SL2COF,ELAM,XMESH(I-IDI),XMESH(I),CNLOC,F,2,IFAIL) IF (IFAIL.NE.0) THEN IFAIL = 2 RETURN END IF CONST = CONST + CNLOC C WRITE (7,FMT=*) 'Here is the fundamental matrix:' C WRITE (7,FMT=*) F(1,1),F(1,2) C WRITE (7,FMT=*) F(2,1),F(2,2) CALL F06SAF('N',2,2,CONE,F,2,Z,1,CZERO,ZLOC,1) CALL F06TFF('G',2,1,ZLOC,2,Z,2) C Normalise to get the Riccati variables to a reasonable scale ZLMAX = 0.D0 DO 10 II = 1,2 ZLMAX = MAX(ZLMAX,ABS(Z(II))) 10 CONTINUE DO 20 II = 1,2 Z(II) = Z(II)/ZLMAX 20 CONTINUE CONST = CONST + LOG(ZLMAX) C WRITE (7,FMT=*) 'X,Z(1),Z(2),CONST:',XMESH(I),Z(1),Z(2),CONST 30 CONTINUE IFAIL = 0 RETURN END C ------------------------------------------------------------------- SUBROUTINE EXPSOL(SL2COF,ELAM,XSTART,XEND,CONST,F,IF,IFAIL) C C The Wronskians in the 2nd order system satisfy a linear system of C the form z' = Az, where A is a 2x2 matrix. Here we C attempt to solve this system. C C The equation being solved is C C -(P y')' + Q1 y' + Q0 y = 0 C C Routine to compute exp(A*(xend-xstart)) using routine from Golub and Van Loan C The output is exp(const).F; const is a real constant, a scaling used C to avoid overflow. C IMPLICIT NONE C .. Parameters .. INTEGER IX,IXOUT,ID,IN PARAMETER (IX=2,IXOUT=IX,ID=IX,IN=IX) C .. C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION CONST,XEND,XSTART INTEGER IF,IFAIL C .. C .. Array Arguments .. DOUBLE COMPLEX F(IF,2) C .. C .. Local Scalars .. DOUBLE COMPLEX DIAG,EQR,FACC,MIN1K,OFDIAG,TERM DOUBLE PRECISION DUMMY,FAC,ERREST INTEGER II,IJ,IQ,J,K C .. C .. Local Arrays .. DOUBLE COMPLEX D(ID,IX),N(IN,IX),X(IX,IX),X1(IX,IX),XOUT(IXOUT,IX) DOUBLE PRECISION WREAL(IX) C .. C .. External Functions .. DOUBLE PRECISION F06UAF,X02AJF EXTERNAL F06UAF,X02AJF C .. C .. External Subroutines .. EXTERNAL F04ADF,F06TFF,F06THF,F06ZAF,GETSIG C .. C .. Intrinsic Functions .. INTRINSIC ABS,CMPLX,DBLE,INT,LOG,MAX C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. OFDIAG = CMPLX(0.D0,0.D0) DIAG = CMPLX(1.D0,0.D0) CALL F06THF('G',2,2,OFDIAG,DIAG,X,IX) CALL GETSIG(SL2COF,ELAM,XSTART,XEND,XOUT,IXOUT,ERREST) C XOUT contains the matrix A (not in sparse form). F06UAF computes its C sup norm. C F06UAF here gets the (infinity) norm of a matrix J = MAX(0,1+INT(LOG(F06UAF('I',2,2,XOUT,IXOUT,WREAL)/LOG(2.D0)))) FAC = 0.5D0**J FACC = CMPLX(FAC,0.D0) C WRITE (8,FMT=*) 'J,FAC:',J,FAC IQ = 1 EQR = CMPLX(1.D0/6.D0,0.D0) 10 EQR = EQR/CMPLX(DBLE(16*IQ*IQ* (2*IQ+1)* (2*IQ+3)),0.D0) C WRITE (8,FMT=*) 'E_q:',EQR IQ = IQ + 1 IF (ABS(EQR).GT.X02AJF(1.D0)) GO TO 10 C WRITE (6,FMT=*) 'Q is',IQ C Initialise matrices D, N and X all to be equal to the identity CALL F06THF('G',2,2,OFDIAG,DIAG,X,IX) CALL F06THF('G',2,2,OFDIAG,DIAG,N,IN) CALL F06THF('G',2,2,OFDIAG,DIAG,D,ID) CONST = 1.D0 DO 40 K = 1,IQ CONST = CONST* (DBLE(IQ-K+1)/DBLE(IQ+IQ-K+1))/DBLE(K) CALL F06ZAF('N','N',2,2,2,FACC,XOUT,IXOUT,X,IX,OFDIAG,X1,IX) CALL F06TFF('G',2,2,X1,IX,X,IX) MIN1K = (-1.D0)**K DO 30 IJ = 1,2 DO 20 II = 1,2 TERM = CONST*X(II,IJ) N(II,IJ) = N(II,IJ) + TERM D(II,IJ) = D(II,IJ) + MIN1K*TERM 20 CONTINUE 30 CONTINUE 40 CONTINUE C Solve DF = N for F using a suitable routine IFAIL = 0 C CALL F04ADF(D,ID,N,IN,2,2,F,IF,WREAL,IFAIL) C For k = 1 to j do F = F*F. N.B. F is almost certainly C full. C Normalise F: CONST = MAX(1.D0,F06UAF('M',2,2,F,IF,WREAL)) DO 60 IJ = 1,2 DO 50 II = 1,2 F(II,IJ) = F(II,IJ)/CONST 50 CONTINUE 60 CONTINUE CALL F06TFF('G',2,2,F,IF,X,IX) CONST = LOG(CONST) C WRITE (8,FMT=*) 'CONST (1):',DBLE(2**J)*CONST DO 90 K = 1,J CALL F06ZAF('N','N',2,2,2,DIAG,F,IF,X,IX,OFDIAG,XOUT,IXOUT) CALL F06TFF('G',2,2,XOUT,IXOUT,F,IF) CALL F06TFF('G',2,2,XOUT,IXOUT,X,IX) CONST = CONST + CONST DUMMY = F06UAF('M',2,2,F,IF,WREAL) C WRITE (8,FMT=*) 'DUMMY=',DUMMY DO 80 IJ = 1,2 DO 70 II = 1,2 F(II,IJ) = F(II,IJ)/DUMMY X(II,IJ) = X(II,IJ)/DUMMY 70 CONTINUE 80 CONTINUE CONST = CONST + LOG(DUMMY) 90 CONTINUE C F now contains an approximation to exp(A(xend-xstart))/exp(CONST). 9000 FORMAT (A) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETSIG(SL2COF,ELAM,XSTART,XEND,SIGMA,ISIGMA,ERREST) IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION XEND,XSTART,ERREST INTEGER ISIGMA C .. C .. Array Arguments .. DOUBLE COMPLEX SIGMA(ISIGMA,2) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE PRECISION H,T1,T2 INTEGER I,J C .. C .. Local Arrays .. DOUBLE COMPLEX A1(2,2),A2(2,2),ADIFF(2,2),C1(2,2),C2(2,2) C .. C .. Parameters .. DOUBLE PRECISION Q1,Q2,W1,W2 PARAMETER (Q1=0.21132486540519D0,Q2=0.78867513459481D0, & W1=0.14433756729741D0,W2=1.D0/80.D0) C .. C .. External Subroutines .. EXTERNAL GETA,LIEB C .. T1 = XSTART + Q1* (XEND-XSTART) T2 = XSTART + Q2* (XEND-XSTART) H = XEND - XSTART CALL GETA(SL2COF,ELAM,T1,A1,2) CALL GETA(SL2COF,ELAM,T2,A2,2) CALL LIEB(A2,2,A1,2,C1,2,2) DO 20 J = 1,2 DO 10 I = 1,2 ADIFF(I,J) = A2(I,J) - A1(I,J) 10 CONTINUE 20 CONTINUE CALL LIEB(ADIFF,2,C1,2,C2,2,2) ERREST = 0.D0 DO 40 J = 1,2 DO 30 I = 1,2 SIGMA(I,J) = H* ((A1(I,J)+A2(I,J))/2.D0+ & H* (W1*C1(I,J)+H*W2*C2(I,J))) ERREST = MAX(ERREST,ABS(W2*C2(I,J))) 30 CONTINUE 40 CONTINUE C H = H/MAX(1.D0,MAX(XSTART**2,XEND**2)) C ERREST = ABS(ERREST*H*H*H) C H = H/MAX(1.D0,MAX(ABS(XSTART),ABS(XEND))**0.25) C ERREST = ABS(ERREST*H*H*H) C H = H/MAX(ABS(H),1.D0) C ERREST = ABS(ERREST*H*H*H) C H = H/MAX(1.D0,MAX(ABS(XSTART),ABS(XEND))) ERREST = ABS(ERREST*H) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETMSH(A,B,XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS, & TOL,IFAIL) IMPLICIT NONE C C This subroutine uses the natural error estimate from the Magnus formula C as a means of automatic mesh generation. C INTEGER NMESH,IMATCH,IFAIL,I,ISIG PARAMETER (ISIG=2) DOUBLE PRECISION A,B,XMESH(0:NMESH),TOL,H,ERREST,ERRMIN DOUBLE COMPLEX ELAM,SIGMA(ISIG,ISIG) EXTERNAL SL2COF,SL2BCS C The direction of integration is not important for setting up the mesh, so C we do it in one big sweep from x=A to x=B. C Set the initial stepsize: IF (B.LE.A) THEN IFAIL = -1 RETURN END IF H = (B-A)/DBLE(MAX(NMESH,100)) I = 0 XMESH(I) = A 100 IF ((I+1).GT.NMESH) THEN IFAIL = 2 RETURN END IF XMESH(I+1) = XMESH(I) + H C & MIN(ABS(H),0.05D0)*SIGN(1.D0,B-A) C WRITE (17,FMT=*) 'XMESH(',I+1,') = ',XMESH(I+1) IF (XMESH(I+1).GT.B) XMESH(I+1) = B C Get error estimate C WRITE (17,FMT=*) 'Calling GETSIG' CALL GETSIG(SL2COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,ISIG,ERREST) C WRITE (17,FMT=*) 'Called GETSIG, ERREST=',ERREST IF (ABS(ERREST).LE.TOL) THEN C The step was successful. We can go to the next step, but first we C may want to increase the stepsize. IF (2.D0*ABS(ERREST).LT.TOL) THEN H = H/MAX(0.2D0,(2.D0*ABS(ERREST)/TOL)**0.2D0) END IF I = I + 1 IF (XMESH(I).NE.B) GOTO 100 C If we are here it means we have finished meshing; re-set the value of C NMESH. NMESH = I ELSE C The step was unsuccessful. We must reduce the stepsize and repeat C the step. H = H*MIN(0.8D0,(TOL/ABS(ERREST))**0.2D0) IF (ABS(H).LT.1.D-8) THEN C The stepsize is too small to make further progress IFAIL = 1 RETURN END IF GOTO 100 END IF C We have now finished meshing. We have re-set the value of NMESH. We must C also choose the best matching point. Choose the point where there is `the C least going on'. IMATCH = 0 ERRMIN = 1.0D38 DO 10 I = 1,NMESH-1 CALL GETSIG(SL2COF,ELAM,XMESH(I),XMESH(I+1),SIGMA,2,ERREST) H = XMESH(I+1)-XMESH(I) H = H/MAX(1.D0,MAX(XMESH(I)**2,XMESH(I+1)**2)) ERREST = ERREST/(H**3) IF (ERREST.LT.ERRMIN) IMATCH = I ERRMIN = MIN(ERRMIN,ERREST) 10 CONTINUE C The matching point has now been reset. IFAIL = 0 RETURN END C --------------------------------------------------------------------------- SUBROUTINE LIEB(A,IA,B,IB,C,IC,N) IMPLICIT NONE C Compute AB-BA and put the answer in C C .. Parameters .. DOUBLE COMPLEX CONE,CMONE,CZERO PARAMETER (CONE= (1.D0,0.D0),CMONE= (-1.D0,0.D0), & CZERO= (0.D0,0.D0)) C .. C .. Scalar Arguments .. INTEGER IA,IB,IC,N C .. C .. Array Arguments .. DOUBLE COMPLEX A(IA,N),B(IB,N),C(IC,N) C .. C .. External Subroutines .. EXTERNAL F06ZAF C .. CALL F06ZAF('N','N',N,N,N,CONE,A,IA,B,IB,CZERO,C,IC) CALL F06ZAF('N','N',N,N,N,CMONE,B,IB,A,IA,CONE,C,IC) RETURN END C --------------------------------------------------------------------------- SUBROUTINE GETA(SL2COF,ELAM,X,A,IA) IMPLICIT NONE C .. Scalar Arguments .. DOUBLE COMPLEX ELAM DOUBLE PRECISION X INTEGER IA C .. C .. Array Arguments .. DOUBLE COMPLEX A(IA,2) C .. C .. Subroutine Arguments .. EXTERNAL SL2COF C .. C .. Local Scalars .. DOUBLE COMPLEX Q0,Q1 DOUBLE PRECISION P INTEGER I,J C .. CALL SL2COF(X,ELAM,P,Q1,Q0) A(1,1) = CMPLX(0.D0,0.D0) A(2,1) = Q0 A(1,2) = CMPLX(1.D0/P,0.D0) A(2,2) = Q1/P RETURN END