share/odepack/fortran/opkdemo1.f

   1       program opkdemo1
   2 c-----------------------------------------------------------------------
   3 c Demonstration program for the DLSODE 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 package is used to solve two simple problems,
   9 c one with a full Jacobian, the other with a banded Jacobian,
  10 c with all 8 of the appropriate values of mf in each case.
  11 c If the errors are too large, or other difficulty occurs,
  12 c a warning message is printed.  All output is on unit lout = 6.
  13 c-----------------------------------------------------------------------
  14       external f1, jac1, f2, jac2
  15       integer i, iopar, iopt, iout, istate, itask, itol, iwork,
  16      1   leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter,
  17      2   ml, mu, neq(1), nerr, nfe, nfea, nje, nout, nqu, nst
  18       double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t,
  19      1   tout, tout1, y
  20       dimension y(25), rwork(697), iwork(45), rtol(1), atol(1)
  21       data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/
  22 c
  23       nerr = 0
  24       itol = 1
  25       rtol(1) = 0.0d0
  26       atol(1) = 1.0d-6
  27       lrw = 697
  28       liw = 45
  29       iopt = 0
  30 c
  31 c First problem
  32 c
  33       neq(1) = 2
  34       nout = 4
  35       write (lout,110) neq(1),itol,rtol(1),atol(1)
  36  110  format(/' Demonstration program for DLSODE package'///
  37      1  ' Problem 1:  Van der Pol oscillator:'/
  38      2  '  xdotdot - 3*(1 - x**2)*xdot + x = 0, ',
  39      3  '   x(0) = 2, xdot(0) = 0'/
  40      4  ' neq =',i2/
  41      5  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
  42 c
  43       do 195 meth = 1,2
  44       do 190 miter = 0,3
  45       mf = 10*meth + miter
  46       write (lout,120) mf
  47  120  format(///' Solution with mf =',i3//
  48      1     5x,'t               x               xdot       nq      h'//)
  49       t = 0.0d0
  50       y(1) = 2.0d0
  51       y(2) = 0.0d0
  52       itask = 1
  53       istate = 1
  54       tout = tout1
  55       ero = 0.0d0
  56       do 170 iout = 1,nout
  57         call dlsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
  58      1     iopt,rwork,lrw,iwork,liw,jac1,mf)
  59         hu = rwork(11)
  60         nqu = iwork(14)
  61         write (lout,140) t,y(1),y(2),nqu,hu
  62  140    format(d15.5,d16.5,d14.3,i5,d14.3)
  63         if (istate .lt. 0) go to 175
  64         iopar = iout - 2*(iout/2)
  65         if (iopar .ne. 0) go to 170
  66         er = abs(y(1))/atol(1)
  67         ero = max(ero,er)
  68         if (er .gt. 1000.0d0) then
  69           write (lout,150)
  70  150      format(//' Warning: error exceeds 1000 * tolerance'//)
  71           nerr = nerr + 1
  72         endif
  73  170    tout = tout + dtout
  74  175  continue
  75       if (istate .lt. 0) nerr = nerr + 1
  76       nst = iwork(11)
  77       nfe = iwork(12)
  78       nje = iwork(13)
  79       lenrw = iwork(17)
  80       leniw = iwork(18)
  81       nfea = nfe
  82       if (miter .eq. 2) nfea = nfe - neq(1)*nje
  83       if (miter .eq. 3) nfea = nfe - nje
  84       write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
  85  180  format(//' Final statistics for this run:'/
  86      1  ' rwork size =',i4,'   iwork size =',i4/
  87      2  ' number of steps =',i5/
  88      3  ' number of f-s   =',i5/
  89      4  ' (excluding J-s) =',i5/
  90      5  ' number of J-s   =',i5/
  91      6  ' error overrun =',d10.2)
  92  190  continue
  93  195  continue
  94 c
  95 c Second problem
  96 c
  97       neq(1) = 25
  98       ml = 5
  99       mu = 0
 100       iwork(1) = ml
 101       iwork(2) = mu
 102       mband = ml + mu + 1
 103       nout = 5
 104       write (lout,210) neq(1),ml,mu,itol,rtol(1),atol(1)
 105  210  format(///70('-')///
 106      1  ' Problem 2: ydot = A * y , where',
 107      2  '  A is a banded lower triangular matrix'/
 108      3     12x, 'derived from 2-D advection PDE'/
 109      4  ' neq =',i3,'   ml =',i2,'   mu =',i2/
 110      5  ' itol =',i3,'   rtol =',d10.1,'   atol =',d10.1//)
 111       do 295 meth = 1,2
 112       do 290 miter = 0,5
 113       if (miter .eq. 1 .or. miter .eq. 2) go to 290
 114       mf = 10*meth + miter
 115       write (lout,220) mf
 116  220  format(///' Solution with mf =',i3//
 117      1       5x,'t             max.err.     nq      h'//)
 118       t = 0.0d0
 119       do 230 i = 2,neq(1)
 120  230    y(i) = 0.0d0
 121       y(1) = 1.0d0
 122       itask = 1
 123       istate = 1
 124       tout = 0.01d0
 125       ero = 0.0d0
 126       do 270 iout = 1,nout
 127         call dlsode(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
 128      1     iopt,rwork,lrw,iwork,liw,jac2,mf)
 129         call edit2(y,t,erm)
 130         hu = rwork(11)
 131         nqu = iwork(14)
 132         write (lout,240) t,erm,nqu,hu
 133  240    format(d15.5,d14.3,i5,d14.3)
 134         if (istate .lt. 0) go to 275
 135         er = erm/atol(1)
 136         ero = max(ero,er)
 137         if (er .gt. 1000.0d0) then
 138           write (lout,150)
 139           nerr = nerr + 1
 140         endif
 141  270    tout = tout*10.0d0
 142  275  continue
 143       if (istate .lt. 0) nerr = nerr + 1
 144       nst = iwork(11)
 145       nfe = iwork(12)
 146       nje = iwork(13)
 147       lenrw = iwork(17)
 148       leniw = iwork(18)
 149       nfea = nfe
 150       if (miter .eq. 5) nfea = nfe - mband*nje
 151       if (miter .eq. 3) nfea = nfe - nje
 152       write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
 153  290  continue
 154  295  continue
 155       write (lout,300) nerr
 156  300  format(////' Number of errors encountered =',i3)
 157 c      stop
 158       end
 159
 160       subroutine f1 (neq, t, y, ydot)
 161       integer neq
 162       double precision t, y, ydot
 163       dimension y(*), ydot(*)
 164       dimension neq(*)
 165       ydot(1) = y(2)
 166       ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
 167       return
 168       end
 169
 170       subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
 171       integer neq, ml, mu, nrowpd
 172       double precision t, y, pd
 173       dimension y(*), pd(nrowpd,*)
 174       dimension neq(*)
 175       pd(1,1) = 0.0d0
 176       pd(1,2) = 1.0d0
 177       pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0
 178       pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1))
 179       return
 180       end
 181
 182       subroutine f2 (neq, t, y, ydot)
 183       integer neq, i, j, k, ng
 184       double precision t, y, ydot, alph1, alph2, d
 185       dimension y(*), ydot(*)
 186       dimension neq(*)
 187       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 188       do 10 j = 1,ng
 189       do 10 i = 1,ng
 190         k = i + (j - 1)*ng
 191         d = -2.0d0*y(k)
 192         if (i .ne. 1) d = d + y(k-1)*alph1
 193         if (j .ne. 1) d = d + y(k-ng)*alph2
 194  10     ydot(k) = d
 195       return
 196       end
 197
 198       subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
 199       integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
 200       double precision t, y, pd, alph1, alph2
 201       dimension y(*), pd(nrowpd,*)
 202       dimension neq(*)
 203       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 204       mband = ml + mu + 1
 205       mu1 = mu + 1
 206       mu2 = mu + 2
 207       do 10 j = 1,neq(1)
 208         pd(mu1,j) = -2.0d0
 209         pd(mu2,j) = alph1
 210  10     pd(mband,j) = alph2
 211       do 20 j = ng,neq(1),ng
 212  20     pd(mu2,j) = 0.0d0
 213       return
 214       end
 215
 216       subroutine edit2 (y, t, erm)
 217       integer i, j, k, ng
 218       double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt
 219       dimension y(25)
 220       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
 221       erm = 0.0d0
 222       if (t .eq. 0.0d0) return
 223       ex = 0.0d0
 224       if (t .le. 30.0d0) ex = exp(-2.0d0*t)
 225       a2 = 1.0d0
 226       do 60 j = 1,ng
 227         a1 = 1.0d0
 228         do 50 i = 1,ng
 229           k = i + (j - 1)*ng
 230           yt = t**(i+j-2)*ex*a1*a2
 231           er = abs(y(k)-yt)
 232           erm = max(erm,er)
 233           a1 = a1*alph1/i
 234  50       continue
 235         a2 = a2*alph2/j
 236  60     continue
 237       return
 238       end