[2696] | 1 | MODULE EXAMPLE2 |
---|
| 2 | |
---|
| 3 | ! DEMONSTRATION PROGRAM FOR THE DVODE_F90 PACKAGE. |
---|
| 4 | |
---|
| 5 | ! The following is a modification of the previous example |
---|
| 6 | ! program to illustrate root finding. The problem is from |
---|
| 7 | ! chemical kinetics, and consists of the following three |
---|
| 8 | ! rate equations: |
---|
| 9 | ! dy1/dt = -.04d0*y1 + 1.d4*y2*y3 |
---|
| 10 | ! dy2/dt = .04d0*y1 - 1.d4*y2*y3 - 3.d7*y2**2 |
---|
| 11 | ! dy3/dt = 3.d7*y2**2 |
---|
| 12 | ! on the interval from t = 0.0d0 to t = 4.d10, with initial |
---|
| 13 | ! conditions y1 = 1.0d0, y2 = y3 = 0.0d0. The problem is stiff. |
---|
| 14 | ! In addition, we want to find the values of t, y1, y2, |
---|
| 15 | ! and y3 at which: |
---|
| 16 | ! (1) y1 reaches the value 1.d-4, and |
---|
| 17 | ! (2) y3 reaches the value 1.d-2. |
---|
| 18 | |
---|
| 19 | ! The following coding solves this problem with DVODE_F90 |
---|
| 20 | ! using an internally generated dense Jacobian and |
---|
| 21 | ! printing results at t = .4, 4., ..., 4.d10, and at the |
---|
| 22 | ! computed roots. It uses ITOL = 2 and ATOL much smaller |
---|
| 23 | ! for y2 than y1 or y3 because y2 has much smaller values. |
---|
| 24 | ! At the end of the run, statistical quantities of interest |
---|
| 25 | ! are printed (see optional outputs in the full description |
---|
| 26 | ! below). Output is written to the file example2.dat. |
---|
| 27 | |
---|
| 28 | CONTAINS |
---|
| 29 | |
---|
| 30 | SUBROUTINE FEX(NEQ,T,Y,YDOT) |
---|
| 31 | IMPLICIT NONE |
---|
| 32 | INTEGER NEQ |
---|
| 33 | DOUBLE PRECISION T, Y, YDOT |
---|
| 34 | DIMENSION Y(3), YDOT(3) |
---|
| 35 | |
---|
| 36 | YDOT(1) = -0.04D0*Y(1) + 1.0D4*Y(2)*Y(3) |
---|
| 37 | YDOT(3) = 3.0D7*Y(2)*Y(2) |
---|
| 38 | YDOT(2) = -YDOT(1) - YDOT(3) |
---|
| 39 | RETURN |
---|
| 40 | END SUBROUTINE FEX |
---|
| 41 | |
---|
| 42 | SUBROUTINE GEX(NEQ,T,Y,NG,GOUT) |
---|
| 43 | IMPLICIT NONE |
---|
| 44 | INTEGER NEQ, NG |
---|
| 45 | DOUBLE PRECISION T, Y, GOUT |
---|
| 46 | DIMENSION Y(3), GOUT(2) |
---|
| 47 | |
---|
| 48 | GOUT(1) = Y(1) - 1.0D-4 |
---|
| 49 | GOUT(2) = Y(3) - 1.0D-2 |
---|
| 50 | RETURN |
---|
| 51 | END SUBROUTINE GEX |
---|
| 52 | |
---|
| 53 | END MODULE EXAMPLE2 |
---|
| 54 | |
---|
| 55 | !****************************************************************** |
---|
| 56 | |
---|
| 57 | PROGRAM RUNEXAMPLE2 |
---|
| 58 | |
---|
| 59 | USE DVODE_F90_M |
---|
| 60 | USE EXAMPLE2 |
---|
| 61 | |
---|
| 62 | IMPLICIT NONE |
---|
| 63 | INTEGER ITASK, ISTATE, NG, NEQ, IOUT, JROOT, ISTATS, IERROR, I |
---|
| 64 | DOUBLE PRECISION ATOL, RTOL, RSTATS, T, TOUT, Y |
---|
| 65 | DIMENSION Y(3), ATOL(3), RSTATS(22), ISTATS(31), JROOT(2) |
---|
| 66 | |
---|
| 67 | TYPE (VODE_OPTS) :: OPTIONS |
---|
| 68 | |
---|
| 69 | OPEN (UNIT=6,FILE='example2.dat') |
---|
| 70 | IERROR = 0 |
---|
| 71 | NEQ = 3 |
---|
| 72 | Y(1) = 1.0D0 |
---|
| 73 | Y(2) = 0.0D0 |
---|
| 74 | Y(3) = 0.0D0 |
---|
| 75 | T = 0.0D0 |
---|
| 76 | TOUT = 0.4D0 |
---|
| 77 | RTOL = 1.0D-4 |
---|
| 78 | ATOL(1) = 1.0D-8 |
---|
| 79 | ATOL(2) = 1.0D-10 |
---|
| 80 | ATOL(3) = 1.0D-8 |
---|
| 81 | ITASK = 1 |
---|
| 82 | ISTATE = 1 |
---|
| 83 | NG = 2 |
---|
| 84 | ! OPTIONS = SET_OPTS(DENSE_J=.TRUE.,RELERR=RTOL,ABSERR_VECTOR=ATOL, & |
---|
| 85 | ! NEVENTS=NG) |
---|
| 86 | OPTIONS = SET_NORMAL_OPTS(DENSE_J=.TRUE.,RELERR=RTOL, & |
---|
| 87 | ABSERR_VECTOR=ATOL,NEVENTS=NG) |
---|
| 88 | DO 20 IOUT = 1, 12 |
---|
| 89 | 10 CONTINUE |
---|
| 90 | CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=GEX) |
---|
| 91 | CALL GET_STATS(RSTATS,ISTATS,NG,JROOT) |
---|
| 92 | WRITE (6,90000) T, Y(1), Y(2), Y(3) |
---|
| 93 | DO I = 1, NEQ |
---|
| 94 | IF (Y(I)<0.0D0) IERROR = 1 |
---|
| 95 | END DO |
---|
| 96 | 90000 FORMAT (' At t =',D12.4,' Y =',3D14.6) |
---|
| 97 | IF (ISTATE<0) GO TO 30 |
---|
| 98 | IF (ISTATE==2) GO TO 20 |
---|
| 99 | WRITE (6,90001) JROOT(1), JROOT(2) |
---|
| 100 | 90001 FORMAT (5X,' The above line is a root, JROOT =',2I5) |
---|
| 101 | ISTATE = 2 |
---|
| 102 | GO TO 10 |
---|
| 103 | 20 TOUT = TOUT*10.0D0 |
---|
| 104 | WRITE (6,90002) ISTATS(11), ISTATS(12), ISTATS(13), ISTATS(10) |
---|
| 105 | IF (IERROR==1) THEN |
---|
| 106 | WRITE (6,90003) |
---|
| 107 | ELSE |
---|
| 108 | WRITE (6,90004) |
---|
| 109 | END IF |
---|
| 110 | 90002 FORMAT (/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4, & |
---|
| 111 | ' No. g-s =',I4/) |
---|
| 112 | STOP |
---|
| 113 | 30 WRITE (6,90005) ISTATE |
---|
| 114 | 90003 FORMAT (/' An error occurred.') |
---|
| 115 | 90004 FORMAT (/' No errors occurred.') |
---|
| 116 | 90005 FORMAT (///' Error halt.. ISTATE =',I3) |
---|
| 117 | STOP |
---|
| 118 | END PROGRAM RUNEXAMPLE2 |
---|