C ALGORITHM 690, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 17, NO. 2, PP. 178-206. JUNE, 1991. c c file containg example programs and PDECHEB software. c c This file contains c example problem 1 c example problem 2 c example problem 3 c example problem 4 c example problem 5 c dassl time integration routine c pdecheb spatial discretisation routine c interface to AND the linpack full and banded routines c c C C0 COLLOCATION PARAMETERS PARAMETER ( IBK = 21, NEL = IBK-1 , NPDE = 1, NV = 1, 1 NPOLY = 2, NPTS = NEL*NPOLY+1, NXI = 1, 2 NEQ = NPTS * NPDE + NV, 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) + 4 NPDE * 8 + 3 + NV + NXI, C DDASSL TIME INTEGRATION PARAMETERS 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2, 6 LIW = 20 + NEQ ) C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, IRESWK, 1 IDEV, ITRACE DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ), 1 WKRES(NWKRES), RWORK(LRW), XI(1), T, TOUT, RTOL, ATOL, 2 ENORM, GERR, VERROR, CTIME, TOL EXTERNAL PDECHB, DGEJAC COMMON /SDEV2/ ITRACE, IDEV COMMON /PROB1/ TOL TOL = 0.1D-5/50.D0 C N.B. CPU TIMER COMMENTED OUT FOR PORTABILITY C CALL TIMER( CTIME, 1) M = 0 T = TOL IDEV = 4 ITRACE = 1 WRITE(IDEV,9)NPOLY, NEL 9 FORMAT(' TEST PROBLEM 1'/' ***********'/' POLY OF DEGREE =',I4, 1 ' NO OF ELEMENTS = ',I4) XI(1) = 1.0D0 DO 10 I = 1,IBK 10 XBK(I) = (I-1.0D0)/(IBK-1.0D0) C INITIALISE THE P.D.E. WORKSPACE ITIME = 1 CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) IF(ITIME .EQ. -1)THEN WRITE(IDEV, 15) 15 FORMAT(' INITCC ROUTINE RETURNED ITIME = -1 - RUN HALTED ') GOTO 100 END IF C SETUP DASSL PARAMETERS RTOL = TOL ATOL = TOL DO 20 I = 1,11 20 INFO(I) = 0 C C BANDED MATRIX OPTION WHEN INFO(6) = 1 IF(INFO(6) .EQ. 1)THEN IWORK(1) = IBAND IWORK(2) = IBAND END IF 30 TOUT = T * 10.0D0 IF(TOUT .GE. 2.D0)TOUT =2.0D0 CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL, 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC) IF( IDID .LT. 0 )THEN C DASSL FAILED TO FINISH INTEGRATION. WRITE(IDEV,40)T,IDID 40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3) GOTO 100 ELSE C DASSL INTEGRATED TO T = TOUT C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION. ITRACE = 1 CALL ERROR( Y, NPDE, NPTS, X, M, ENORM, GERR, T, RTOL, ATOL, 1 ITRACE, WKRES, NWKRES) ITRACE = 0 VERROR = Y(NEQ) - T WRITE(IDEV,50)Y(NEQ),VERROR 50 FORMAT(' MOVING BOUNDARY IS AT ',D12.4,' WITH ERROR=',D12.4) IF(TOUT .LT. 1.99D0 ) GOTO 30 END IF 100 CONTINUE C CALL TIMER(CTIME, 2) WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13), CTIME 110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4,' CPU=',D11.3) STOP END C********************************************************************** C EXAMPLE PROBLEM 1 C SOLUTION OF MOVING BOUNDARY PROBLEM BY CO-ORDINATE TRANSFORMATION. C******************************************************************** C THIS PROBLEM IS THE ONE PHASE STEFAN PROBLEM (HOFFMAN (1977) ) SEE C FURZELAND R.M. A COMPARATIVE STUDY OF NUMERICAL METHODS FOR MOVING C BOUNDARY PROBLEMS. J.I.M.A. (1977) ,26, PP 411 - 429. C THE PROBLEM HAS MELTING DUE TO HEAT INPUT AT THE FIXED C BOUNDARY . THE P.D.E. IS DEFINED BY THE EQUATIONS C U = U 0 < Y < S(T) , 0.1 < T < 1 C T YY C U = - EXP(T) , Y = 0 C Y . C U = 0 AND S(T) = - U ON THE MOVING BOUNDARY Y = S(T). C Y C AND THE INITIAL SOLUTION VALUES AT T = 0.1 ARE GIVEN BY THE ANALYTIC C SOLUTION C U = EXP(T-Y) - 1 , S(T) = T. C THE PROBLEM IS REWRITTEN BY USING THE CO-ORDINATE TRANSFORMATION C GIVEN BY X(T) = Y / S(T) . THE EQUATIONS THEN READ C . C S * S * U - S * S * X * U = U , X IN (0,1). C T X XX C WITH THE NEUMANN TYPE BOUNDARY CONDITIONS C . C U = - EXP(T) AT X=0 AND U = - S(T) * S(T) AT X = 1 C X X C AND THE O.D.E. COUPLING POINT EQUATION AT X = 1 WHICH IMPLICITLY C DEFINES S(T) IS GIVEN BY C U(1,T) = 0 C THE EXACT SOLUTION IS NOW DEFINED BY C U(X,T) = EXP((T - X*S(T)) , S(T) = T C C WE SHALL NOW DETAIL THE ROUTINES NEEDED TO DESCRIBE THIS PROBLEM. C PROBLEM DESCRIPTION ROUTINES C ****************************** C EXACT SOLUTION SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U) C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN) INTEGER NPDE, NPTS DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME DO 10 I = 1,NPTS 10 U(1,I) = DEXP( TIME * (1 - X(I))) - 1.0D0 RETURN END SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV, V) C ROUTINE FOR O.D.E. AND P.D.E. INITIAL VALUES. INTEGER NPDE, NPTS, NV DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME, V(NV) COMMON /PROB1/ TOL TIME= TOL V(1)= TOL CALL EXACT(TIME,NPDE,NPTS,X,U) RETURN END C SUBROUTINE SPDEFN(T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R, 1 NV, V, VDOT, IRES) C PROBLEM INTERFACE FOR THE MOVING BOUNDARY PROBLEM. INTEGER NPTL, NPDE, NV, I, IRES DOUBLE PRECISION X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), T, 1 V(1), VDOT(1), Q(NPDE,NPTL) ,R(NPDE,NPTL), 2 UDOT(NPDE,NPTL), UTDX(NPDE,NPTL) DO 10 I = 1,NPTL R(1,I) = DUDX(1,I) Q(1,I) = V(1)*V(1)*UDOT(1,I) -X(I)*VDOT(1)*DUDX(1,I) * V(1) 10 CONTINUE RETURN END SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, 1 LEFT, NV, V, VDOT, IRES) C THIS SUBROUTINE PROVIDES THE LEFT AND RIGHT BOUNDARY VALUES C FOR THE MOVING BOUNDARY PROBLEM IN THE FORM. C BETA(I) * DU/DX(I) = GAMMA(I) C WHERE I = 1,NPDE AND GAMMA IS A FUNCTION OF U,X AND T C INTEGER NPDE, NV, IRES LOGICAL LEFT DOUBLE PRECISION BETA(NPDE), GAMMA(NPDE), U(NPDE), UX(NPDE) - ,T, V(1), VDOT(1), UDOT(NPDE), UTDX(NPDE) BETA(1) = 1.0D0 IF(LEFT)THEN GAMMA(1) = -V(1)*DEXP(T) ELSE GAMMA(1) = -V(1)*VDOT(1) END IF RETURN END C SUBROUTINE SODEFN(T, NV, V, VDOT, NPDE, NXI, XI, UI, UXI, RI, 1 UTI, UTXI, VRES, IRES) C ROUTINE TO PROVIDE RESIDUAL OF COUPLED ODE SYSTEM FOR THE C MOVING BOUNDARY PROBLEM. C NOTE HOW IRES CAN BE RESET TO COPE WIH ILLEGAL VALUES OF THE C MOVING BOUNDARY POSITION V(1). INTEGER NPDE, NXI, NV, IRES DOUBLE PRECISION T, XI(NXI), UI(NPDE,NXI), UXI(NPDE,NXI), 1 RI(NPDE,NXI), UTI(NPDE,NXI), UTXI(NPDE,NXI), VRES(NV), 2 V(NV), VDOT(NV) VRES(1) = UI(1,1) IF(V(1) .LT. 0.0D0)IRES = -1 RETURN END C C0 COLLOCATION PARAMETERS PARAMETER ( IBK = 2, NEL = IBK-1 , NPDE = 1, NV = 0, 1 NPOLY = 10, NPTS = NEL*NPOLY+1, NXI = 0, 2 NEQ = NPTS * NPDE + NV, C C NWKRES= 2*(NPOLY+1)*(NPOLY+NEL+2) + 2 + NV + 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) + 4 NPDE * 8 + 3 + NV + NXI, C C NPDE * (7 * (NPOLY+1+NXI) + 8), C DDASSL TIME INTEGRATION PARAMETERS 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2, 6 LIW = 20 + NEQ ) C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, IRESWK, 1 IDEV, ITRACE DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ), 1 WKRES(NWKRES), RWORK(LRW), XI(1), T, TOUT, RTOL, ATOL, 2 ENORM, GERR, CTIME EXTERNAL PDECHB, DGEJAC COMMON /SDEV2/ ITRACE, IDEV C CPU TIMER COMMENTED OUT FOR PORTABILITY C CALL TIMER(CTIME ,1) M = 2 T = 0.0D0 IDEV = 4 ITRACE = 1 WRITE(IDEV,9)NPOLY, NEL 9 FORMAT(' TEST PROBLEM 1'/' ***********'/' POLY OF DEGREE =',I4, 1 ' NO OF ELEMENTS = ',I4) DO 10 I = 1,IBK 10 XBK(I) = (I-1.0D0) / (IBK - 1.0D0) C INITIALISE THE P.D.E. WORKSPACE ITIME = 1 CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) IF(ITIME .EQ. -1)THEN WRITE(IDEV, 15) 15 FORMAT(' INITCC ROUTINE RETURNED ITIME = -1 - RUN HALTED ') GOTO 100 END IF C SETUP DASSL PARAMETERS RTOL = 1.0D-8 ATOL = 1.0D-8 DO 20 I = 1,11 20 INFO(I) = 0 C INFO(11)= 1 C BANDED MATRIX OPTION WHEN INFO(6) = 1 IF(INFO(6) .EQ. 1)THEN IWORK(1) = IBAND IWORK(2) = IBAND END IF 30 TOUT = T + 0.1D0 CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL, 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC) IF( IDID .LT. 0 )THEN C DASSL FAILED TO FINISH INTEGRATION. WRITE(IDEV,40)T,IDID 40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3) GOTO 100 ELSE C DASSL INTEGRATED TO T = TOUT C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION. CALL ERROR( Y, NPDE, NPTS, X, M, ENORM, GERR, T, RTOL, ATOL, 1 ITRACE, WKRES, NWKRES) IF(TOUT .LT. 0.99D0 ) GOTO 30 END IF 100 CONTINUE C CALL TIMER(CTIME, 2) WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13), CTIME 110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4,' CPU=',D11.3) STOP END C EXAMPLE PROBLEM TWO C ******************** C THIS PROBLEM IS DEFINED BY C -2 2 2 C U U = X ( X U U ) + 5 U + 4 X U U , X IN (0,1) C T X X X C C THE LEFT BOUNDARY CONDITION AT X = 0 (LEFT = .TRUE. ) IS GIVEN BY C U (0,T) = 0.0 C X C THE RIGHT BOUNDARY CONDITION IS (LEFT = .FALSE.) C U( 1,T) = EXP ( -T ) C C THE INITIAL CONDITION IS GIVEN BY THE EXACT SOLUTION ; C U( X, T ) = EXP ( 1 - X*X - T ) , X IN ( 0,1) C 2 C********************************************************************** SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V) C ROUTINE FOR P.D.E. INITIAL VALUES. INTEGER NPDE, NPTS, NV DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), V(1), T T = 0.0D0 C V(1) IS A DUMMY VARIABLE AS THERE ARE NO COUPLED O.D.E.S CALL EXACT( T, NPDE, NPTS, X, U ) RETURN END C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R, 1 NV, V, VDOT, IRES) C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E. C THE P.D.E. IS WRITEN AS -M M C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U )) C X T TX X T TX X C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE. C DEFINITIONS FOR THE MODEL PROBLEM ARE GIVEN C NOTE NV = 0 : THERE IS NO O.D.E PART. INTEGER NPDE, NPTL, I, J, NV, IRES DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), 1 UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V, VDOT, 2 UTDX(NPDE,NPTL) DO 10 I = 1,NPTL R(1,I) = U(1,I) * DUDX(1,I) Q(1,I) = U(1,I) * UDOT(1,I) - 5.0D0 * U(1,I)**2 1 - 4.0D0 * U(1,I)*DUDX(1,I)*X(I) 10 CONTINUE RETURN END C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT, 1 NV, V, VDOT, IRES) C BOUNDARY CONDITIONS ROUTINE INTEGER NPDE, NV, IRES DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), C2, 1 UX(NPDE), V, VDOT, UDOT(NPDE), UTDX(NPDE) LOGICAL LEFT IF(LEFT) THEN BETA (1) = 1.0D0 GAMMA(1) = 0.0D0 ELSE C BETA (1) = 0.0D0 C GAMMA(1) = U(1) - DEXP( -T ) BETA (1) = 1.0D0 GAMMA(1) = - 2.D0 *U(1)**2 END IF RETURN END C C DUMMY O.D.E. ROUTINE AS NV IS ZERO SUBROUTINE SODEFN RETURN END C EXACT SOLUTION SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U) C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN) INTEGER NPDE, NPTS, I DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME DO 10 I = 1,NPTS 10 U(1,I) = DEXP( 1.0D0 - X(I)**2 - TIME) RETURN END c problem 3 C C0 COLLOCATION PARAMETERS PARAMETER ( IBK = 3, NEL = IBK-1 , NPDE = 1, NV = 0, 1 NPOLY = 6, NPTS = NEL*NPOLY+1, NXI = 0, 2 NEQ = NPTS * NPDE + NV, 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) + 4 NPDE * 8 + 3 + NV + NXI, C 3 NWKRES= 2*(NPOLY+1)*(NPOLY+NEL+2) + 2 + NV + C 4 NPDE * (7 * (NPOLY+1+NXI) + 8), C DDASSL TIME INTEGRATION PARAMETERS 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2, 6 LIW = 20 + NEQ ) C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, IRESWK, 1 IDEV, ITRACE, IDERIV, IFL, ITYPE DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ), Z(NPTS), 1 WKRES(NWKRES), RWORK(LRW), XI(1), T, TOUT, RTOL, ATOL, 2 ENORM, GERR, CTIME, DYDX(NEQ), DYCALC(NPDE,NPTS,2) EXTERNAL PDECHB, DGEJAC COMMON /SDEV2/ ITRACE, IDEV COMMON /PROB3/IDERIV C CPU TIMER COMMENTED OUT FOR PORTABILITY. C CALL TIMER ( CTIME, 1) M = 0 T = 0.0D0 IDEV = 4 ITRACE = 1 WRITE(IDEV,9)NPOLY, NEL 9 FORMAT(' TEST PROBLEM 3'/' ***********'/' POLY OF DEGREE =',I4, 1 ' NO OF ELEMENTS = ',I4) DO 10 I = 1,IBK 10 XBK(I) = -1.0D0 + 2.0D0 * (I-1.0D0)/(IBK -1.0D0) C INITIALISE THE P.D.E. WORKSPACE ITIME = 1 CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) IF(ITIME .EQ. -1)THEN WRITE(IDEV, 15) 15 FORMAT(' INITCC ROUTINE RETURNED ITIME = -1 - RUN HALTED ') GOTO 100 END IF C SETUP DASSL PARAMETERS RTOL = 1.0D-5 ATOL = 1.0D-5 DO 20 I = 1,11 20 INFO(I) = 0 C INFO(11)= 1 C BANDED MATRIX OPTION WHEN INFO(6) = 1 IF(INFO(6) .EQ. 1)THEN IWORK(1) = IBAND IWORK(2) = IBAND END IF T = 0.0D0 30 TOUT = T + 0.1D0 CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL, 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC) IF( IDID .LT. 0 )THEN C DASSL FAILED TO FINISH INTEGRATION. WRITE(IDEV,40)T,IDID 40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3) GOTO 100 ELSE C DASSL INTEGRATED TO T = TOUT C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION. IDERIV = 0 CALL ERROR( Y, NPDE, NPTS, X, M, ENORM, GERR, T, RTOL, ATOL, 1 ITRACE, WKRES, NWKRES) IFL = 0 ITYPE = 2 DO 45 I = 1,NPTS 45 Z(I) = X(I) CALL INTERC(Z,DYCALC,NPTS,Y,NEQ,NPDE,IFL,ITYPE,WKRES,NWKRES) IDERIV = 1 CALL EXACT(T, NPDE, NPTS, X, DYDX) DO 50 I = 1,NPTS GERRDX = ABS( DYDX(I) - DYCALC(1,I,2)) WRITE(IDEV,49)X(I),DYDX(I),DYCALC(1,I,2),GERRDX 49 FORMAT(' X =',D11.3,' TRUE = ',D11.3,' CALC= ',D11.3,' ERR=', 1 D11.3) 50 CONTINUE IF(TOUT .LT. 0.99D0 ) GOTO 30 END IF 100 CONTINUE C CALL TIMER(CTIME, 2) WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13), CTIME 110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4,' CPU=',D11.3) STOP END C EXAMPLE PROBLEM THREE C ********************* C THIS PROBLEM IS DEFINED BY C -1 C U = ( C U ) - C * EXP(-2U) + EXP(-U) , X IN (-1,0) C T 1 X X 1 C AND C -1 C U = ( C U ) - C * EXP(-2U) + EXP(-U) , X IN (0,1) C T 2 X X 2 C WHERE C C = 0.1 AND C = 1.0 C 1 2 C C THE LEFT BOUNDARY CONDITION AT X =-1 (LEFT = .TRUE. ) IS GIVEN BY C U(-1,T) = LOG ( - C + T + P) C 1 C THE RIGHT BOUNDARY CONDITION IS (LEFT = .FALSE.) C U( 1,T) + (C + T + P ) U = LOG ( - C + T + P) + 1.0D0 C X C C THE INITIAL CONDITION IS GIVEN BY THE EXACT SOLUTION ; C U( X, T ) = LOG ( C X + T + P ) , X IN ( -1, 0) C 1 C U( X, T ) = LOG ( C X + T + P ) , X IN ( 0, 1) C 2 C********************************************************************** SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V) C ROUTINE FOR P.D.E. INITIAL VALUES. INTEGER NPDE, NPTS, NV DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), V(1), T T = 0.0D0 C V(1) IS A DUMMY VARIABLE AS THERE ARE NO COUPLED O.D.E.S CALL EXACT( T, NPDE, NPTS, X, U ) RETURN END C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R, 1 NV, V, VDOT, IRES) C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E. C THE P.D.E. IS WRITEN AS -M M C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U )) C X T TX X T TX X C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE. C DEFINITIONS FOR THE MODEL PROBLEM ARE GIVEN C NOTE NV = 0 : THERE IS NO O.D.E PART. INTEGER NPDE, NPTL, I, J, NV, IRES DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), 1 UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V, VDOT, 2 UTDX(NPDE,NPTL), C IF(X(1) .LT. 0.0D0 .AND. X(NPTL) .LE. 0.0D0)THEN C ELEMENT TO LEFT OF THE INTERFACE AT 0.0 C = 0.1D0 ELSE C = 1.0D0 END IF DO 10 I = 1,NPTL R(1,I) = DUDX(1,I) /C Q(1,I) = UDOT(1,I) - DEXP(-U(1,I))- DEXP(-2.0D0*U(1,I))* C 10 CONTINUE RETURN END C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT, 1 NV, V, VDOT, IRES) C BOUNDARY CONDITIONS ROUTINE INTEGER NPDE, NV, IRES DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), C2, 1 UX(NPDE), V, VDOT, UDOT(NPDE), UTDX(NPDE) LOGICAL LEFT IF(LEFT) THEN BETA (1) = 0.0D0 GAMMA(1) = U(1) - DLOG( -0.1 + T + 1.0D0) ELSE C2 = 1.0D0 BETA (1) = C2 * ( C2 + T + 1.0D0) GAMMA(1) = U(1) - DLOG( C2 + T + 1.0D0) + 1.0D0 END IF RETURN END C C DUMMY O.D.E. ROUTINE AS NV IS ZERO SUBROUTINE SODEFN RETURN END C EXACT SOLUTION SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U) C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN) INTEGER NPDE, NPTS, I, IDERIV DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME, C COMMON /PROB3/ IDERIV IF(IDERIV .EQ. 0)THEN DO 10 I = 1,NPTS C = 1.0D0 IF(X(I) .LT. 0.0D0)C = 0.1D0 10 U(1,I) = DLOG( TIME + 1.0D0 + C * X(I)) ELSE DO 20 I = 1,NPTS C = 1.0D0 IF(X(I) .LT. 0.0D0)C = 0.1D0 U(1,I) = C / ( TIME + 1.0D0 + C * X(I)) IF(X(I) .EQ. 0.0D0) U(1,I) = 0.55D0 / ( TIME + 1.0D0 ) 20 CONTINUE END IF RETURN END c problem 4 C *********************************************************** C BP PROBLEM - VAPOUR EVAPORATION OVER POOL C REGION OF INTEGRATION CONSISTS OF 2 AREAS A VISCOUS SUB-LAYER AND C A TURBULENT REGION, (THE DIVISION OCCURS AT X=0.508D-03). C THE PDE IS DIFFERENT IN EACH REGION. C *********************************************************** C C0 COLLOCATION PARAMETERS PARAMETER ( IBK = 8, NEL = IBK-1 , NPDE = 1, NV = 3, 1 NPOLY = 03, NPTS = NEL*NPOLY+1, NXI = NPTS, 2 NEQ = NPTS * NPDE + NV, 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) + 4 NPDE * 8 + 3 + NV + NXI, C DDASSL TIME INTEGRATION PARAMETERS 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2, 6 LIW = 20 + NEQ ) C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, 1 IDEV, ITRACE, GRNPTS, IFL, NOUT, KTIME, ITYPE DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ), TINC(11), 1 WKRES(NWKRES), RWORK(LRW), XI(NXI), T, TOUT, RTOL, ATOL, 2 U0,VM,DTX1,DTX2,DM1,DM2,K,SCM,PE,MW,RHO,RT,Q3 3 ,TEND,W Q1,Q2,TEMP, XOUT(100), UOUT(100,1), CPU, XBAR REAL GRX(800), GRY(800), GRZ(800) EXTERNAL PDECHB, DGEJAC C C COMMON BLOCKS TO PASS ACROSS PROBLEM DEPENDENT CONSTANTS. COMMON /C0/ PE,MW,RHO,RT,W COMMON /PDES/ U0,VM,DTX1,DTX2,DM1,DM2,K COMMON /SDEV2/ ITRACE, IDEV C IBM CALL TO SWITCH OFF UNDERFLOW COMMENTED OUT C CALL ERRSET(208, 256, -1, -1, 0) C CPU TIMER COMMENTED OUT FOR PORTABILITY C CALL TIMER (CPU, 1) PE = 0.39005D+4 MW = 0.92142D+2 RHO = 0.3767D+1 RT = 0.8317D+4*0.29815D+3 U0 = 0.3164D+0 VM = 0.147D-04 DTX1 = 0.0D+0 SCM = 1.7D+0 DM1 = VM/SCM K = 0.41D+0 DM2 = 0.0D+0 DTX2 = U0*K W = 0.25D0 GRNPTS = 1 WRITE(IDEV,9)NPOLY, NEL 9 FORMAT(' TEST PROBLEM 4'/' ***********'/' POLY OF DEGREE =',I4, 1 ' NO OF ELEMENTS = ',I4) RTOL = 0.1D-4 ATOL = 0.1D-4 ITRACE = 0 IDEV = 4 WRITE(IDEV,104)RTOL, ATOL, ITRACE, IDEV 104 FORMAT(//' RTOL=',D12.3,' ATOL=',D12.3,' ITRACE AND IDEV=',2I4) C WRITE(4,55)ATOL, RTOL, NPTS 55 FORMAT(//' SOLUTION TO B.P. POOL EVAPORATION PROBLEM USING 1 DASSL INTEGRATOR WITH FULL MATRIX ROUTINES '/ 2 ' ATOL = ',D11.3,' RTOL = ',D11.3,' NPTS = ',I5/) NOUT = 20 XOUT(1) = 0.0D0 XOUT(2) = 0.127D-3 XOUT(3) = 0.254D-3 XOUT(4) = 0.381D-3 XOUT(5) = 0.508D-3 XOUT(6) = 0.635D-3 XOUT(7) = 0.762D-3 XOUT(8) = 0.889D-3 XOUT(9) = 0.1D-2 XOUT(10)= 0.3D-2 XOUT(11)= 0.5D-2 XOUT(12)= 0.75D-2 XOUT(13)= 0.1D-1 XOUT(14)= 0.3D-1 XOUT(15)= 0.5D-1 XOUT(16)= 0.75D-1 XOUT(17)= 0.1D0 XOUT(18)= 0.15D0 XOUT(19)= 0.2D0 XOUT(20)= 0.22D0 XBAR = XOUT(5) DO 1000 I = 1,NOUT TEMP = DLOG10( 1.0D0 + XOUT(I)/XBAR *2.0D0) WRITE(IDEV,999)I,XOUT(I),TEMP 999 FORMAT(' I=',I3,' XOUT=',D13.5,' LOG10=',D13.5) 1000 CONTINUE C C TEMPORARY VALUES OF XI FOR FIRST CALL TO INICHB DO 291 I = 1,NPTS 291 XI(I) =(I-1.0D0) /(NPTS-1.0D0) C XBK(1) = 0.0D0 XBK(2) = XBAR* 0.5D0 XBK(3) = XBAR XBK(4) = XBAR * 1.5D0 XBK(5) = XBAR * 2.0D0 XBK(6) = XBAR*11.0 XBK(7) = XBAR * 121 XBK(8) = 1.0D0 ITIME = 1 C INITIALISE THE P.D.E. WORKSPACE CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) DO 292 I = 1,NPTS C FINAL VALUES OF XI 292 XI(I) = X(I) CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) IF(ITIME .EQ. -1)THEN WRITE(IDEV, 15) 15 FORMAT(' INITCC ROUTINE RETURNED ITIME = -1 - RUN HALTED ') GOTO 900 ELSE WRITE(IDEV,16)(Y(I), I = 1,NPTS) 16 FORMAT(' INITIAL VALUES ARE =',5D11.3) END IF C SETUP DASSL PARAMETERS DO 20 I = 1,11 20 INFO(I) = 0 C INFO(11)= 1 C BANDED MATRIX OPTION WHEN INFO(6) = 1 IF(INFO(6) .EQ. 1)THEN IWORK(1) = IBAND IWORK(2) = IBAND END IF T = 0.0D0 TINC(1) = 0.0001D0 TINC(2) = 0.0010 TINC(3) = 0.01D0 TINC(4) = 0.050D0 TINC(5) = 0.1D0 TINC(6) = 0.15D0 TINC(7) = 0.25D0 TINC(8) = 0.50D0 TINC(9) = 0.65D0 TINC(10)= 0.80D0 TINC(11)= 1.00D0 TEND = 1.0D0 KTIME = 1 ITYPE = 1 CALL INTERC(XOUT,UOUT,NOUT,Y,NEQ,NPDE,IFL,ITYPE,WKRES,NWKRES) WRITE (IDEV,82) T, (UOUT(I,1),I=1,NOUT,3) GRNPTS = 0 DO 800 I = 1,NOUT GRNPTS = GRNPTS + 1 GRX(GRNPTS) = T GRZ(GRNPTS) = UOUT(I,1)/UOUT(1,1) GRY(GRNPTS) = DLOG10( 1.D0+XOUT(I)/XBAR * 2.0D0) IF(ITRACE .GE.0)WRITE(IDEV,899)GRY(GRNPTS),GRZ(GRNPTS) 800 CONTINUE WRITE(IDEV,81)(XOUT(I), I = 1,NOUT,2) 81 FORMAT (/' T/X', 4X,9D11.3) C TIME LOOP: 100 TOUT = TINC(KTIME) CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL, 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, NWKRES, DGEJAC) C DASSL FAILED TO FINISH INTEGRATION. WRITE(IDEV,40)T,IDID 40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3) IF( IDID .LT. 0 )GOTO 900 C DASSL INTEGRATED TO T = TOUT C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION. CALL INTERC(XOUT,UOUT,NOUT,Y,NEQ,NPDE,IFL,ITYPE,WKRES,NWKRES) 82 FORMAT(1X,F3.1,' U ',9D11.3/) WRITE (IDEV,82) TOUT, (UOUT(I,1),I=1,NOUT,3) WRITE (6,82) TOUT, (UOUT(I,1),I=1,NOUT,3) C C COMPARE RATE OF EVAPORATION Q1 AT SURFACE OF POOL WITH QUANTITY OF C VAPOUR Q2 WHICH PASSES ABOVE END OF POOL Q1 = Y(NEQ-2) Q2 = Y(NEQ-1) Q3 = Y(NEQ) WRITE(IDEV,83) Q1,Q2,Q3 83 FORMAT(' Q1 , Q2 AND Q3 ARE ',3D13.5) C C PUT INTERPOLATED RESULTS IN ARRAY. C I =(KTIME/2) * 2 IF (I .EQ. KTIME)GOTO 91 DO 90 I = 1,NOUT GRNPTS = GRNPTS + 1 GRX(GRNPTS) = TOUT GRZ(GRNPTS) = UOUT(I,1)/UOUT(1,1) GRY(GRNPTS) = DLOG10( 1.D0+XOUT(I)/XBAR * 2.0D0) IF(ITRACE .GE.0)WRITE(IDEV,899)GRY(GRNPTS),GRZ(GRNPTS) 899 FORMAT(' X AND Y VALUES ARE ',2E12.4) 90 CONTINUE 91 KTIME = KTIME + 1 C C CHECK IF INTEGRATION WAS SUCCESSFUL AND WHETHER FURTHER TIME C STEPS NEEDED IF(TOUT.LT.TEND.AND.(IDID.EQ.2 .OR. IDID .EQ. 3)) GO TO 100 WRITE(IDEV,2112)Q1,Q2,DABS(Q3) 2112 FORMAT(' RATE OF EVAPORATION AT SURFACE OF POOL Q1 = ',D14.7,/ - ' QUANTITY OF VAPOUR ABOVE END OF POOL Q2 = ',D14.7,/ - ' ABSOLUTE DIFFERENCE Q3 = ',D11.4,/ - '********************************************************',/) 80 CONTINUE C CALL TIMER(CPU,2) 900 WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13), CPU 110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4,' CPU=',D11.3) STOP END C EXAMPLE PROBLEM FOUR C ********************* C THIS PROBLEM IS DEFINED BY C C C X U = ( C U ) , X IN (0 , XBAR) C 1 T 2 X X C AND C C (C LOG(X) +C ) U = ( C X U ) , X IN (XBAR , 1) C 3 4 T 5 X X C C WHERE -6 C C = 6810.0 C = 8.65 10 C =0.7717 C = 9.313 C = 0.1297 C 1 2 3 4 5 C C THE LEFT BOUNDARY CONDITION AT X =-1 (LEFT = .TRUE. ) IS GIVEN BY C U(0,T) = 0.038475 C C THE RIGHT BOUNDARY CONDITION IS (LEFT = .FALSE.) C U (1,T) = 0 C X C C THE INITIAL CONDITION IS GIVEN BY C U(X,0) = 0 C C THE ALGEBRAIC VARIABLES Q (T) , Q (T) AND Q (T) ARE DEFINED BY C 1 2 3 C C . -7 C Q = -7.983 10 U (0 , T) C 1 X C C -2 1 C Q = 9.4175 10 I P(X) U(X,T) DX C 2 0 C C WHERE P(X) = C X FOR X IN (0, XBAR) C 1 C C P(X) = C LOG(X) + C FOR X IN (XBAR, 1) C 3 4 C AND THE VALUES OF THE CONSTANTS ARE GIVEN ABOVE. C C Q (T) = Q (T) - Q (T) C 3 2 1 C C********************************************************************** SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V) C ROUTINE FOR P.D.E. INITIAL VALUES. INTEGER NPDE, NPTS, NV, I DOUBLE PRECISION X(NPTS), U(NPDE,NPTS),PE,MW,RHO,RT,W,V(3) COMMON/C0/PE,MW,RHO,RT,W DO 10 I= 2,NPTS 10 U(1,I) = 0.0D+0 U(1,1) = (PE*MW)/(RHO*RT) V(1) = 0.0D0 V(2) = 0.0D0 V(3) = 0.0D0 RETURN END C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R, 1 NV, V, VDOT, IRES) C********************************************************************** C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E. C THE P.D.E. IS WRITEN AS -M M C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U )) C X T TX X T TX X C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE. C********************************************************************** INTEGER NPDE, NPTL, I, NV, IRES DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), 1 DM2, UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V(3), 2 K, UTDX(NPDE,NPTL), U0, VM, DTX1, DTX2, DM1, VDOT(3) COMMON /PDES/ U0,VM,DTX1,DTX2,DM1,DM2,K DO 100 I = 1,NPTL IF(X(1) .LT. 0.506D-3 .AND. X(NPTL) .LT. 0.600D-3)THEN C ELEMENT TO LEFT OF THE INTERFACE AT 0.508D-3 Q(1,I) = (X(I)*U0**2)/VM * UDOT(1,I) R(1,I) = (DTX1 + DM1)*DUDX(1,I) ELSE Q(1,I) = ((U0/K)*DLOG(U0*X(I)/VM) + 5.1*U0) * UDOT(1,I) R(1,I) = ((DTX2*X(I)) + DM2)*DUDX(1,I) ENDIF 100 CONTINUE RETURN END C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT, 1 NV, V, VDOT, IRES) C BOUNDARY CONDITIONS ROUTINE INTEGER NPDE, NV, IRES DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), PE,MW,RHO, 1 UX(NPDE), V(3), VDOT(3), UDOT(NPDE), UTDX(NPDE), RT, W LOGICAL LEFT COMMON/C0/PE,MW,RHO,RT,W IF(LEFT) THEN GAMMA(1) = U(1)- (PE*MW)/(RHO*RT) BETA(1) = 0.0D+0 ELSE GAMMA(1) = 0.0D0 BETA(1) = 1.0D0 END IF RETURN END SUBROUTINE SODEFN(T, NV, V, VDOT, NPDE, NXI, X, Y, UXI, RI, 1 UTI, UTXI, VRES, IRES) C ROUTINE FOR AUXILIARY O.D.E.S (IF ANY) IN MASTER EQN. FORM (4.3) INTEGER NPDE, NXI, NV, IRES, NPTL, L, J, I DOUBLE PRECISION T, X(NXI), Y(NXI), UXI(NPDE,NXI), 1 RI(NPDE,NXI), UTI(NPDE,NXI), UTXI(NPDE,NXI), VRES(NV), 2 V(3), VDOT(3), PE,MW,RHO,RT,W,U0,VM,DTX1,DTX2,DM2,K 3 ,DM1, Q2, H, CCRULE COMMON /C0/ PE,MW,RHO,RT,W COMMON /PDES/ U0,VM,DTX1,DTX2,DM1,DM2,K COMMON /SCHSZ5/ NPTL COMMON /SCHSZ6/ CCRULE(50) C VRES(1) = VDOT(1) + W*RHO*DM1*UXI(1,1) Q2 = 0.0D0 DO 3 I = 1,2 J = (NPTL-1) * (I-1) + 1 L = (NPTL-1) * I + 1 H = ( X(L) - X(J)) * 0.5D0 DO 3 II = 1,NPTL IK = J + II - 1 C CLENSHAW - CURTIS QUADRATURE UP TO INTERFACE POINT. Q2 = Q2 + (W*RHO*U0**2)/VM * X(IK) * Y(IK) * CCRULE(II) * H 3 CONTINUE C C CLENSHAW - CURTIS QUADRATURE BEYOND THE INTERFACE POINT. DO 5 I = 3,7 J = (NPTL-1) * (I-1) + 1 L = (NPTL-1) * I + 1 H = ( X(L) - X(J)) * 0.5D0 DO 5 II = 1, NPTL IK = J + II - 1 Q2=Q2 + H* ((U0/K)*DLOG(U0*X(IK)/VM)+5.1*U0) * Y(IK)*CCRULE(II) 1 * W * RHO 5 CONTINUE VRES(2) = V(2) - Q2 VRES(3) = V(3) - (V(2)-V(1)) RETURN END c problem 5 C *********************************************************** C FOURTH ORDER P.D.E. PROBLEM WRITTEN AS ELLIPTIC-PARABOLIC SYSTEM. C C U = K U + UU - U U C XXT XXXX XXX X XX C C *********************************************************** C C C0 COLLOCATION PARAMETERS PARAMETER ( IBK = 21, NEL = IBK-1 , NPDE = 2, NV = 0, 1 NPOLY = 02, NPTS = NEL*NPOLY+1, NXI = 0, 2 NEQ = NPTS * NPDE + NV, 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) + 4 NPDE * 8 + 3 + NV + NXI, C DDASSL TIME INTEGRATION PARAMETERS 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2, 6 LIW = 20 + NEQ ) C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, 1 IDEV, ITRACE, GRNPTS, IFL, NOUT, KTIME, ITYPE, NP DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ), TINC(15), 1 WKRES(NWKRES), RWORK(LRW), XI, T, TOUT, RTOL, ATOL, 3 TEND, K, XOUT(6), UOUT(2,6) EXTERNAL PDECHB, DGEJAC C C COMMON BLOCKS TO PASS ACROSS PROBLEM DEPENDENT CONSTANTS. COMMON /PDES/ K COMMON /SDEV2/ ITRACE, IDEV DATA XOUT(1)/-1.0D+0/, XOUT(2)/-0.6D+0/, XOUT(3)/-0.2D+0/, * XOUT(4)/0.2D+0/, XOUT(5)/0.6D+0/, XOUT(6)/1.0D+0/ WRITE(IDEV,9)NPOLY, NEL 9 FORMAT(' TEST PROBLEM 4'/' ***********'/' POLY OF DEGREE =',I4, 1 ' NO OF ELEMENTS = ',I4) RTOL = 0.1D-4 ATOL = 0.1D-4 ITRACE = 0 IDEV = 4 WRITE(IDEV,104)RTOL, ATOL, ITRACE, IDEV 104 FORMAT(//' RTOL=',D12.3,' ATOL=',D12.3,' ITRACE AND IDEV=',2I4) C WRITE(4,55)ATOL, RTOL, NPTS 55 FORMAT(//' SOLUTION TO FOURTH ORDER P.D.E. PROBLEM USING 1 DASSL INTEGRATOR WITH BANDED MATRIX ROUTINES '/ 2 ' ATOL = ',D11.3,' RTOL = ',D11.3,' NPTS = ',I5/) C C EQUALLY SPACED BREAKPOINTS. C DO 105 I = 1,IBK XBK(I) = -1.0D0 + (I -1.0D0)* 2.D0 / (IBK-1.D0) 105 CONTINUE K = 1.00D0 ITIME = 1 T = 0.0D0 C INITIALISE THE P.D.E. WORKSPACE CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND, 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV) IF(ITIME .EQ. -1)THEN WRITE(IDEV, 15) 15 FORMAT(' INICHB ROUTINE RETURNED ITIME = -1 - RUN HALTED ') GOTO 900 END IF C SETUP DASSL PARAMETERS DO 20 I = 1,11 20 INFO(I) = 0 INFO(6)= 1 INFO(9)= 1 INFO(7)= 1 IWORK(3)= 4 C BANDED MATRIX OPTION WHEN INFO(6) = 1 IF(INFO(6) .EQ. 1)THEN IWORK(1) = IBAND IWORK(2) = IBAND END IF T = 0.0D0 TINC(1) = 0.0001D0 RWORK(2)= TINC(1) * 0.1D0 TINC(2) = 0.0010 TINC(3) = 0.01D0 TINC(4) = 0.1D0 TINC(5) = 1.0D0 TINC(6) = 1.00D1 TINC(7) = 2.00D1 TINC(8) = 4.00D1 TINC(9 )= 6.00D1 TINC(10)= 8.00D1 TINC(11)= 1.00D2 TINC(12)= 1.00D3 TEND = 1.0D3 KTIME = 1 WRITE(IDEV,83)(XOUT(I),I = 1,6) C TIME LOOP: 100 TOUT = TINC(KTIME) IF(KTIME.GT.1)RWORK(2) = 0.05D0 *(TOUT- TINC(KTIME-1)) IF(KTIME .EQ.12)THEN INFO(4) = 1 RWORK(1) = TEND END IF C CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL, 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, NWKRES, DGEJAC) C DASSL FAILED TO FINISH INTEGRATION. WRITE(IDEV,40)T,IDID,Y(1),Y(2),Y(NEQ-1), Y(NEQ) 40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3/ 1 ' LEFT SOL=',2D11.3,' RIGHT SOL=',2D11.3) IF( IDID .LT. 0 )GOTO 900 C DASSL INTEGRATED TO T = TOUT C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION. ITYPE = 1 NP = 6 CALL INTERC( XOUT, UOUT, NP, Y, NEQ, NPDE, IFLAG, 1 ITYPE, WKRES, NWKRES) WRITE(IDEV,82)(UOUT(1,I),I = 1,6) WRITE(IDEV,84)(UOUT(2,I),I = 1,6) 83 FORMAT(1X,'X',9D11.3/) 82 FORMAT(1X,'U',9D11.3/) 84 FORMAT(1X,'V',9D11.3/) C 91 KTIME = KTIME + 1 C C CHECK IF INTEGRATION WAS SUCCESSFUL AND WHETHER FURTHER TIME C STEPS NEEDED IF(TOUT.LT.TEND.AND.(IDID.EQ.2 .OR. IDID .EQ. 3)) GO TO 100 80 CONTINUE 900 WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13) 110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4) STOP END C EXAMPLE PROBLEM FIVE C ********************* C THIS PROBLEM IS DEFINED BY C C V = U C XX C AND C C V = ( K V ) + U V - U V , X IN (-1 , 1) C T X X X X C WHERE C K = 0.15 C C THE LEFT BOUNDARY CONDITION AT X =-1 (LEFT = .TRUE. ) ARE GIVEN BY C C U = 1 U = 0.0 C X C THE RIGHT BOUNDARY CONDITION ARE (LEFT = .FALSE.) C C U = -1 U = 0.0D0 C X C THE INITIAL CONDITION IS GIVEN BY C C U(X,0) = -SIN ( PI /2 X ) C********************************************************************** SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V) C ROUTINE FOR P.D.E. INITIAL VALUES. PARAMETER (PIBY2 = 1.5707963D0) INTEGER NPDE, NPTS, NV, I DOUBLE PRECISION X(NPTS), U(NPDE,NPTS),V DO 10 I= 1,NPTS U(1,I) = -SIN( PIBY2 * X(I) ) 10 U(2,I) = - PIBY2**2 * U(1,I) RETURN END C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R, 1 NV, V, VDOT, IRES) C********************************************************************** C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E. C THE P.D.E. IS WRITEN AS -M M C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U )) C X T TX X T TX X C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE. C********************************************************************** INTEGER NPDE, NPTL, I, NV, IRES DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), 1 UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V, 2 K, UTDX(NPDE,NPTL), VDOT COMMON /PDES/ K DO 100 I = 1,NPTL Q(1,I) = U(2,I) R(1,I) = DUDX(1,I) Q(2,I) = UDOT(2,I) - U(1,I)*DUDX(2,I) + DUDX(1,I)*U(2,I) R(2,I) = K*DUDX(2,I) 100 CONTINUE RETURN END C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT, 1 NV, V, VDOT, IRES) C BOUNDARY CONDITIONS ROUTINE INTEGER NPDE, NV, IRES DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), 1 UX(NPDE), V, VDOT, UDOT(NPDE), UTDX(NPDE) LOGICAL LEFT IF(LEFT) THEN GAMMA(1) = 0.0D0 BETA(1) = 1.0D+0 GAMMA(2) = U(1) - 1.0D0 BETA(2) = 0.0D+0 ELSE GAMMA(1) = 0.0D0 BETA(1) = 1.0D+0 GAMMA(2) = U(1) + 1.0D0 BETA(2) = 0.0D+0 END IF RETURN END SUBROUTINE SODEFN(T, NV, V, VDOT, NPDE, NXI, X, Y, UXI, RI, 1 UTI, UTXI, VRES, IRES) C ROUTINE FOR AUXILIARY O.D.E.S (IF ANY) IN MASTER EQN. FORM (4.3) INTEGER NPDE, NXI, NV, IRES, NPTL, L, J, I DOUBLE PRECISION T, X(NXI), Y(NXI), UXI(NPDE,NXI), 1 RI(NPDE,NXI), UTI(NPDE,NXI), UTXI(NPDE,NXI), VRES(NV), 2 V, VDOT C DUMMY ROUTINE AS THERE ARE NO O.D.E.S RETURN END c c main body of the software follows................ c SUBROUTINE DDASSL (RES,NEQ,T,Y,YPRIME,TOUT, * INFO,RTOL,ATOL,IDID, * RWORK,LRW,IWORK,LIW,RPAR,IPAR, * JAC) C C***BEGIN PROLOGUE DDASSL C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***CATEGORY NO. D2A2 C***KEYWORDS DIFFERENTIAL/ALGEBRAIC,BACKWARD DIFFERENTIATION FORMULAS C IMPLICIT DIFFERENTIAL SYSTEMS C***AUTHOR PETZOLD,LINDA R. C APPLIED MATHEMATICS DIVISION 8331 C SANDIA NATIONAL LABORATORIES C LIVERMORE, CA. 94550 C***PURPOSE DIFFERENTIAL/ALGEBRAIC SYSTEM SOLVER C***DESCRIPTION C --------------------------------------------------------------------- C C this code solves a system of differential/ C algebraic equations of the form C g(t,y,yprime) = 0. C C subroutine ddassl uses the backward C differentiation formulas of orders one C through five to solve a system of the above C form for y and yprime. values for y C and yprime at the initial time must C be given as input. these values must C be consistent, (that is. if t,y,yprime C are the given initial values, they must C satisfy g(t,y,yprime) = 0.) C the subroutine solves the system from t to tout. it is C easy to continue the solution to get results C at additional tout. this is the interval C mode of operation. intermediate results can C also be obtained easily by using the intermediate- C output capability. C C C ------------description of arguments to ddassl----------------------- C ------------(an overview)-------------------------------------------- C C the parameters are C C res -- this is a subroutine which you provide C to define the differential/algebraic C system C C neq -- this is the number of equations C to be solved C C t -- this is the current value of the C independent variable. C C tout -- this is a point at which a solution C is desired. C C info(*) -- the basic task of the code is C to solve the system from t to C tout and return an answer at tout. C info(*) is an integer array which is C used to communicate exactly how you C want this task to be carried out. C C y(*) -- this array contains the solution C components at t C C yprime(*) -- this array contains the derivatives C of the solution components at t C C rtol,atol -- these quantities represent C absolute and relative error C tolerances which you provide to indicate C how accurately you wish the solution C to be computed. you may choose them C to be both scalars or else both C vectors. C C idid -- this scalar quantity is an indicator reporting C what the code did. you must monitor this C integer variable to decide what action to C take next. C C rwork(*),lrw -- rwork(*) is a real work array of C length lrw which provides the code C with needed storage space. C C iwork(*),liw -- iwork(*) is an integer work array C of length liw which provides the code C with needed storage space. C C rpar,ipar -- these are real and integer parameter C arrays which you can use for C communication between your calling C program and the res subroutine C (and the jac subroutine) C C jac -- this is the name of a subroutine which you C may choose to provide for defining C a matrix of partial derivatives C described below. C C quantities which are used as input items are C neq,t,y(*),yprime(*),tout,info(*), C rtol,atol,rwork(1),rwork(2),rwork(3),lrw,iwork(1), C iwork(2),iwork(3),and liw. C C quantities which may be altered by the code are C t,y(*),yprime(*),info(1),rtol,atol, C idid,rwork(*) and iwork(*) C C ----------input-what to do on the first call to ddassl--------------- C C C the first call of the code is defined to be the start of each new C problem. read through the descriptions of all the following items, C provide sufficient storage space for designated arrays, set C appropriate variables for the initialization of the problem, and C give information about how you want the problem to be solved. C C C res -- provide a subroutine of the form C subroutine res(t,y,yprime,delta,ires,rpar,ipar) C to define the system of differential/algebraic C equations which is to be solved. for the given values C of t,y and yprime, the subroutine should C return the residual of the differential/algebraic C system C delta = g(t,y,yprime) C (delta(*) is a vector of length neq which is C output for res.) C C subroutine res must not alter t,y or yprime. C you must declare the name res in an external C statement in your program that calls ddassl. C you must dimension y,yprime and delta in res. C C ires is an integer flag which is always equal to C zero on input. subroutine res should alter ires C only if it encounters an illegal value of y or C a stop condition. set ires = -1 if an input value C is illegal, and ddassl will try to solve the problem C without getting ires = -1. if ires = -2, ddassl C will return control to the calling program C with idid = -11. C C rpar and ipar are real and integer parameter arrays which C you can use for communication between your calling program C and subroutine res. they are not altered by ddassl. if you C do not need rpar or ipar, ignore these parameters by treat- C ing them as dummy arguments. if you do choose to use them, C dimension them in your calling program and in res as arrays C of appropriate length. C C neq -- set it to the number of differential equations. C (neq .ge. 1) C C t -- set it to the initial point of the integration. C t must be defined as a variable. C C y(*) -- set this vector to the initial values of the neq solution C components at the initial point. you must dimension y of C length at least neq in your calling program. C C yprime(*) -- set this vector to the initial values of C the neq first derivatives of the solution C components at the initial point. you C must dimension yprime at least neq C in your calling program. if you do not C know initial values of some of the solution C components, see the explanation of info(11). C C tout - set it to the first point at which a solution C is desired. you can not take tout = t. C integration either forward in t (tout .gt. t) or C backward in t (tout .lt. t) is permitted. C C the code advances the solution from t to tout using C step sizes which are automatically selected so as to C achieve the desired accuracy. if you wish, the code will C return with the solution and its derivative at C intermediate steps (intermediate-output mode) so that C you can monitor them, but you still must provide tout in C accord with the basic aim of the code. C C the first step taken by the code is a critical one C because it must reflect how fast the solution changes near C the initial point. the code automatically selects an C initial step size which is practically always suitable for C the problem. by using the fact that the code will not step C past tout in the first step, you could, if necessary, C restrict the length of the initial step size. C C for some problems it may not be permissable to integrate C past a point tstop because a discontinuity occurs there C or the solution or its derivative is not defined beyond C tstop. when you have declared a tstop point (see info(4) C and rwork(1)), you have told the code not to integrate C past tstop. in this case any tout beyond tstop is invalid C input. C C info(*) - use the info array to give the code more details about C how you want your problem solved. this array should be C dimensioned of length 15, though ddassl uses C only the first nine entries. you must respond to all of C the following items which are arranged as questions. the C simplest use of the code corresponds to answering all C questions as yes ,i.e. setting all entries of info to 0. C C info(1) - this parameter enables the code to initialize C itself. you must set it to indicate the start of every C new problem. C C **** is this the first call for this problem ... C yes - set info(1) = 0 C no - not applicable here. C see below for continuation calls. **** C C info(2) - how much accuracy you want of your solution C is specified by the error tolerances rtol and atol. C the simplest use is to take them both to be scalars. C to obtain more flexibility, they can both be vectors. C the code must be told your choice. C C **** are both error tolerances rtol, atol scalars ... C yes - set info(2) = 0 C and input scalars for both rtol and atol C no - set info(2) = 1 C and input arrays for both rtol and atol **** C C info(3) - the code integrates from t in the direction C of tout by steps. if you wish, it will return the C computed solution and derivative at the next C intermediate step (the intermediate-output mode) or C tout, whichever comes first. this is a good way to C proceed if you want to see the behavior of the solution. C if you must have solutions at a great many specific C tout points, this code will compute them efficiently. C C **** do you want the solution only at C tout (and not at the next intermediate step) ... C yes - set info(3) = 0 C no - set info(3) = 1 **** C C info(4) - to handle solutions at a great many specific C values tout efficiently, this code may integrate past C tout and interpolate to obtain the result at tout. C sometimes it is not possible to integrate beyond some C point tstop because the equation changes there or it is C not defined past tstop. then you must tell the code C not to go past. C C **** can the integration be carried out without any C restrictions on the independent variable t ... C yes - set info(4)=0 C no - set info(4)=1 C and define the stopping point tstop by C setting rwork(1)=tstop **** C C info(5) - to solve differential/algebraic problems it is C necessary to use a matrix of partial derivatives of the C system of differential equations. if you do not C provide a subroutine to evaluate it analytically (see C description of the item jac in the call list), it will C be approximated by numerical differencing in this code. C although it is less trouble for you to have the code C compute partial derivatives by numerical differencing, C the solution will be more reliable if you provide the C derivatives via jac. sometimes numerical differencing C is cheaper than evaluating derivatives in jac and C sometimes it is not - this depends on your problem. C C **** do you want the code to evaluate the partial C derivatives automatically by numerical differences .. C yes - set info(5)=0 C no - set info(5)=1 C and provide subroutine jac for evaluating the C matrix of partial derivatives **** C C info(6) - ddassl will perform much better if the matrix of C partial derivatives, dg/dy + cj*dg/dyprime, C (here cj is a scalar determined by ddassl) C is banded and the code is told this. in this C case, the storage needed will be greatly reduced, C numerical differencing will be performed much cheaper, C and a number of important algorithms will execute much C faster. the differential equation is said to have C half-bandwidths ml (lower) and mu (upper) if equation i C involves only unknowns y(j) with C i-ml .le. j .le. i+mu C for all i=1,2,...,neq. thus, ml and mu are the widths C of the lower and upper parts of the band, respectively, C with the main diagonal being excluded. if you do not C indicate that the equation has a banded matrix of partial C derivatives C the code works with a full matrix of neq**2 elements C (stored in the conventional way). computations with C banded matrices cost less time and storage than with C full matrices if 2*ml+mu .lt. neq. if you tell the C code that the matrix of partial derivatives has a banded C structure and you want to provide subroutine jac to C compute the partial derivatives, then you must be careful C to store the elements of the matrix in the special form C indicated in the description of jac. C C **** do you want to solve the problem using a full C (dense) matrix (and not a special banded C structure) ... C yes - set info(6)=0 C no - set info(6)=1 C and provide the lower (ml) and upper (mu) C bandwidths by setting C iwork(1)=ml C iwork(2)=mu **** C C C info(7) -- you can specify a maximum (absolute value of) C stepsize, so that the code C will avoid passing over very C large regions. C C **** do you want the code to decide C on its own maximum stepsize? C yes - set info(7)=0 C no - set info(7)=1 C and define hmax by setting C rwork(2)=hmax **** C C info(8) -- differential/algebraic problems C may occaisionally suffer from C severe scaling difficulties on the C first step. if you know a great deal C about the scaling of your problem, you can C help to alleviate this problem by C specifying an initial stepsize ho. C C **** do you want the code to define C its own initial stepsize? C yes - set info(8)=0 C no - set info(8)=1 C and define ho by setting C rwork(3)=ho **** C C info(9) -- if storage is a severe problem, C you can save some locations by C restricting the maximum order maxord. C the default value is 5. for each C order decrease below 5, the code C requires neq fewer locations, however C it is likely to be slower. in any C case, you must have 1 .le. maxord .le. 5 C **** do you want the maximum order to C default to 5? C yes - set info(9)=0 C no - set info(9)=1 C and define maxord by setting C iwork(3)=maxord **** C C info(10) --if you know that the solutions to your equations wil C always be nonnegative, it may help to set this C parameter. however, it is probably best to C try the code without using this option first, C and only to use this option if that doesn't C work very well. C **** do you want the code to solve the problem without C invoking any special nonnegativity constraints? C yes - set info(10)=0 C no - set info(10)=1 C C info(11) --ddassl normally requires the initial t, C y, and yprime to be consistent. that is, C you must have g(t,y,yprime) = 0 at the initial C time. if you do not know the initial C derivative precisely, you can let ddassl try C to compute it. C **** are the initial t, y, yprime consistent? C yes - set info(11) = 0 C no - set info(11) = 1, C and set yprime to an initial approximation C to yprime. (if you have no idea what C yprime should be, set it to zero. note C that the initial y should be such C that there must exist a yprime so that C g(t,y,yprime) = 0.) C C rtol, atol -- you must assign relative (rtol) and absolute (atol C error tolerances to tell the code how accurately you wan C the solution to be computed. they must be defined as C variables because the code may change them. you have two C choices -- C both rtol and atol are scalars. (info(2)=0) C both rtol and atol are vectors. (info(2)=1) C in either case all components must be non-negative. C C the tolerances are used by the code in a local error tes C at each step which requires roughly that C abs(local error) .le. rtol*abs(y)+atol C for each vector component. C (more specifically, a root-mean-square norm is used to C measure the size of vectors, and the error test uses the C magnitude of the solution at the beginning of the step.) C C the true (global) error is the difference between the tr C solution of the initial value problem and the computed C approximation. practically all present day codes. C including this one, control the local error at each step C and do not even attempt to control the global error C directly. C usually, but not always, the true accuracy of C the computed y is comparable to the error tolerances. th C code will usually, but not always, deliver a more accura C solution if you reduce the tolerances and integrate agai C by comparing two such solutions you can get a fairly C reliable idea of the true error in the solution at the C bigger tolerances. C C setting atol=0. results in a pure relative error test on C that component. setting rtol=0. results in a pure absolu C error test on that component. a mixed test with non-zero C rtol and atol corresponds roughly to a relative error C test when the solution component is much bigger than ato C and to an absolute error test when the solution componen C is smaller than the threshold atol. C C the code will not attempt to compute a solution at an C accuracy unreasonable for the machine being used. it wil C advise you if you ask for too much accuracy and inform C you as to the maximum accuracy it believes possible. C C rwork(*) -- dimension this real work array of length lrw in your C calling program. C C lrw -- set it to the declared length of the rwork array. C you must have C lrw .ge. 40+(maxord+4)*neq+neq**2 C for the full (dense) jacobian case (when info(6)=0), or C lrw .ge. 40+(maxord+4)*neq+(2*ml+mu+1)*neq C for the banded user-defined jacobian case C (when info(5)=1 and info(6)=1), or C lrw .ge. 40+(maxord+4)*neq+(2*ml+mu+1)*neq C +2*(neq/(ml+mu+1)+1) C for the banded finite-difference-generated jacobian case C (when info(5)=0 and info(6)=1) C C iwork(*) -- dimension this integer work array of length liw in C your calling program. C C liw -- set it to the declared length of the iwork array. C you must have liw .ge. 20+neq C C rpar, ipar -- these are parameter arrays, of real and integer C type, respectively. you can use them for communication C between your program that calls ddassl and the C res subroutine (and the jac subroutine). they are not C altered by ddassl. if you do not need rpar or ipar, igno C these parameters by treating them as dummy arguments. if C you do choose to use them, dimension them in your callin C program and in res (and in jac) as arrays of appropriate C length. C C jac -- if you have set info(5)=0, you can ignore this parameter C by treating it as a dummy argument. otherwise, you must C provide a subroutine of the form C jac(t,y,yprime,pd,cj,rpar,ipar) C to define the matrix of partial derivatives C pd=dg/dy+cj*dg/dyprime C cj is a scalar which is input to jac. C for the given values of t,y,yprime, the C subroutine must evaluate the non-zero partial C derivatives for each equation and each solution C compowent, and store these values in the C matrix pd. the elements of pd are set to zero C before each call to jac so only non-zero elements C need to be defined. C C subroutine jac must not alter t,y,(*),yprime(*),or cj. C you must declare the name jac in an C external statement in your program that calls C ddassl. you must dimension y, yprime and pd C in jac. C C the way you must store the elements into the pd matrix C depends on the structure of the matrix which you C indicated by info(6). C *** info(6)=0 -- full (dense) matrix *** C when you evaluate the (non-zero) partial derivative C of equation i with respect to variable j, you must C store it in pd according to C pd(i,j) = * dg(i)/dy(j)+cj*dg(i)/dyprime(j)* C *** info(6)=1 -- banded jacobian with ml lower and mu C upper diagonal bands (refer to info(6) description o C ml and mu) *** C when you evaluate the (non-zero) partial derivative C of equation i with respect to variable j, you must C store it in pd according to C irow = i - j + ml + mu + 1 C pd(irow,j) = *dg(i)/dy(j)+cj*dg(i)/dyprime(j)* C rpar and ipar are real and integer parameter arrays whic C you can use for communication between your calling C program and your jacobian subroutine jac. they are not C altered by ddassl. if you do not need rpar or ipar, igno C these parameters by treating them as dummy arguments. if C you do choose to use them, dimension them in your callin C program and in jac as arrays of appropriate length. C C C C optionally replaceable norm routine: C ddassl uses a weighted norm ddanrm to measure the size C of vectors such as the estimated error in each step. C a function subprogram C double precision function ddanrm(neq,v,wt,rpar,ipar) C dimension v(neq),wt(neq) C is used to define this norm. here, v is the vector C whose norm is to be computed, and wt is a vector of C weights. a ddanrm routine has been included with ddassl C which computes the weighted root-mean-square norm C given by C ddanrm=sqrt((1/neq)*sum(v(i)/wt(i))**2) C this norm is suitable for most problems. in some C special cases, it may be more convenient and/or C efficient to define your own norm by writing a function C subprogram to be called instead of ddanrm. this should C however, be attempted only after careful thought and C consideration. C C C------output-after any return from ddassl---- C C the principal aim of the code is to return a computed solution at C tout, although it is also possible to obtain intermediate results C along the way. to find out whether the code achieved its goal C or if the integration process was interrupted before the task was C completed, you must check the idid parameter. C C C t -- the solution was successfully advanced to the C output value of t. C C y(*) -- contains the computed solution approximation at t. C C yprime(*) -- contains the computed derivative C approximation at t C C idid -- reports what the code did C C *** task completed *** C reported by positive values of idid C C idid = 1 -- a step was successfully taken in the C intermediate-output mode. the code has not C yet reached tout. C C idid = 2 -- the integration to tout was successfully C completed (t=tout) by stepping exactly to tout. C C idid = 3 -- the integration to tout was successfully C completed (t=tout) by stepping past tout. C y(*) is obtained by interpolation. C yprime(*) is obtained by interpolation. C C *** task interrupted *** C reported by negative values of idid C C idid = -1 -- a large amount of work has been expended. C (about 500 steps) C C idid = -2 -- the error tolerances are too stringent. C C idid = -3 -- the local error test cannot be satisfied C because you specified a zero component in atol C and the corresponding computed solution C component is zero. thus, a pure relative error C test is impossible for this component. C C idid = -6 -- ddassl had repeated error test C failures on the last attempted step. C C idid = -7 -- the corrector could not converge. C C idid = -8 -- the matrix of partial derivatives C is singular. C C idid = -9 -- the corrector could not converge. C there were repeated error test failures C in this step. C C idid =-10 -- the corrector could not converge C because ires was equal to minus one. C C idid =-11 -- ires equal to -2 was encountered C and control is being returned to the C calling program. C C idid =-12 -- ddassl failed to compute the initial C yprime. C C C C idid = -13,..,-32 -- not applicable for this code C C *** task terminated *** C reported by the value of idid=-33 C C idid = -33 -- the code has encountered trouble from which C it cannot recover. a message is printed C explaining the trouble and control is returned C to the calling program. for example, this occurs C when invalid input is detected. C C rtol, atol -- these quantities remain unchanged except when C idid = -2. in this case, the error tolerances have been C increased by the code to values which are estimated to b C appropriate for continuing the integration. however, the C reported solution at t was obtained using the input valu C of rtol and atol. C C rwork, iwork -- contain information which is usually of no C interest to the user but necessary for subsequent calls. C however, you may find use for C C rwork(3)--which contains the step size h to be C attempted on the next step. C C rwork(4)--which contains the current value of the C independent variable, i.e. the farthest point C integration has reached. this will be different C from t only when interpolation has been C performed (idid=3). C C rwork(7)--which contains the stepsize used C on the last successful step. C C iwork(7)--which contains the order of the method to C be attempted on the next step. C C iwork(8)--which contains the order of the method used C on the last step. C C iwork(11)--which contains the number of steps taken so f C C iwork(12)--which contains the number of calls to res C so far. C C iwork(13)--which contains the number of evaluations of C the matrix of partial derivatives needed so far C C iwork(14)--which contains the total number C of error test failures so far. C C iwork(15)--which contains the total number C of convergence test failures so far. C (includes singular iteration matrix C failures.) C C C C input -- what to do to continue the integration C (calls after the first) ** C C this code is organized so that subsequent calls to continue the C integration involve little (if any) additional effort on your C part. you must monitor the idid parameter in order to determine C what to do next. C C recalling that the principal task of the code is to integrate C from t to tout (the interval mode), usually all you will need C to do is specify a new tout upon reaching the current tout. C C do not alter any quantity not specifically permitted below, C in particular do not alter neq,t,y(*),yprime(*),rwork(*),iwork(*) C or the differential equation in subroutine res. any such C alteration constitutes a new problem and must be treated as such, C i.e. you must start afresh. C C you cannot change from vector to scalar error control or vice C versa (info(2)) but you can change the size of the entries of C rtol, atol. increasing a tolerance makes the equation easier C to integrate. decreasing a tolerance will make the equation C harder to integrate and should generally be avoided. C C you can switch from the intermediate-output mode to the C interval mode (info(3)) or vice versa at any time. C C if it has been necessary to prevent the integration from going C past a point tstop (info(4), rwork(1)), keep in mind that the C code will not integrate to any tout beyound the currently C specified tstop. once tstop has been reached you must change C the value of tstop or set info(4)=0. you may change info(4) C or tstop at any time but you must supply the value of tstop in C rwork(1) whenever you set info(4)=1. C C do not change info(5), info(6), iwork(1), or iwork(2) C unless you are going to restart the code. C C *** following a completed task *** C if C idid = 1, call the code again to continue the integration C another step in the direction of tout. C C idid = 2 or 3, define a new tout and call the code again. C tout must be different from t. you cannot change C the direction of integration without restarting. C C *** following an interrupted task *** C to show the code that you realize the task was C interrupted and that you want to continue, you C must take appropriate action and set info(1) = 1 C if C idid = -1, the code has taken about 500 steps. C if you want to continue, set info(1) = 1 and C call the code again. an additional 500 steps C will be allowed. C C C idid = -2, the error tolerances rtol, atol have been C increased to values the code estimates appropriate C for continuing. you may want to change them C yourself. if you are sure you want to continue C with relaxed error tolerances, set info(1)=1 and C call the code again. C C idid = -3, a solution component is zero and you set the C corresponding component of atol to zero. if you C are sure you want to continue, you must first C alter the error criterion to use positive values C for those components of atol corresponding to zero C solution components, then set info(1)=1 and call C the code again. C C idid = -4,-5 --- cannot occur with this code C C idid = -6, repeated error test failures occurred on the C last attempted step in ddassl. a singularity in the C solution may be present. if you are absolutely C certain you want to continue, you should restart C the integration.(provide initial values of y and C yprime which are consistent) C C idid = -7, repeated convergence test failures occurred C on the last attempted step in ddassl. an inaccurate o C illconditioned jacobian may be the problem. if you C are absolutely certain you want to continue, you C should restart the integration. C C idid = -8, the matrix of partial derivatives is singular. C some of your equations may be redundant. C ddassl cannot solve the problem as stated. C it is possible that the redundant equations C could be removed, and then ddassl could C solve the problem. it is also possible C that a solution to your problem either C does not exist or is not unique. C C idid = -9, ddassl had multiple convergence test C failures, preceeded by multiple error C test failures, on the last attempted step. C it is possible that your problem C is ill-posed, and cannot be solved C using this code. or, there may be a C discontinuity or a singularity in the C solution. if you are absolutely certain C you want to continue, you should restart C the integration. C C idid =-10, ddassl had multiple convergence test failures C because ires was equal to minus one. C if you are absolutely certain you want C to continue, you should restart the C integration. C C idid =-11, ires=-2 was encountered, and control is being C returned to the calling program. C C idid =-12, ddassl failed to compute the initial yprime. C this could happen because the initial C approximation to yprime was not very good, or C if a yprime consistent with the initial y C does not exist. the problem could also be caused C by an inaccurate or singular iteration matrix. C C C C idid = -13,..,-32 --- cannot occur with this code C C *** following a terminated task *** C if idid= -33, you cannot continue the solution of this C problem. an attempt to do so will result in your C run being terminated. C C --------------------------------------------------------------------- C C***REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC C SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637, C SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982. C***ROUTINES CALLED DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,XERRWV,D1MACH C***COMMON BLOCKS DDA001 C***END PROLOGUE DDASSL C C IMPLICIT REAL*8 (A-H,O-Z) LOGICAL DONE EXTERNAL RES,JAC DIMENSION Y(1),YPRIME(1) DIMENSION INFO(15) DIMENSION RWORK(1),IWORK(1) DIMENSION RTOL(1),ATOL(1) DIMENSION RPAR(1),IPAR(1) COMMON/SDEV2/ ITRACE, IDEV COMMON /SDDTR/ TERKP1, TERK, TERKM1, TERKM2 COMMON/DDA001/NPD,NTEMP, * LML,LMU,LMXORD,LMTYPE, * LNST,LNRE,LNJE,LETF,LCTF,LIPVT DATA LTSTOP,LHMAX,LH,LTN, * LCJ,LCJOLD,LHOLD,LS,LROUND, * LALPHA,LBETA,LGAMMA, * LPSI,LSIGMA,LDELTA * /1,2,3,4, * 5,6,7,8,9, * 11,17,23, * 29,35,41/ IF(INFO(1).NE.0)GO TO 100 C C----------------------------------------------------------------------- C this block is executed for the initial call only. C it contains checking of inputs and initializations. C----------------------------------------------------------------------- C C first check info array to make sure all elements of info C are either zero or one. DO 10 I=2,11 IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701 10 CONTINUE C IF(NEQ.LE.0)GO TO 702 C C set pointers into iwork LML=1 LMU=2 LMXORD=3 LMTYPE=4 LJCALC=5 LPHASE=6 LK=7 LKOLD=8 LNS=9 LNSTL=10 LNST=11 LNRE=12 LNJE=13 LETF=14 LCTF=15 LIPVT=21 LIWM=1 C C check and compute maximum order MXORD=5 IF(INFO(9).EQ.0)GO TO 20 MXORD=IWORK(LMXORD) IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703 20 IWORK(LMXORD)=MXORD C C compute mtype,lenpd,lenrw.check ml and mu. IF(INFO(6).NE.0)GO TO 40 LENPD=NEQ**2 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD IF(INFO(5).NE.0)GO TO 30 IWORK(LMTYPE)=2 GO TO 60 30 IWORK(LMTYPE)=1 GO TO 60 40 IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717 IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718 LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ IF(INFO(5).NE.0)GO TO 50 IWORK(LMTYPE)=5 MBAND=IWORK(LML)+IWORK(LMU)+1 MSAVE=(NEQ/MBAND)+1 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE GO TO 60 50 IWORK(LMTYPE)=4 LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD C C check lengths of rwork and iwork 60 LENIW=20+NEQ IF(LRW.LT.LENRW)GO TO 704 IF(LIW.LT.LENIW)GO TO 705 C C check to see that tout is different from t IF(TOUT .EQ. T)GO TO 719 C C check hmax IF(INFO(7).EQ.0)GO TO 70 HMAX=RWORK(LHMAX) IF(HMAX.LE.0.0D0)GO TO 710 70 CONTINUE C C initialize counters IWORK(LNST)=0 IWORK(LNRE)=0 IWORK(LNJE)=0 C IWORK(LNSTL)=0 IDID=1 GO TO 200 C C----------------------------------------------------------------------- C this block is for continuation calls C only. here we check info(1),and if the C last step was interrupted we check whether C appropriate action was taken. C----------------------------------------------------------------------- C 100 CONTINUE IF(INFO(1).EQ.1)GO TO 110 IF(INFO(1).NE.-1)GO TO 701 C if we are here, the last step was interrupted C by an error condition from ddastp,and C appropriate action was not taken. this C is a fatal error. CALL XERRWV( *49HDASSL-- THE LAST STEP TERMINATED WITH A NEGATIVE, *49,201,0,0,0,0,0,0.0D0,0.0D0) CALL XERRWV( *47HDASSL-- VALUE (=I1) OF IDID AND NO APPROPRIATE, *47,202,0,1,IDID,0,0,0.0D0,0.0D0) CALL XERRWV( *41HDASSL-- ACTION WAS TAKEN. RUN TERMINATED, *41,203,1,0,0,0,0,0.0D0,0.0D0) RETURN 110 CONTINUE IWORK(LNSTL)=IWORK(LNST) C C----------------------------------------------------------------------- C this block is executed on all calls. C the error tolerance parameters are C checked, and the work array pointers C are set. C----------------------------------------------------------------------- C 200 CONTINUE C check rtol,atol NZFLG=0 RTOLI=RTOL(1) ATOLI=ATOL(1) DO 210 I=1,NEQ IF(INFO(2).EQ.1)RTOLI=RTOL(I) IF(INFO(2).EQ.1)ATOLI=ATOL(I) IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1 IF(RTOLI.LT.0.0D0)GO TO 706 IF(ATOLI.LT.0.0D0)GO TO 707 210 CONTINUE IF(NZFLG.EQ.0)GO TO 708 C C set up rwork storage.iwork storage is fixed C in data statement. LE=LDELTA+NEQ LWT=LE+NEQ LPHI=LWT+NEQ LPD=LPHI+(IWORK(LMXORD)+1)*NEQ LWM=LPD NPD=1 NTEMP=NPD+LENPD IF(INFO(1).EQ.1)GO TO 400 C C----------------------------------------------------------------------- C this block is executed on the initial call C only. set the initial step size, and C the error weight vector, and phi. C compute initial yprime, if necessary. C----------------------------------------------------------------------- C 300 CONTINUE TN=T IDID=1 C C set error weight vector wt CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) DO 305 I = 1,NEQ IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713 305 CONTINUE C C compute unit roundoff and hmin UROUND = D1MACH(4) RWORK(LROUND) = UROUND HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT)) C C check initial interval to see that it is long enough TDIST = DABS(TOUT - T) IF(TDIST .LT. HMIN) GO TO 714 C C check ho, if this was input IF (INFO(8) .EQ. 0) GO TO 310 HO = RWORK(LH) IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711 IF (HO .EQ. 0.0D0) GO TO 712 GO TO 320 310 CONTINUE C C compute initial stepsize, to be used by either C ddastp or ddaini, depending on info(11) HO = 0.001D0*TDIST YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR) IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM HO = DSIGN(HO,TOUT-T) C adjust ho if necessary to meet hmax bound 320 IF (INFO(7) .EQ. 0) GO TO 330 RH = DABS(HO)/HMAX IF (RH .GT. 1.0D0) HO = HO/RH C compute tstop, if applicable 330 IF (INFO(4) .EQ. 0) GO TO 340 TSTOP = RWORK(LTSTOP) IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715 IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709 C C compute initial derivative, if applicable 340 IF (INFO(11) .EQ. 0) GO TO 350 CALL DDAINI(T,Y,YPRIME,NEQ, * RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR, * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), * RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),INFO(10)) IF(ITRACE .GE. 1)WRITE(IDEV,349)IDID 349 FORMAT(' IDID FROM INIT SOLVER IS ',I3) IF (IDID .LT. 0) GO TO 390 C C load h with ho. store h in rwork(lh) 350 H = HO RWORK(LH) = H C C load y and h*yprime into phi(*,1) and phi(*,2) 360 ITEMP = LPHI + NEQ DO 370 I = 1,NEQ RWORK(LPHI + I - 1) = Y(I) 370 RWORK(ITEMP + I - 1) = H*YPRIME(I) C 390 GO TO 500 C C------------------------------------------------------- C this block is for continuation calls only. its C purpose is to check stop conditions before C taking a step. C adjust h if necessary to meet hmax bound C------------------------------------------------------- C 400 CONTINUE DONE = .FALSE. TN=RWORK(LTN) H=RWORK(LH) IF(INFO(7) .EQ. 0) GO TO 410 RH = DABS(H)/HMAX IF(RH .GT. 1.0D0) H = H/RH 410 CONTINUE IF(T .EQ. TOUT) GO TO 719 IF((T - TOUT)*H .GT. 0.0D0) GO TO 711 IF(INFO(4) .EQ. 1) GO TO 430 IF(INFO(3) .EQ. 1) GO TO 420 IF((TN-TOUT)*H.LT.0.0D0)GO TO 490 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T=TOUT IDID = 3 DONE = .TRUE. GO TO 490 420 IF((TN-T)*H .LE. 0.0D0) GO TO 490 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T = TN IDID = 1 DONE = .TRUE. GO TO 490 425 CONTINUE CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T = TOUT IDID = 3 DONE = .TRUE. GO TO 490 430 IF(INFO(3) .EQ. 1) GO TO 440 TSTOP=RWORK(LTSTOP) IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715 IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709 IF((TN-TOUT)*H.LT.0.0D0)GO TO 450 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T=TOUT IDID = 3 DONE = .TRUE. GO TO 490 440 TSTOP = RWORK(LTSTOP) IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715 IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709 IF((TN-T)*H .LE. 0.0D0) GO TO 450 IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445 CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T = TN IDID = 1 DONE = .TRUE. GO TO 490 445 CONTINUE CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD), * RWORK(LPHI),RWORK(LPSI)) T = TOUT IDID = 3 DONE = .TRUE. GO TO 490 450 CONTINUE C check whether we are with in roundoff of tstop IF(DABS(TN-TSTOP).GT.100.0D0*UROUND* * (DABS(TN)+DABS(H)))GO TO 460 IDID=2 T=TSTOP DONE = .TRUE. GO TO 490 460 TNEXT=TN+H*(1.0D0+4.0D0*UROUND) IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490 H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND) RWORK(LH)=H C 490 IF (DONE) GO TO 590 C C------------------------------------------------------- C the next block contains the call to the C one-step integrator ddastp. C this is a looping point for the integration C steps. C check for too many steps. C update wt. C check for too much accuracy requested. C compute minimum stepsize. C------------------------------------------------------- C 500 CONTINUE C check for failure to compute initial yprime IF (IDID .EQ. -12) GO TO 527 C C check for too many steps IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) * GO TO 510 IDID=-1 GO TO 527 C C update wt 510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI), * RWORK(LWT),RPAR,IPAR) DO 520 I=1,NEQ IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520 IDID=-3 GO TO 527 520 CONTINUE C C test for too much accuracy requested. R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)* * 100.0D0*UROUND C IF(ITRACE.GE.1 .AND. R.GT.1.0D0)WRITE(IDEV,521)RTOL(1),ATOL(1),R, C * UROUND C521 FORMAT(' TOLS ',2D11.3,' MUST BE MULT BY R=',D11.3,' UR =',D11.3) IF(R.LE.1.0D0 )GO TO 525 C multiply rtol and atol by r and return IF(INFO(2).EQ.1)GO TO 523 RTOL(1)=R*RTOL(1) ATOL(1)=R*ATOL(1) IDID=-2 GO TO 527 523 DO 524 I=1,NEQ RTOL(I)=R*RTOL(I) 524 ATOL(I)=R*ATOL(I) IDID=-2 GO TO 527 525 CONTINUE C C compute minimum stepsize HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT)) C CALL DDASTP(TN,Y,YPRIME,NEQ, * RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR, * RWORK(LPHI),RWORK(LDELTA),RWORK(LE), * RWORK(LWM),IWORK(LIWM), * RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA), * RWORK(LPSI),RWORK(LSIGMA), * RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD), * RWORK(LS),HMIN,RWORK(LROUND), * IWORK(LPHASE),IWORK(LJCALC),IWORK(LK), * IWORK(LKOLD),IWORK(LNS),INFO(10)) 527 IF(IDID.LT.0)GO TO 600 C C------------------------------------------------------ C this block handles the case of a successful C return from ddastp (idid=1) test for C stop conditions. C-------------------------------------------------------- C IF(ITRACE .GE. 1)WRITE(IDEV,528)TN,H, IWORK(LKOLD) C IF(ITRACE .GE. 1)WRITE(IDEV,5281)TERKP1, TERK, TERKM1, TERKM2 528 FORMAT(' AT T= ',D11.3,' H=',D11.3,' ORDER=',I3) 5281 FORMAT(' ERRORS FOR DESCENDING ORDERS ARE ',4D11.3) C CALL DASMON(NEQ, TN, H, Y, YPRIME, IWORK(LKOLD)) IF(INFO(4).NE.0)GO TO 540 IF(INFO(3).NE.0)GO TO 530 IF((TN-TOUT)*H.LT.0.0D0)GO TO 500 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) IDID=3 T=TOUT GO TO 580 530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535 T=TN IDID=1 GO TO 580 535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) IDID=3 T=TOUT GO TO 580 540 IF(INFO(3).NE.0)GO TO 550 IF((TN-TOUT)*H.LT.0.0D0)GO TO 542 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) T=TOUT IDID=3 GO TO 580 542 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND* * (DABS(TN)+DABS(H)))GO TO 545 TNEXT=TN+H*(1.0D0+4.0D0*UROUND) IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500 H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND) GO TO 500 545 IDID=2 T=TSTOP GO TO 580 550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555 IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552 T=TN IDID=1 GO TO 580 552 IDID=2 T=TSTOP GO TO 580 555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ, * IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI)) T=TOUT IDID=3 580 CONTINUE C C-------------------------------------------------------- C all successful returns from ddassl are made from C this block. C-------------------------------------------------------- C 590 CONTINUE RWORK(LTN)=TN RWORK(LH)=H RETURN C C----------------------------------------------------------------------- C this block handles all unsuccessful C returns other than for illegal input. C----------------------------------------------------------------------- C 600 CONTINUE ITEMP=-IDID GO TO (610,620,630,690,690,640,650,660,670,675, * 680,685), ITEMP C C the maximum number of steps was taken before C reaching tout 610 CALL XERRWV( *38HDASSL-- AT CURRENT T (=R1) 500 STEPS, *38,610,0,0,0,0,1,TN,0.0D0) CALL XERRWV(48HDASSL-- TAKEN ON THIS CALL BEFORE REACHING TOUT, *48,611,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C too much accuracy for machine precision 620 CALL XERRWV( *47HDASSL-- AT T (=R1) TOO MUCH ACCURACY REQUESTED, *47,620,0,0,0,0,1,TN,0.0D0) CALL XERRWV( *48HDASSL-- FOR PRECISION OF MACHINE. RTOL AND ATOL, *48,621,0,0,0,0,0,0.0D0,0.0D0) CALL XERRWV( *45HDASSL-- WERE INCREASED TO APPROPRIATE VALUES, *45,622,0,0,0,0,0,0.0D0,0.0D0) C GO TO 690 C wt(i) .le. 0.0d0 for some i (not at start of problem) 630 CALL XERRWV( *38HDASSL-- AT T (=R1) SOME ELEMENT OF WT, *38,630,0,0,0,0,1,TN,0.0D0) CALL XERRWV(28HDASSL-- HAS BECOME .LE. 0.0, *28,631,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C error test failed repeatedly or with h=hmin 640 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,640,0,0,0,0,2,TN,H) CALL XERRWV( *57HDASSL-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN, *57,641,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C corrector convergence failed repeatedly or with h=hmin 650 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,650,0,0,0,0,2,TN,H) CALL XERRWV( *48HDASSL-- CORRECTOR FAILED TO CONVERGE REPEATEDLY, *48,651,0,0,0,0,0,0.0D0,0.0D0) CALL XERRWV( *28HDASSL-- OR WITH ABS(H)=HMIN, *28,652,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C the iteration matrix is singular 660 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,660,0,0,0,0,2,TN,H) CALL XERRWV( *37HDASSL-- ITERATION MATRIX IS SINGULAR, *37,661,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C corrector failure preceeded by error test failures. 670 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,670,0,0,0,0,2,TN,H) CALL XERRWV( *49HDASSL-- CORRECTOR COULD NOT CONVERGE. ALSO, THE, *49,671,0,0,0,0,0,0.0D0,0.0D0) CALL XERRWV( *38HDASSL-- ERROR TEST FAILED REPEATEDLY., *38,672,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C corrector failure because ires = -1 675 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,675,0,0,0,0,2,TN,H) CALL XERRWV( *45HDASSL-- CORRECTOR COULD NOT CONVERGE BECAUSE, *455,676,0,0,0,0,0,0.0D0,0.0D0) CALL XERRWV( *36HDASSL-- IRES WAS EQUAL TO MINUS ONE, *36,677,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C failure because ires = -2 680 CALL XERRWV( *40HDASSL-- AT T (=R1) AND STEPSIZE H (=R2), *40,680,0,0,0,0,2,TN,H) CALL XERRWV( *36HDASSL-- IRES WAS EQUAL TO MINUS TWO, *36,681,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 C C failed to compute initial yprime 685 CALL XERRWV( *44HDASSL-- AT T (=R1) AND STEPSIZE H (=R2) THE, *44,685,0,0,0,0,2,TN,HO) CALL XERRWV( *45HDASSL-- INITIAL YPRIME COULD NOT BE COMPUTED, *45,686,0,0,0,0,0,0.0D0,0.0D0) GO TO 690 690 CONTINUE INFO(1)=-1 T=TN RWORK(LTN)=TN RWORK(LH)=H RETURN C----------------------------------------------------------------------- C this block handles all error returns due C to illegal input, as detected before calling C ddastp. first the error message routine is C called. if this happens twice in C succession, execution is terminated C C----------------------------------------------------------------------- 701 CALL XERRWV( *55HDASSL-- SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE, *55,1,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 702 CALL XERRWV(25HDASSL-- NEQ (=I1) .LE. 0, *25,2,0,1,NEQ,0,0,0.0D0,0.0D0) GO TO 750 703 CALL XERRWV(34HDASSL-- MAXORD (=I1) NOT IN RANGE, *34,3,0,1,MXORD,0,0,0.0D0,0.0D0) GO TO 750 704 CALL XERRWV( *60HDASSL-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2), *60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0) GO TO 750 705 CALL XERRWV( *60HDASSL-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2), *60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0) GO TO 750 706 CALL XERRWV( *39HDASSL-- SOME ELEMENT OF RTOL IS .LT. 0, *39,6,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 707 CALL XERRWV( *39HDASSL-- SOME ELEMENT OF ATOL IS .LT. 0, *39,7,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 708 CALL XERRWV( *47HDASSL-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO, *47,8,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 709 CALL XERRWV( *54HDASSL-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2), *54,9,0,0,0,0,2,TSTOP,TOUT) GO TO 750 710 CALL XERRWV(28HDASSL-- HMAX (=R1) .LT. 0.0, *28,10,0,0,0,0,1,HMAX,0.0D0) GO TO 750 711 CALL XERRWV(34HDASSL-- TOUT (=R1) BEHIND T (=R2), *34,11,0,0,0,0,2,TOUT,T) GO TO 750 712 CALL XERRWV(29HDASSL-- INFO(8)=1 AND H0=0.0, *29,12,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 713 CALL XERRWV(39HDASSL-- SOME ELEMENT OF WT IS .LE. 0.0, *39,13,0,0,0,0,0,0.0D0,0.0D0) GO TO 750 714 CALL XERRWV( *61HDASSL-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION, *61,14,0,0,0,0,2,TOUT,T) GO TO 750 715 CALL XERRWV( *49HDASSL-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2), *49,15,0,0,0,0,2,TSTOP,T) GO TO 750 717 CALL XERRWV( *52HDASSL-- ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ, *52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0) GO TO 750 718 CALL XERRWV( *52HDASSL-- MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ, *52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0) GO TO 750 719 CALL XERRWV( *39HDASSL-- TOUT (=R1) IS EQUAL TO T (=R2), *39,19,0,0,0,0,2,TOUT,T) GO TO 750 750 IF(INFO(1).EQ.-1) GO TO 760 INFO(1)=-1 IDID=-33 RETURN 760 CALL XERRWV( *46HDASSL-- REPEATED OCCURRENCES OF ILLEGAL INPUT, *46,801,0,0,0,0,0,0.0D0,0.0D0) 770 CALL XERRWV( *47HDASSL-- RUN TERMINATED. APPARENT INFINITE LOOP, *47,802,1,0,0,0,0,0.0D0,0.0D0) RETURN C-----------end of subroutine ddassl------------------------------------ END SUBROUTINE DDAWTS(NEQ,IWT,RTOL,ATOL,Y,WT,RPAR,IPAR) C C***BEGIN PROLOGUE DDAWTS C***REFER TO DDASSL C***ROUTINES CALLED (NONE) C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDAWTS C----------------------------------------------------------------------- C this subroutine sets the error weight vector C wt according to wt(i)=rtol(i)*abs(y(i))+atol(i), C i=1,-,n. C rtol and atol are scalars if iwt = 0, C and vectors if iwt = 1. C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RTOL(1),ATOL(1),Y(1),WT(1) DIMENSION RPAR(1),IPAR(1) RTOLI=RTOL(1) ATOLI=ATOL(1) DO 20 I=1,NEQ IF (IWT .EQ.0) GO TO 10 RTOLI=RTOL(I) ATOLI=ATOL(I) 10 WT(I)=RTOLI*DABS(Y(I))+ATOLI 20 CONTINUE RETURN C-----------end of subroutine ddawts------------------------------------ END SUBROUTINE DDASTP(X,Y,YPRIME,NEQ, * RES,JAC,H,WT,JSTART,IDID,RPAR,IPAR, * PHI,DELTA,E,WM,IWM, * ALPHA,BETA,GAMMA,PSI,SIGMA, * CJ,CJOLD,HOLD,S,HMIN,UROUND, * IPHASE,JCALC,K,KOLD,NS,NONNEG) C C***BEGIN PROLOGUE DDASTP C***REFER TO DDASSL C***ROUTINES CALLED DDANRM,DDAJAC,DDASLV,DDATRP C***COMMON BLOCKS DDA001 C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDASTP C C C----------------------------------------------------------------------- C dastep solves a system of differential/ C algebraic equations of the form C g(x,y,yprime) = 0, for one step (normally C from x to x+h). C C the methods used are modified divided C difference,fixed leading coefficient C forms of backward differentiation C formulas. the code adjusts the stepsize C and order to control the local error per C step. C C C the parameters represent C x -- independent variable C y -- solution vector at x C yprime -- derivative of solution vector C after successful step C neq -- number of equations to be integrated C res -- external user-supplied subroutine C to evaluate the residual. the call is C call res(x,y,yprime,delta,ires,rpar,ipar) C x,y,yprime are input. delta is output. C on input, ires=0. res should alter ires only C if it encounters an illegal value of y or a C stop condition. set ires=-1 if an input value C of y is illegal, and dastep will try to solve C the problem without getting ires = -1. if C ires=-2, dastep returns control to the calling C program with idid = -11. C jac -- external user-supplied routine to evaluate C the iteration matrix (this is optional) C the call is of the form C call jac(x,y,yprime,pd,cj,rpar,ipar) C pd is the matrix of partial derivatives, C pd=dg/dy+cj*dg/dyprime C h -- appropriate step size for next step. C normally determined by the code C wt -- vector of weights for error criterion. C jstart -- integer variable set 0 for C first step, 1 otherwise. C idid -- completion code with the following meanings% C idid= 1 -- the step was completed successfully C idid=-6 -- the error test failed repeatedly C idid=-7 -- the corrector could not converge C idid=-8 -- the iteration matrix is singular C idid=-9 -- the corrector could not converge. C there were repeated error test C failures on this step. C idid=-10-- the corrector could not converge C because ires was equal to minus one C idid=-11-- ires equal to -2 was encountered, C and control is being returned to C the calling program C rpar,ipar -- real and integer parameter arrays that C are used for communication between the C calling program and external user routines C they are not altered by dastep C phi -- array of divided differences used by C dastep. the length is neq*(k+1),where C k is the maximum order C delta,e -- work vectors for dastep of length neq C wm,iwm -- real and integer arrays storing C matrix information such as the matrix C of partial derivatives,permutation C vector,and various other information. C C the other parameters are information C which is needed internally by dastep to C continue from step to step. C C----------------------------------------------------------------------- C C C IMPLICIT REAL*8(A-H,O-Z) LOGICAL CONVGD DIMENSION Y(1),YPRIME(1),WT(1) DIMENSION PHI(NEQ,1),DELTA(1),E(1) DIMENSION WM(1),IWM(1) DIMENSION PSI(1),ALPHA(1),BETA(1),GAMMA(1),SIGMA(1) DIMENSION RPAR(1),IPAR(1) EXTERNAL RES,JAC COMMON /SDEV2/ ITRACE, IDEV COMMON /SDDTR/ TERKP1, TERK, TERKM1, TERKM2 COMMON/DDA001/NPD,NTEMP, * LML,LMU,LMXORD,LMTYPE, * LNST,LNRE,LNJE,LETF,LCTF,LIPVT COMMON /ERRCNT/ IEFAIL DATA MAXIT/4/ DATA XRATE/0.25D0/ C C C C C C----------------------------------------------------------------------- C block 1. C initialize. on the first call,set C the order to 1 and initialize C other variables. C----------------------------------------------------------------------- C C initializations for all calls IDID=1 XOLD=X NCF=0 NSF=0 NEF=0 IF(JSTART .NE. 0) GO TO 120 C C if this is the first step,perform C other initializations IWM(LETF) = 0 IWM(LCTF) = 0 K=1 KOLD=0 HOLD=0.0D0 JSTART=1 PSI(1)=H CJOLD = 1.0D0/H CJ = CJOLD S = 100.D0 JCALC = -1 DELNRM=1.0D0 IPHASE = 0 NS=0 120 CONTINUE C C C C C C----------------------------------------------------------------------- C block 2 C compute coefficients of formulas for C this step. C----------------------------------------------------------------------- 200 CONTINUE KP1=K+1 KP2=K+2 KM1=K-1 XOLD=X IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0 NS=MIN0(NS+1,KOLD+2) NSP1=NS+1 IF(KP1 .LT. NS)GO TO 230 C BETA(1)=1.0D0 ALPHA(1)=1.0D0 TEMP1=H GAMMA(1)=0.0D0 SIGMA(1)=1.0D0 DO 210 I=2,KP1 TEMP2=PSI(I-1) PSI(I-1)=TEMP1 BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2 TEMP1=TEMP2+H ALPHA(I)=H/TEMP1 SIGMA(I)=DFLOAT(I-1)*SIGMA(I-1)*ALPHA(I) GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H 210 CONTINUE PSI(KP1)=TEMP1 230 CONTINUE C C compute alphas, alpha0 ALPHAS = 0.0D0 ALPHA0 = 0.0D0 DO 240 I = 1,K ALPHAS = ALPHAS - 1.0D0/DFLOAT(I) ALPHA0 = ALPHA0 - ALPHA(I) 240 CONTINUE C C compute leading coefficient cj CJLAST = CJ CJ = -ALPHAS/H C C compute variable stepsize error coefficient ck CK = DABS(ALPHA(KP1) + ALPHAS - ALPHA0) CK = DMAX1(CK,ALPHA(KP1)) C C decide whether new jacobian is needed TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE) TEMP2 = 1.0D0/TEMP1 IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1 IF (CJ .NE. CJLAST) S = 100.D0 C C change phi to phi star IF(KP1 .LT. NSP1) GO TO 280 DO 270 J=NSP1,KP1 DO 260 I=1,NEQ 260 PHI(I,J)=BETA(J)*PHI(I,J) 270 CONTINUE 280 CONTINUE C C update time X=X+H C C C C C C----------------------------------------------------------------------- C block 3 C predict the solution and derivative, C and solve the corrector equation C----------------------------------------------------------------------- C C first,predict the solution and derivative 300 CONTINUE DO 310 I=1,NEQ Y(I)=PHI(I,1) 310 YPRIME(I)=0.0D0 DO 330 J=2,KP1 DO 320 I=1,NEQ Y(I)=Y(I)+PHI(I,J) 320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) 330 CONTINUE PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR) C C C C solve the corrector equation using a C modified newton scheme. CONVGD= .TRUE. M=0 IWM(LNRE)=IWM(LNRE)+1 IRES = 0 CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR) IF (IRES .LT. 0) GO TO 380 C C C if indicated,reevaluate the C iteration matrix pd = dg/dy + cj*dg/dyprime C (where g(x,y,yprime)=0). set C jcalc to 0 as an indicator that C this has been done. IF(JCALC .NE. -1)GO TO 340 IWM(LNJE)=IWM(LNJE)+1 JCALC=0 CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H, * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR) IF(ITRACE .GE. 1)WRITE(IDEV,331) 331 FORMAT(' JAC EVAL') CJOLD=CJ S = 100.D0 IF (IRES .LT. 0) GO TO 380 IF(IER .NE. 0)GO TO 380 NSF=0 C C C initialize the error accumulation vector e. 340 CONTINUE DO 345 I=1,NEQ 345 E(I)=0.0D0 C S = 100.E0 C C C corrector loop. 350 CONTINUE C C multiply residual by temp1 to accelerate convergence TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD) DO 355 I = 1,NEQ 355 DELTA(I) = DELTA(I) * TEMP1 C C compute a new iterate (back-substitution). C store the correction in delta. CALL DDASLV(NEQ,DELTA,WM,IWM) C C update y,e,and yprime DO 360 I=1,NEQ Y(I)=Y(I)-DELTA(I) E(I)=E(I)-DELTA(I) 360 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) C C test for convergence of the iteration DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375 IF (M .GT. 0) GO TO 365 OLDNRM = DELNRM GO TO 367 365 RATE = (DELNRM/OLDNRM)**(1.0D0/DFLOAT(M)) IF (RATE .GT. 0.90D0) GO TO 370 S = RATE/(1.0D0 - RATE) 367 IF (S*DELNRM .LE. 0.33D0) GO TO 375 C C the corrector has not yet converged. C update m and test whether the C maximum number of iterations have C been tried. M=M+1 IF(M.GE.MAXIT)GO TO 370 C C evaluate the residual C and go back to do another iteration IWM(LNRE)=IWM(LNRE)+1 IRES = 0 CALL RES(X,Y,YPRIME,DELTA,IRES, * RPAR,IPAR) IF (IRES .LT. 0) GO TO 380 GO TO 350 C C C the corrector failed to converge in maxit C iterations. if the iteration matrix C is not current,re-do the step with C a new iteration matrix. 370 CONTINUE IF(ITRACE .GE. 1)WRITE(IDEV,371) 371 FORMAT(' CONVERGENCE FAILURE 1') IF(JCALC.EQ.0)GO TO 380 JCALC=-1 GO TO 300 C C C the iteration has converged. if nonnegativity of solution is C required, set the solution nonnegative, if the perturbation C to do it is small enough. if the change is too large, then C consider the corrector iteration to have failed. 375 IF(NONNEG .EQ. 0) GO TO 390 DO 377 I = 1,NEQ 377 DELTA(I) = DMIN1(Y(I),0.0D0) DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR) IF(DELNRM .GT. 0.33D0) GO TO 380 DO 378 I = 1,NEQ 378 E(I) = E(I) - DELTA(I) GO TO 390 C C C exits from block 3 C no convergence with current iteration C matrix,or singular iteration matrix 380 CONVGD= .FALSE. IF(ITRACE .GE. 1)WRITE(IDEV,381) 381 FORMAT(' CONVERGENCE FAILURE 2') 390 JCALC = 1 IF(.NOT.CONVGD)GO TO 600 C C C C C C----------------------------------------------------------------------- C block 4 C estimate the errors at orders k,k-1,k-2 C as if constant stepsize was used. estimate C the local error at order k and test C whether the current step is successful. C----------------------------------------------------------------------- C C estimate errors at orders k,k-1,k-2 ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR) ERK = SIGMA(K+1)*ENORM TERK = FLOAT(K+1)*ERK EST = ERK KNEW=K IF(K .EQ. 1)GO TO 430 DO 405 I = 1,NEQ 405 DELTA(I) = PHI(I,KP1) + E(I) ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) TERKM1 = FLOAT(K)*ERKM1 IF(K .GT. 2)GO TO 410 IF(TERKM1 .LE. 0.5*TERK)GO TO 420 GO TO 430 410 CONTINUE DO 415 I = 1,NEQ 415 DELTA(I) = PHI(I,K) + DELTA(I) ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) TERKM2 = FLOAT(K-1)*ERKM2 IF(DMAX1(TERKM1,TERKM2).GT.TERK)GO TO 430 C lower the order 420 CONTINUE KNEW=K-1 EST = ERKM1 C C C calculate the local error for the current step C to see if the step was successful 430 CONTINUE ERR = CK * ENORM IF(ITRACE .GE. 1)WRITE(IDEV,431)ERR 431 FORMAT(' SCALED LOCAL ERROR IS ',D12.3) IF(ERR .GT. 1.0D0)GO TO 600 C C C C C C----------------------------------------------------------------------- C block 5 C the step is successful. determine C the best order and stepsize for C the next step. update the differences C for the next step. C----------------------------------------------------------------------- IDID=1 IWM(LNST)=IWM(LNST)+1 KDIFF=K-KOLD KOLD=K HOLD=H C C C estimate the error at order k+1 unless% C already decided to lower order, or C already using maximum order, or C stepsize not constant, or C order raised in previous step IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1 IF(IPHASE .EQ. 0)GO TO 545 IF(KNEW.EQ.KM1)GO TO 540 IF(K.EQ.IWM(LMXORD)) GO TO 550 IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550 DO 510 I=1,NEQ 510 DELTA(I)=E(I)-PHI(I,KP2) ERKP1 = (1.0D0/DFLOAT(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR) TERKP1 = FLOAT(K+2)*ERKP1 IF(K.GT.1)GO TO 520 IF(TERKP1.GE.0.5D0*TERK)GO TO 550 GO TO 530 520 IF(TERKM1.LE.DMIN1(TERK,TERKP1))GO TO 540 IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550 TEMP2=K+1 R= (2.0D0*ERK +0.0001D0)**(-1.0D0/TEMP2) 1 / (2.0D0*ERKP1+0.0001D0)**(-1.0D0/(TEMP2+1.D0)) IF(R .GE. 0.9D0)GOTO 550 C C raise order 530 K=KP1 IF(ITRACE .GE. 1)WRITE(IDEV,531) 531 FORMAT(' ORDER RAISE CONSIDERED') EST = ERKP1 GO TO 550 C C lower order 540 K=KM1 EST = ERKM1 GO TO 550 C C if iphase = 0, increase order by one and multiply stepsize by C factor two 545 K = KP1 IF(ITRACE .GE. 1)WRITE(IDEV,546) 546 FORMAT(' ORDER RAISE WITH IPHASE =0') HNEW = H*2.0D0 H = HNEW GO TO 575 C C C determine the appropriate stepsize for C the next step. 550 HNEW=H TEMP2=K+1 R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) IF(R .LT. 2.0D0) GO TO 555 HNEW = 2.0D0*H GO TO 560 555 IF(R .GT. 1.0D0) GO TO 560 R = DMAX1(0.5D0,DMIN1(0.9D0,R)) HNEW = H*R 560 H=HNEW C C C update differences for next step 575 CONTINUE IF(KOLD.EQ.IWM(LMXORD))GO TO 585 DO 580 I=1,NEQ 580 PHI(I,KP2)=E(I) 585 CONTINUE DO 590 I=1,NEQ 590 PHI(I,KP1)=PHI(I,KP1)+E(I) DO 595 J1=2,KP1 J=KP1-J1+1 DO 595 I=1,NEQ 595 PHI(I,J)=PHI(I,J)+PHI(I,J+1) RETURN C C C C C C----------------------------------------------------------------------- C block 6 C the step is unsuccessful. restore x,psi,phi C determine appropriate stepsize for C continuing the integration, or exit with C an error flag if there have been many C failures. C----------------------------------------------------------------------- 600 IPHASE = 1 C C restore x,phi,psi X=XOLD IF(KP1.LT.NSP1)GO TO 630 DO 620 J=NSP1,KP1 TEMP1=1.0D0/BETA(J) DO 610 I=1,NEQ 610 PHI(I,J)=TEMP1*PHI(I,J) 620 CONTINUE 630 CONTINUE DO 640 I=2,KP1 640 PSI(I-1)=PSI(I)-H C C C test whether failure is due to corrector iteration C or error test IF(CONVGD)GO TO 660 IWM(LCTF)=IWM(LCTF)+1 C C C the newton iteration failed to converge with C a current iteration matrix. determine the cause C of the failure and take appropriate action. IF(IER.EQ.0)GO TO 650 C C the iteration matrix is singular. reduce C the stepsize by a factor of 4. if C this happens three times in a row on C the same step, return with an error flag NSF=NSF+1 R = 0.25D0 H=H*R IF (NSF .LT. 3 .AND. DABS(H) .GE. HMIN) GO TO 690 IDID=-8 GO TO 675 C C C the newton iteration failed to converge for a reason C other than a singular iteration matrix. if ires = -2, then C return. otherwise, reduce the stepsize and try again, unless C too many failures have occured. 650 CONTINUE IF (IRES .GT. -2) GO TO 655 IDID = -11 GO TO 675 655 NCF = NCF + 1 R = 0.25D0 H = H*R IF (NCF .LT. 10 .AND. DABS(H) .GE. HMIN) GO TO 690 IDID = -7 IF (IRES .LT. 0) IDID = -10 IF (NEF .GE. 3) IDID = -9 GO TO 675 C C C the newton scheme converged,and the cause C of the failure was the error estimate C exceeding the tolerance. 660 NEF=NEF+1 IEFAIL = IEFAIL + 1 IF(ITRACE .GE.1)WRITE(IDEV,661) 661 FORMAT(' ERROR TEST FAILED') IWM(LETF)=IWM(LETF)+1 IF (NEF .GT. 1) GO TO 665 C C on first error test failure, keep current order or lower C order by one. compute new stepsize based on differences C of the solution. K = KNEW TEMP2 = K + 1 R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) R = DMAX1(0.25D0,DMIN1(0.9D0,R)) H = H*R IF (DABS(H) .GE. HMIN) GO TO 690 IDID = -6 GO TO 675 C C on second error test failure, use the current order or C decrease order by one. reduce the stepsize by a factor of C one quarter. 665 IF (NEF .GT. 2) GO TO 670 K = KNEW H = 0.25D0*H IF (DABS(H) .GE. HMIN) GO TO 690 IDID = -6 GO TO 675 C C on third and subsequent error test failures, set the order to C one and reduce the stepsize by a factor of one quarter 670 K = 1 H = 0.25D0*H IF (DABS(H) .GE. HMIN) GO TO 690 IDID = -6 GO TO 675 C C C C C for all crashes, restore y to its last value, C interpolate to find yprime at last x, and return 675 CONTINUE CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI) RETURN C C C go back and try this step again 690 GO TO 200 C C------end of subroutine dastep------ END SUBROUTINE DDASLV(NEQ,DELTA,WM,IWM) C C***BEGIN PROLOGUE DDASLV C***REFER TO DDASSL C***ROUTINES CALLED DGESL,DGBSL C***COMMON BLOCKS DDA001 C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDASLV C----------------------------------------------------------------------- C this routine manages the solution of the linear C system arising in the newton iteration. C matrices and real temporary storage and C real information are stored in the array wm. C integer matrix information is stored in C the array iwm. C for a dense matrix, the linpack routine C dgesl is called. C for a banded matrix,the linpack routine C dgbsl is called. C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DELTA(1),WM(1),IWM(1) COMMON/DDA001/NPD,NTEMP,LML,LMU, * LMXORD,LMTYPE, * LNST,LNRE,LNJE,LETF,LCTF,LIPVT COMMON /PSTATS/ ID(3), IS IS = IS + 1 C MTYPE=IWM(LMTYPE) GO TO(100,100,300,400,400),MTYPE C C dense matrix 100 CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0) RETURN C C dummy section for mtype=3 300 CONTINUE RETURN C C banded matrix 400 MEBAND=2*IWM(LML)+IWM(LMU)+1 CALL DGBSL(WM(NPD),MEBAND,NEQ,IWM(LML), * IWM(LMU),IWM(LIPVT),DELTA,0) RETURN C------end of subroutine ddaslv------ END SUBROUTINE DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H, * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR) C C***BEGIN PROLOGUE DDAJAC C***REFER TO DDASSL C***ROUTINES CALLED DGEFA,DGBFA C***COMMON BLOCKS DDA001 C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDAJAC C----------------------------------------------------------------------- C this routine computes the iteration matrix C pd=dg/dy+cj*dg/dyprime (where g(x,y,yprime)=0). C here pd is computed by the user-supplied C routine jac if iwm(mtype) is 1 or 4, and C it is computed by numerical finite differencing C if iwm(mtype)is 2 or 5 C the parameters have the following meanings. C y = array containing predicted values C yprime = array containing predicted derivatives C delta = residual evaluated at (x,y,yprime) C (used only if iwm(mtype)=2 or 5) C cj = scalar parameter defining iteration matrix C h = current stepsize in integration C ier = variable which is .ne. 0 C if iteration matrix is singular, C and 0 otherwise. C wt = vector of weights for computing norms C e = work space (temporary) of length neq C wm = real work space for matrices. on C output it contains the lu decomposition C of the iteration matrix. C iwm = integer work space containing C matrix information C res = name of the external user-supplied routine C to evaluate the residual function g(x,y,yprime) C ires = flag which is equal to zero if no illegal values C in res, and less than zero otherwise. (if ires C is less than zero, the matrix was not completed) C in this case (if ires .lt. 0), then ier = 0. C uround = the unit roundoff error of the machine being used. C jac = name of the external user-supplied routine C to evaluate the iteration matrix (this routine C is only used if iwm(mtype) is 1 or 4) C----------------------------------------------------------------------- C IMPLICIT REAL*8(A-H,O-Z) EXTERNAL RES,JAC DIMENSION Y(1),YPRIME(1),DELTA(1),WT(1),E(1) DIMENSION WM(1),IWM(1),RPAR(1),IPAR(1) COMMON/DDA001/NPD,NTEMP, * LML,LMU,LMXORD,LMTYPE, * LNST,LNRE,LNJE,LETF,LCTF,LIPVT C IER = 0 NPDM1=NPD-1 MTYPE=IWM(LMTYPE) GO TO (100,200,300,400,500),MTYPE C C C dense user-supplied matrix 100 LENPD=NEQ*NEQ DO 110 I=1,LENPD 110 WM(NPDM1+I)=0.0D0 CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR) GO TO 230 C C C dense finite-difference-generated matrix 200 IRES=0 NROW=NPDM1 SQUR = DSQRT(UROUND) DO 210 I=1,NEQ DEL=SQUR*DMAX1(DABS(Y(I)),DABS(H*YPRIME(I)), * DABS(WT(I))) DEL=DSIGN(DEL,H*YPRIME(I)) DEL=(Y(I)+DEL)-Y(I) YSAVE=Y(I) YPSAVE=YPRIME(I) Y(I)=Y(I)+DEL YPRIME(I)=YPRIME(I)+CJ*DEL CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR) IF (IRES .LT. 0) RETURN DELINV=1.0D0/DEL DO 220 L=1,NEQ 220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV NROW=NROW+NEQ Y(I)=YSAVE YPRIME(I)=YPSAVE 210 CONTINUE C C C do dense-matrix lu decomposition on pd 230 CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER) RETURN C C C dummy section for iwm(mtype)=3 300 RETURN C C C banded user-supplied matrix 400 LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ DO 410 I=1,LENPD 410 WM(NPDM1+I)=0.0D0 CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR) MEBAND=2*IWM(LML)+IWM(LMU)+1 GO TO 550 C C C banded finite-difference-generated matrix 500 MBAND=IWM(LML)+IWM(LMU)+1 MBA=MIN0(MBAND,NEQ) MEBAND=MBAND+IWM(LML) MEB1=MEBAND-1 MSAVE=(NEQ/MBAND)+1 ISAVE=NTEMP-1 IPSAVE=ISAVE+MSAVE IRES=0 SQUR=DSQRT(UROUND) DO 540 J=1,MBA DO 510 N=J,NEQ,MBAND K= (N-J)/MBAND + 1 WM(ISAVE+K)=Y(N) WM(IPSAVE+K)=YPRIME(N) DEL=SQUR*DMAX1(DABS(Y(N)),DABS(H*YPRIME(N)), * DABS(WT(N))) DEL=DSIGN(DEL,H*YPRIME(N)) DEL=(Y(N)+DEL)-Y(N) Y(N)=Y(N)+DEL 510 YPRIME(N)=YPRIME(N)+CJ*DEL CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR) IF (IRES .LT. 0) RETURN DO 530 N=J,NEQ,MBAND K= (N-J)/MBAND + 1 Y(N)=WM(ISAVE+K) YPRIME(N)=WM(IPSAVE+K) DEL=SQUR*DMAX1(DABS(Y(N)),DABS(H*YPRIME(N)), * DABS(WT(N))) DEL=DSIGN(DEL,H*YPRIME(N)) DEL=(Y(N)+DEL)-Y(N) DELINV=1.0D0/DEL I1=MAX0(1,(N-IWM(LMU))) I2=MIN0(NEQ,(N+IWM(LML))) II=N*MEB1-IWM(LML)+NPDM1 DO 520 I=I1,I2 520 WM(II+I)=(E(I)-DELTA(I))*DELINV 530 CONTINUE 540 CONTINUE C C C do lu decomposition of banded pd 550 CALL DGBFA(WM(NPD),MEBAND,NEQ, * IWM(LML),IWM(LMU),IWM(LIPVT),IER) RETURN C------end of subroutine ddajac------ END SUBROUTINE DDATRP(X,XOUT,YOUT,YPOUT,NEQ,KOLD,PHI,PSI) C C***BEGIN PROLOGUE DDATRP C***REFER TO DDASSL C***ROUTINES CALLED (NONE) C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDATRP C C----------------------------------------------------------------------- C the methods in subroutine dastep use polynomials C to approximate the solution. ddatrp approximates the C solution and its derivative at time xout by evaluating C one of these polynomials,and its derivative,there. C information defining this polynomial is passed from C dastep, so ddatrp cannot be used alone. C C the parameters are% C x the current time in the integration. C xout the time at which the solution is desired C yout the interpolated approximation to y at xout C (this is output) C ypout the interpolated approximation to yprime at xout C (this is output) C neq number of equations C kold order used on last successful step C phi array of scaled divided differences of y C psi array of past stepsize history C----------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION YOUT(1),YPOUT(1) DIMENSION PHI(NEQ,1),PSI(1) KOLDP1=KOLD+1 TEMP1=XOUT-X DO 10 I=1,NEQ YOUT(I)=PHI(I,1) 10 YPOUT(I)=0.0D0 C=1.0D0 D=0.0D0 GAMMA=TEMP1/PSI(1) DO 30 J=2,KOLDP1 D=D*GAMMA+C/PSI(J-1) C=C*GAMMA GAMMA=(TEMP1+PSI(J-1))/PSI(J) DO 20 I=1,NEQ YOUT(I)=YOUT(I)+C*PHI(I,J) 20 YPOUT(I)=YPOUT(I)+D*PHI(I,J) 30 CONTINUE RETURN C C------end of subroutine ddatrp------ END SUBROUTINE DDAINI(X,Y,YPRIME,NEQ, * RES,JAC,H,WT,IDID,RPAR,IPAR, * PHI,DELTA,E,WM,IWM, * HMIN,UROUND,NONNEG) C C***BEGIN PROLOGUE DDAINI C***REFER TO DDASSL C***ROUTINES CALLED DDANRM,DDAJAC,DDASLV C***COMMON BLOCKS DDA001 C***DATE WRITTEN 830315 (YYMMDD) C***REVISION DATE 830315 (YYMMDD) C***END PROLOGUE DDAINI C C------------------------------------------------------- C ddaini takes one step of size h or smaller C with the backward euler method, to C find yprime at the initial time x. a modified C damped newton iteration is used to C solve the corrector iteration. C C the initial guess yprime is used in the C prediction, and in forming the iteration C matrix, but is not involved in the C error test. this may have trouble C converging if the initial guess is no C good, or if g(xy,yprime) depends C nonlinearly on yprime. C C the parameters represent: C x -- independent variable C y -- solution vector at x C yprime -- derivative of solution vector C neq -- number of equations C h -- stepsize. imder may use a stepsize C smaller than h. C wt -- vector of weights for error C criterion C idid -- completion code with the following meanings C idid= 1 -- yprime was found successfully C idid=-12 -- ddaini failed to find yprime C rpar,ipar -- real and integer parameter arrays C that are not altered by ddaini C phi -- work space for ddaini C delta,e -- work space for ddaini C wm,iwm -- real and integer arrays storing C matrix information C C----------------------------------------------------------------- C C C IMPLICIT REAL*8 (A-H,O-Z) LOGICAL CONVGD DIMENSION Y(1),YPRIME(1),WT(1) DIMENSION PHI(NEQ,1),DELTA(1),E(1) DIMENSION WM(1),IWM(1) DIMENSION RPAR(1),IPAR(1) EXTERNAL RES,JAC COMMON/SDEV2/ ITRACE, IDEV COMMON/DDA001/NPD,NTEMP, * LML,LMU,LMXORD,LMTYPE, * LNST,LNRE,LNJE,LETF,LCTF,LIPVT C DATA MAXIT/12/,MJAC/8/ DATA DAMP/0.75D0/ C C C C--------------------------------------------------- C block 1. C initializations. C--------------------------------------------------- C IDID=1 NEF=0 NCF=0 NSF=0 YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR) C C save y and yprime in phi DO 100 I=1,NEQ PHI(I,1)=Y(I) 100 PHI(I,2)=YPRIME(I) C C C C---------------------------------------------------- C block 2. C do one backward euler step. C---------------------------------------------------- C C set up for start of corrector iteration 200 CJ=1.0D0/H XNEW=X+H C C predict solution and derivative C DO 250 I=1,NEQ 250 Y(I)=Y(I)+H*YPRIME(I) C JCALC=-1 M=0 CONVGD=.TRUE. C C C corrector loop. 300 IWM(LNRE)=IWM(LNRE)+1 IRES=0 C CALL RES(XNEW,Y,YPRIME,DELTA,IRES,RPAR,IPAR) IF (IRES.LT.0) GO TO 430 C C C evaluate the iteration matrix IF (JCALC.NE.-1) GO TO 310 IWM(LNJE)=IWM(LNJE)+1 JCALC=0 CALL DDAJAC(NEQ,XNEW,Y,YPRIME,DELTA,CJ,H, * IER,WT,E,WM,IWM,RES,IRES, * UROUND,JAC,RPAR,IPAR) IF(ITRACE .GE. 1)WRITE(IDEV,301) 301 FORMAT(' JAC EVAL IN DDAINI ') C S=1000000.D0 IF (IRES.LT.0) GO TO 430 IF (IER.NE.0) GO TO 430 NSF=0 C C C C C multiply residual by damping factor 310 CONTINUE DO 320 I=1,NEQ 320 DELTA(I)=DELTA(I)*DAMP C C C compute a new iterate (back substitution) C store the correction in delta C CALL DDASLV(NEQ,DELTA,WM,IWM) C C C update y and yprime C DO 330 I=1,NEQ Y(I)=Y(I)-DELTA(I) 330 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) C C C test for convergence of the iteration. C DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) IF (DELNRM.LE.100.D0*UROUND*YNORM) * GO TO 400 C IF (M.GT.0) GO TO 340 OLDNRM=DELNRM GO TO 350 C 340 RATE=(DELNRM/OLDNRM)**(1.0D0/DFLOAT(M)) IF (RATE.GT.0.90D0) GO TO 430 S=RATE/(1.0D0-RATE) IF(ITRACE .GE. 2)WRITE(IDEV,341)RATE, DELNRM, S 341 FORMAT(' RATE= ',D11.3,' DELNRM=',D11.3,' S = ',D11.3) C 350 IF (S*DELNRM .LE. 0.33D0) GO TO 400 C C C the corrector has not yet converged. update C m and and test whether the maximum C number of iterations have been tried. C every mjac iterations, get a new C iteration matrix. C M=M+1 IF (M.GE.MAXIT) GO TO 430 C IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1 C GO TO 300 C C C C the iteration has converged. C check nonnegativity constraints 400 IF (NONNEG.EQ.0) GO TO 450 DO 410 I=1,NEQ 410 DELTA(I)=DMIN1(Y(I),0.0D0) C DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) IF (DELNRM.GT.0.33D0) GO TO 430 C DO 420 I=1,NEQ Y(I)=Y(I)-DELTA(I) 420 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) GO TO 450 C C C exits from corrector loop. 430 CONVGD=.FALSE. 450 IF (.NOT.CONVGD) GO TO 600 C C C C----------------------------------------------------- C block 3. C the corrector iteration converged. C do error test. C----------------------------------------------------- C C DO 510 I=1,NEQ E(I)=Y(I)-PHI(I,1) 510 CONTINUE C ERR=DDANRM(NEQ,E,WT,RPAR,IPAR) C IF (ERR.LE.1.0D0) RETURN IF(ITRACE .GE. 1)WRITE(IDEV,511)ERR 511 FORMAT(' SCALED LOCAL ERROR IN INITIAL STEP IS',D11.3) DO 512 I = 1,NEQ 512 YPRIME(I) = 0.0D0 RETURN C C-------------------------------------------------------- C block 4. C the backward euler step failed. restore y C and yprime to their original values. C reduc