C C Main program to execute DDARDM demonstration program. C C***BEGIN PROLOGUE DDARDM C***DATE WRITTEN 900214 (YYMMDD) C***REVISION DATE 900214 (YYMMDD) C***CATEGORY NO. D2A2 C***KEYWORDS LIBRARY=SLATEC,QUICK CHECK,DDARDM, C TYPE=DOUBLE C***AUTHOR PETZOLD,LINDA R.,COMPUTING AND MATHEMATICS RESEARCH DIVISION C LAWRENCE LIVERMORE NATIONAL LABORATORY C L - 316, P.O. BOX 808, C LIVERMORE, CA. 94550 C***PURPOSE Quick check for SLATEC routine DDASRT. c***DESCRIPTION C Demonstration program for DDASRT. C This version is in double precision. C C DDASRT is used to solve two simple problems, C one nonstiff and one intermittently stiff. C If the errors are too large, or other difficulty occurs, C a warning message is printed. All output is on unit LUN. C C To run the demonstration problems with full printing, C set KPRINT = 3. C C***REFERENCES (NONE) C***ROUTINES CALLED DDASRT,RES1,GR1,RES2,JAC2,GR2 C***END PROLOGUE DDARDM C C----------------------------------------------------------------------- C IMPLICIT DOUBLE PRECISION (A-H,O-Z) EXTERNAL RES1, GR1, RES2, JAC2, GR2 INTEGER IDID, IOUT, IWORK, JROOT, INFO, JTYPE, * KROOT, LENIW, LENRW, LIW, LRW, LUN, NEQ, NERR, NG, * NRE, NREA, NGE, NJE, NST DOUBLE PRECISION ATOL, ER, ERO, ERRT, RTOL, RWORK, * T, TOUT, TZERO, Y, YT DIMENSION Y(2), YPRIME(2), RTOL(2), ATOL(2), * RWORK(100), IWORK(100), JROOT(2), INFO(15) COMMON /LOCAL/NEQ C LUN = 6 KPRINT = 3 IPASS = 1 NERR = 0 C----------------------------------------------------------------------- C First problem. C The initial value problem is.. C DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y, Y(1) = 1, 1 .LE. T .LE. 6 C The solution is Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3 C The two root functions are.. C G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT) (with root at T = 2.5), C G2 = LOG(Y) - 2.2491 (with roots at T = 2.47 and 2.53) C----------------------------------------------------------------------- C SET ALL INPUT PARAMETERS AND PRINT HEADING. DO 10 I = 1,15 10 INFO(I) = 0 NEQ = 1 T = 1.0D0 Y(1) = 1.0D0 TOUT = 2.0D0 RTOL(1) = 0.0D0 ATOL(1) = 1.0D-6 LRW = 100 LIW = 100 IDID = 0 C C SET INFO(11) TO ONE IF DDASRT IS TO COMPUTE THE INITIAL C YPRIME, AND GENERATE AN INITIAL GUESS FOR YPRIME. C OTHERWISE, SET INFO(11) = 0, AND SUPPLY THE CORRECT INITIAL C VALUE FOR YPRIME. C YPRIME(1) = 3.0D0 INFO(11) = 0 C C NOTE: JTYPE INDICATES THE JACOBIAN TYPE: C JTYPE = 1 ==> JACOBIAN IS DENSE AND USER-SUPPLIED C JTYPE = 2 ==> JACOBIAN IS DENSE AND COMPUTED INTERNALLY C JTYPE = 2 INFO(5) = 2 - JTYPE NG = 2 IF (KPRINT .GE. 2) THEN WRITE (LUN,110) RTOL(1),ATOL(1),JTYPE ENDIF 110 FORMAT('1'/1X,1P,' DEMONSTRATION PROGRAM FOR DDASRT'///// 1 1X,' PROBLEM 1..'/// 2 1X,' problem is DY/DT = ((2*LOG(Y)+8)/T - 5)*Y, Y(1) = 1'// 3 1X,' solution is Y(T) = EXP(-T**2 + 5*T - 4)'// 4 1X,' root functions are..'/ 5 10X,' G1 = DY/DT (root at T = 2.5).'/ 6 10X,' G2 = LOG(Y) - 6.2491 (roots at T = 2.47 and T = 2.53).'// 7 1X,' RTOL =',E10.1,' ATOL =',E10.1/ 8 1X,' JTYPE =',I3/////) C C CALL DDASRT IN LOOP OVER TOUT VALUES 2,3,4,5,6 ERO = 0.0D0 DO 180 IOUT = 1,5 120 CONTINUE C CALL DDASRT (RES1,NEQ,T,Y,YPRIME,TOUT, * INFO,RTOL,ATOL,IDID, * RWORK,LRW,IWORK,LIW,RPAR,IPAR,JDUM, * GR1,NG,JROOT) C C PRINT Y AND ERROR IN Y, AND PRINT WARNING IF ERROR TOO LARGE. YT = EXP(-T*T + 5.0D0*T - 4.0D0) ER = Y(1) - YT IF (KPRINT .GT. 2) THEN WRITE (LUN,130) T,Y(1),ER ENDIF 130 FORMAT(1X,1P,' at t =',E15.7,5X,'y =',E15.7,5X,'error =',E12.4) IF (IDID .LT. 0) GO TO 185 ER = ABS(ER)/ATOL(1) ERO = DMAX1(ERO,ER) IF (ER .LT. 1000.0D0) GO TO 140 IF (KPRINT .GE. 2) THEN IPASS = 0 WRITE (LUN,135) ENDIF 135 FORMAT(//1X,' WARNING.. ERROR EXCEEDS 1000 * TOLERANCE'//) NERR = NERR + 1 140 CONTINUE IF (IDID .NE. 4) GO TO 175 C C IF A ROOT WAS FOUND, WRITE RESULTS AND CHECK ROOT LOCATION. C THEN RETURN TO DDASRT TO CONTINUE THE INTEGRATION. IF (KPRINT .GT. 2) THEN WRITE (LUN,150) T,JROOT(1),JROOT(2) ENDIF 150 FORMAT(/1X,1P,' root found at t =',E15.7,5X,'jroot =',2I5) IF (JROOT(1) .EQ. 1) ERRT = T - 2.5D0 IF (JROOT(2) .EQ. 1 .AND. T .LE. 2.5D0) ERRT = T - 2.47D0 IF (JROOT(2) .EQ. 1 .AND. T .GT. 2.5D0) ERRT = T - 2.53D0 IF (KPRINT .GT. 2) THEN WRITE (LUN,160) ERRT ENDIF 160 FORMAT(1X,1P,' error in t location of root is',E12.4//) IF (ABS(ERRT) .LT. 1.0D-3) GO TO 170 IF (KPRINT .GE. 2) THEN IPASS = 0 WRITE (LUN,165) ENDIF 165 FORMAT(//1X,' WARNING.. ROOT ERROR EXCEEDS 1.0D-3'//) NERR = NERR + 1 170 CONTINUE GO TO 120 C C IF NO ROOT FOUND, INCREMENT TOUT AND LOOP BACK. 175 TOUT = TOUT + 1.0D0 180 CONTINUE C C PROBLEM COMPLETE. PRINT FINAL STATISTICS. 185 CONTINUE IF (IDID .LT. 0) NERR = NERR + 1 NST = IWORK(11) NJE = IWORK(13) NRE = IWORK(12) NGE = IWORK(16) LENRW = 0 LENIW = 0 NREA = NRE IF (JTYPE .EQ. 2) NRE = NRE + NEQ*NJE C IF (KPRINT .GT. 2) THEN WRITE (LUN,190) NST,NRE,NREA,NJE,NGE,ERO ENDIF 190 FORMAT(//1X,1P,' FINAL STATISTICS FOR THIS RUN..'/ * 1X,' number of steps =',I5/ * 1X,' number of r-s =',I5/ * 1X,' (excluding j-s) =',I5/ * 1X,' number of j-s =',I5/ * 1X,' number of g-s =',I5/ * 1X,' error overrun =',E10.2) C C----------------------------------------------------------------------- C Second problem (Van Der Pol oscillator). C The initial value problem is.. C DY1/DT = Y2, DY2/DT = 100*(1 - Y1**2)*Y2 - Y1, C Y1(0) = 2, Y2(0) = 0, 0 .LE. T .LE. 200 C Y1PRIME(0) = 0, Y2PRIME(0) = -2 C The root function is G = Y1. C An analytic solution is not known, but the zeros of Y1 are known C to 15 figures for purposes of checking the accuracy. C----------------------------------------------------------------------- C C RESET INFO ARRAY C DO 195 I = 1,15 195 INFO(I) = 0 C C SET TOLERANCE PARAMETERS AND PRINT HEADING. C NOTE THAT INFO(2) IS SET TO ONE, INDICATING THAT RTOL AND ATOL C ARE ARRAYS. EACH ENTRY OF RTOL AND ATOL MUST THEN BE DEFINED. C INFO(2) = 1 RTOL(1) = 1.0D-6 RTOL(2) = 1.0D-6 ATOL(1) = 1.0D-6 ATOL(2) = 1.0D-4 IF (KPRINT .GE. 2) THEN WRITE (LUN,200) RTOL(1),ATOL(1),ATOL(2) ENDIF 200 FORMAT('1'/1X,1p,' DEMONSTRATION PROGRAM FOR DDASRT'/// * 1X,' PROBLEM 2.. (Van Der Pol oscillator)'/// * 1X,' problem is DY1/DT = Y2, DY2/DT = 100*(1-Y1**2)*Y2 - Y1'/ * 1X,' Y1(0) = 2, Y2(0) = 0'// * 1X,' root function is G = Y1'// * 1X,' RTOL =',E10.1,' ATOL =',2E10.1///) C C NOTE: JTYPE INDICATES THE JACOBIAN TYPE: C JTYPE = 1 ==> JACOBIAN IS DENSE AND USER-SUPPLIED C JTYPE = 2 ==> JACOBIAN IS DENSE AND COMPUTED INTERNALLY C C LOOP OVER JTYPE = 1, 2. SET REMAINING PARAMETERS AND PRINT JTYPE. DO 290 JTYPE = 1,2 C C SET INFO(1) = 0 TO INDICATE START OF A NEW PROBLEM C SET INFO(5) = 2-JTYPE TO TELL DASSL THE JACOBIAN TYPE. C INFO(1) = 0 INFO(5) = 2-JTYPE NEQ = 2 T = 0.0D0 Y(1) = 2.0D0 Y(2) = 0.0D0 YPRIME(1) = 0.D0 YPRIME(2) = -2.0D0 TOUT = 20.0D0 LRW = 100 LIW = 100 NG = 1 IF (KPRINT .GT. 2) THEN WRITE (LUN,210) JTYPE ENDIF 210 FORMAT(////////1X,' SOLUTION WITH JTYPE =',I2////) C C CALL DDASRT IN LOOP OVER TOUT VALUES 20,40,...,200. DO 270 IOUT = 1,10 220 CONTINUE C CALL DDASRT (RES2,NEQ,T,Y,YPRIME,TOUT, * INFO,RTOL,ATOL,IDID, * RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC2, * GR2,NG,JROOT) C C PRINT Y1 AND Y2. IF (KPRINT .GT. 2) THEN WRITE (LUN,230) T,Y(1),Y(2) ENDIF 230 FORMAT(1X,1P,' at t =',E15.7,5X,'y1 =',E15.7,5X,'y2 =',E15.7) IF (IDID .LT. 0) GO TO 275 IF (IDID .NE. 4) GO TO 265 C C IF A ROOT WAS FOUND, WRITE RESULTS AND CHECK ROOT LOCATION. C THEN RETURN TO DDASRT TO CONTINUE THE INTEGRATION. IF (KPRINT .GT. 2) THEN WRITE (LUN,240) T ENDIF 240 FORMAT(/1X,1P,' root found at t =',E15.7) KROOT = INT(T/81.2D0 + 0.5D0) TZERO = 81.17237787055D0 + DFLOAT(KROOT-1)*81.41853556212D0 ERRT = T - TZERO IF (KPRINT .GT. 2) THEN WRITE (LUN,250) ERRT ENDIF 250 FORMAT(1X,1P,' error in t location of root is',E12.4//) IF (ERRT .LT. 1.0D0) GO TO 260 IF (KPRINT .GE. 2) THEN IPASS = 0 WRITE (LUN,255) ENDIF 255 FORMAT(//1X,' WARNING.. ROOT ERROR EXCEEDS 1.0'//) NERR = NERR + 1 260 CONTINUE GO TO 220 C C IF NO ROOT FOUND, INCREMENT TOUT AND LOOP BACK. 265 TOUT = TOUT + 20.0D0 270 CONTINUE C C PROBLEM COMPLETE. PRINT FINAL STATISTICS. 275 CONTINUE IF (IDID .LT. 0) NERR = NERR + 1 NST = IWORK(11) NJE = IWORK(13) NRE = IWORK(12) NGE = IWORK(16) LENRW = 0 LENIW = 0 NREA = NRE IF (JTYPE .EQ. 2) NRE = NRE + NEQ*NJE IF (KPRINT .GE. 2) THEN WRITE (LUN,280) NST,NRE,NREA,NJE,NGE ENDIF 280 FORMAT(//1X,1P,' FINAL STATISTICS FOR THIS RUN..'/ * 1X,' number of steps =',I5/ * 1X,' number of r-s =',I5/ * 1X,' (excluding j-s) =',I5/ * 1X,' number of j-s =',I5/ * 1X,' number of g-s =',I5) 290 CONTINUE C C IF (KPRINT .GE. 2) THEN WRITE (LUN,300) NERR ENDIF 300 FORMAT(////1X,' NUMBER OF ERRORS ENCOUNTERED =',I3) C IF (NERR .GT. 0) THEN IPASS = 0 ENDIF IF ((IPASS .EQ. 1) .AND. (KPRINT .GT. 1)) WRITE (LUN,700) IF ((IPASS .EQ. 0) .AND. (KPRINT .NE. 0)) WRITE (LUN,800) 700 FORMAT (/,' **********DDASRT PASSED ALL TESTS**********') 800 FORMAT (/,' **********DDASRT FAILED SOME TESTS*********') STOP END C SUBROUTINE RES1(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /LOCAL/NEQ DIMENSION Y(NEQ), YPRIME(NEQ), DELTA(NEQ) C C CHECK Y TO MAKE SURE THAT IT IS VALID INPUT. C IF Y IS LESS THAN OR EQUAL TO ZERO, THIS IS INVALID INPUT. C IF (Y(1) .LE. 0.0D0) THEN IRES = -1 ELSE C C CALL F TO OBTAIN F(T,Y) C CALL F1(NEQ,T,Y,DELTA) C C FORM F = F'-F(T,Y) C DO 10 I = 1,NEQ DELTA(I) = YPRIME(I) - DELTA(I) 10 CONTINUE ENDIF C RETURN END C SUBROUTINE F1 (NEQ, T, Y, YDOT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEQ DOUBLE PRECISION T, Y, YDOT DIMENSION Y(1), YDOT(1) YDOT(1) = ((2.0D0*LOG(Y(1)) + 8.0D0)/T - 5.0D0)*Y(1) RETURN END C SUBROUTINE GR1 (NEQ, T, Y, NG, GROOT, RPAR, IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEQ, NG DOUBLE PRECISION T, Y, GROOT DIMENSION Y(1), GROOT(2) GROOT(1) = ((2.0D0*LOG(Y(1)) + 8.0D0)/T - 5.0D0)*Y(1) GROOT(2) = LOG(Y(1)) - 2.2491D0 RETURN END C SUBROUTINE RES2(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /LOCAL/NEQ DIMENSION Y(NEQ), YPRIME(NEQ), DELTA(NEQ) C C CALL F TO OBTAIN F(T,Y) C CALL F2(NEQ,T,Y,DELTA) C C FORM F = F'-F(T,Y) C DO 10 I = 1,NEQ DELTA(I) = YPRIME(I) - DELTA(I) 10 CONTINUE C RETURN END C SUBROUTINE F2 (NEQ, T, Y, YDOT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEQ DOUBLE PRECISION T, Y, YDOT DIMENSION Y(2), YDOT(2) YDOT(1) = Y(2) YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1) RETURN END C SUBROUTINE JAC2 (T, Y, YPRIME, PD, CJ, RPAR, IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEQ, NROWPD DOUBLE PRECISION T, Y, PD PARAMETER (NROWPD=2) DIMENSION Y(2), PD(NROWPD,2) COMMON /LOCAL/NEQ C C FIRST DEFINE THE JACOBIAN MATRIX FOR THE RIGHT HAND SIDE C OF THE ODE: F' = F(T,Y) , I.E. DF/DY) C PD(1,1) = 0.0D0 PD(1,2) = 1.0D0 PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0 PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1)) C C NEXT UPDATE THE JACOBIAN WITH THE RIGHT HAND SIDE TO FORM THE C DAE JACOBIAN: D(F'-F)/DY = DF'/DY - DF/DY = I - DF/DY C PD(1,1) = CJ - PD(1,1) PD(1,2) = - PD(1,2) PD(2,1) = - PD(2,1) PD(2,2) = CJ - PD(2,2) C RETURN END C SUBROUTINE GR2 (NEQ, T, Y, NG, GROOT, RPAR, IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEQ, NG DOUBLE PRECISION T, Y, GROOT DIMENSION Y(2), GROOT(1) GROOT(1) = Y(1) RETURN END