share/odepack/fortran/opkdemo4.f

   1       program opkdemo4
   2 c-----------------------------------------------------------------------
   3 c Demonstration program for the DLSODAR package.
   4 c This is the version of 14 June 2001.
   5 c
   6 c This version is in double precision.
   7 c
   8 c The DLSODAR package is used to solve two simple problems,
   9 c one nonstiff and one intermittently stiff.
  10 c If the errors are too large, or other difficulty occurs,
  11 c a warning message is printed.  All output is on unit lout = 6.
  12 c-----------------------------------------------------------------------
  13       external f1, gr1, f2, jac2, gr2
  14       integer iopt, iout, istate, itask, itol, iwork, jroot, jt,
  15      1   kroot, leniw, lenrw, liw, lrw, lout, neq, nerr, ng,
  16      2   nfe, nfea, nge, nje, nst
  17       double precision atol, er, ero, errt, rtol, rwork,
  18      1   t, tout, tzero, y, yt
  19       dimension y(2), atol(2), rwork(57), iwork(22), jroot(2)
  20       dimension neq(1), rtol(1)
  21       data lout/6/
  22 c
  23       nerr = 0
  24 c-----------------------------------------------------------------------
  25 c First problem.
  26 c The initial value problem is:
  27 c   dy/dt = ((2*log(y) + 8)/t - 5)*y,  y(1) = 1,  1 .le. t .le. 6
  28 c The solution is  y(t) = exp(-t**2 + 5*t - 4)
  29 c The two root functions are:
  30 c   g1 = ((2*log(y)+8)/t - 5)*y (= dy/dt)  (with root at t = 2.5),
  31 c   g2 = log(y) - 2.2491  (with roots at t = 2.47 and 2.53)
  32 c-----------------------------------------------------------------------
  33 c Set all input parameters and print heading.
  34       neq(1) = 1
  35       y(1) = 1.0d0
  36       t = 1.0d0
  37       tout = 2.0d0
  38       itol = 1
  39       rtol(1) = 1.0d-6
  40       atol(1) = 1.0d-6
  41       itask = 1
  42       istate = 1
  43       iopt = 0
  44       lrw = 44
  45       liw = 21
  46       jt = 2
  47       ng = 2
  48       write (lout,110) itol,rtol(1),atol(1),jt
  49  110  format(/' Demonstration program for DLSODAR package'////
  50      1  ' First problem'///
  51      2  ' Problem is  dy/dt = ((2*log(y)+8)/t - 5)*y,  y(1) = 1'//
  52      3  ' Solution is  y(t) = exp(-t**2 + 5*t - 4)'//
  53      4  ' Root functions are:'/
  54      5  10x,' g1 = dy/dt  (root at t = 2.5)'/
  55      6  10x,' g2 = log(y) - 2.2491  (roots at t = 2.47 and t = 2.53)'//
  56      7  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//
  57      8  ' jt =',i3///)
  58 c
  59 c Call DLSODAR in loop over tout values 2,3,4,5,6.
  60       ero = 0.0d0
  61       do 180 iout = 1,5
  62  120    continue
  63         call dlsodar(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
  64      1     iopt,rwork,lrw,iwork,liw,jdum,jt,gr1,ng,jroot)
  65 c
  66 c Print y and error in y, and print warning if error too large.
  67         yt = exp(-t*t + 5.0d0*t - 4.0d0)
  68         er = y(1) - yt
  69         write (lout,130) t,y(1),er
  70  130    format(' At t =',d15.7,5x,'y =',d15.7,5x,'error =',d12.4)
  71         if (istate .lt. 0) go to 185
  72         er = abs(er)/(rtol(1)*abs(y(1)) + atol(1))
  73         ero = max(ero,er)
  74         if (er .gt. 1000.0d0) then
  75           write (lout,140)
  76  140      format(//' Warning: error exceeds 1000 * tolerance'//)
  77           nerr = nerr + 1
  78           endif
  79         if (istate .ne. 3) go to 175
  80 c
  81 c If a root was found, write results and check root location.
  82 c Then reset istate to 2 and return to DLSODAR call.
  83         write (lout,150) t,jroot(1),jroot(2)
  84  150    format(/' Root found at t =',d15.7,5x,'jroot =',2i5)
  85         if (jroot(1) .eq. 1) errt = t - 2.5d0
  86         if (jroot(2) .eq. 1 .and. t .le. 2.5d0) errt = t - 2.47d0
  87         if (jroot(2) .eq. 1 .and. t .gt. 2.5d0) errt = t - 2.53d0
  88         write (lout,160) errt
  89  160    format(' Error in t location of root is',d12.4/)
  90         if (abs(errt) .gt. 1.0d-3) then
  91           write (lout,170)
  92  170      format(//' Warning: root error exceeds 1.0d-3'//)
  93           nerr = nerr + 1
  94           endif
  95         istate = 2
  96         go to 120
  97 c
  98 c If no root found, increment tout and loop back.
  99  175    tout = tout + 1.0d0
 100  180    continue
 101 c
 102 c Problem complete.  Print final statistics.
 103  185  continue
 104       if (istate .lt. 0) nerr = nerr + 1
 105       nst = iwork(11)
 106       nfe = iwork(12)
 107       nje = iwork(13)
 108       nge = iwork(10)
 109       lenrw = iwork(17)
 110       leniw = iwork(18)
 111       nfea = nfe
 112       if (jt .eq. 2) nfea = nfe - neq(1)*nje
 113       write (lout,190) lenrw,leniw,nst,nfe,nfea,nje,nge,ero
 114  190  format(//' Final statistics for this run:'/
 115      1  ' rwork size =',i4,'   iwork size =',i4/
 116      2  ' number of steps =',i5/
 117      3  ' number of f-s   =',i5/
 118      4  ' (excluding j-s) =',i5/
 119      5  ' number of j-s   =',i5/
 120      6  ' number of g-s   =',i5/
 121      7  ' error overrun =',d10.2)
 122 c
 123 c-----------------------------------------------------------------------
 124 c Second problem (Van der Pol oscillator).
 125 c The initial value problem is (after reduction of 2nd order ODE):
 126 c   dy1/dt = y2,  dy2/dt = 100*(1 - y1**2)*y2 - y1,
 127 c   y1(0) = 2,  y2(0) = 0,  0 .le. t .le. 200
 128 c The root function is  g = y1.
 129 c An analytic solution is not known, but the zeros of y1 are known
 130 c to 15 figures for purposes of checking the accuracy.
 131 c-----------------------------------------------------------------------
 132 c Set tolerance parameters and print heading.
 133       itol = 2
 134       rtol(1) = 1.0d-6
 135       atol(1) = 1.0d-6
 136       atol(2) = 1.0d-4
 137       write (lout,200) itol,rtol(1),atol(1),atol(2)
 138  200  format(////80('*')//' Second problem (Van der Pol oscillator)'//
 139      1  ' Problem is dy1/dt = y2,  dy2/dt = 100*(1-y1**2)*y2 - y1'/
 140      2  '            y1(0) = 2,  y2(0) = 0'//
 141      3  ' Root function is  g = y1'//
 142      4  ' itol =',i3,'   rtol =',d10.1,'   atol =',2d10.1)
 143 c
 144 c Loop over jt = 1, 2.  Set remaining parameters and print jt.
 145       do 290 jt = 1,2
 146       neq(1) = 2
 147       y(1) = 2.0d0
 148       y(2) = 0.0d0
 149       t = 0.0d0
 150       tout = 20.0d0
 151       itask = 1
 152       istate = 1
 153       iopt = 0
 154       lrw = 57
 155       liw = 22
 156       ng = 1
 157       write (lout,210) jt
 158  210  format(///' Solution with jt =',i2//)
 159 c
 160 c Call DLSODAR in loop over tout values 20,40,...,200.
 161       do 270 iout = 1,10
 162  220    continue
 163         call dlsodar(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
 164      1     iopt,rwork,lrw,iwork,liw,jac2,jt,gr2,ng,jroot)
 165 c
 166 c Print y1 and y2.
 167         write (lout,230) t,y(1),y(2)
 168  230    format(' At t =',d15.7,5x,'y1 =',d15.7,5x,'y2 =',d15.7)
 169         if (istate .lt. 0) go to 275
 170         if (istate .ne. 3) go to 265
 171 c
 172 c If a root was found, write results and check root location.
 173 c Then reset istate to 2 and return to DLSODAR call.
 174         write (lout,240) t
 175  240    format(/' Root found at t =',d15.7)
 176         kroot = int(t/81.2d0 + 0.5d0)
 177         tzero = 81.17237787055d0 + (kroot-1)*81.41853556212d0
 178         errt = t - tzero
 179         write (lout,250) errt
 180  250    format(' Error in t location of root is',d12.4//)
 181         if (abs(errt) .gt. 1.0d-1) then
 182           write (lout,260)
 183  260      format(//' Warning: root error exceeds 1.0d-1'//)
 184           nerr = nerr + 1
 185           endif
 186         istate = 2
 187         go to 220
 188 c
 189 c If no root found, increment tout and loop back.
 190  265    tout = tout + 20.0d0
 191  270    continue
 192 c
 193 c Problem complete.  Print final statistics.
 194  275  continue
 195       if (istate .lt. 0) nerr = nerr + 1
 196       nst = iwork(11)
 197       nfe = iwork(12)
 198       nje = iwork(13)
 199       nge = iwork(10)
 200       lenrw = iwork(17)
 201       leniw = iwork(18)
 202       nfea = nfe
 203       if (jt .eq. 2) nfea = nfe - neq(1)*nje
 204       write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,nge
 205  280  format(//' Final statistics for this run:'/
 206      1  '  rwork size =',i4,'   iwork size =',i4/
 207      2  '  number of steps =',i5/
 208      3  '  number of f-s   =',i5/
 209      4  '  (excluding j-s) =',i5/
 210      5  '  number of j-s   =',i5/
 211      6  '  number of g-s   =',i5)
 212  290  continue
 213 c
 214       write (lout,300) nerr
 215  300  format(///' Total number of errors encountered =',i3)
 216 c      stop
 217       end
 218
 219       subroutine f1 (neq, t, y, ydot)
 220       integer neq
 221       double precision t, y, ydot
 222       dimension y(1), ydot(1)
 223       dimension neq(*)
 224       ydot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1)
 225       return
 226       end
 227
 228       subroutine gr1 (neq, t, y, ng, groot)
 229       integer neq, ng
 230       double precision t, y, groot
 231       dimension y(1), groot(2)
 232       dimension neq(*)
 233       groot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1)
 234       groot(2) = log(y(1)) - 2.2491d0
 235       return
 236       end
 237
 238       subroutine f2 (neq, t, y, ydot)
 239       integer neq
 240       double precision t, y, ydot
 241       dimension y(2), ydot(2)
 242       dimension neq(*)
 243       ydot(1) = y(2)
 244       ydot(2) = 100.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
 245       return
 246       end
 247
 248       subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
 249       integer neq, ml, mu, nrowpd
 250       double precision t, y, pd
 251       dimension y(2), pd(nrowpd,2)
 252       dimension neq(*)
 253       pd(1,1) = 0.0d0
 254       pd(1,2) = 1.0d0
 255       pd(2,1) = -200.0d0*y(1)*y(2) - 1.0d0
 256       pd(2,2) = 100.0d0*(1.0d0 - y(1)*y(1))
 257       return
 258       end
 259
 260       subroutine gr2 (neq, t, y, ng, groot)
 261       integer neq, ng
 262       double precision t, y, groot
 263       dimension y(2), groot(1)
 264       dimension neq(*)
 265       groot(1) = y(1)
 266       return
 267       end