share/odepack/fortran/opkdemo7.f

   1       program opkdemo7
   2 c-----------------------------------------------------------------------
   3 c Demonstration program for the DLSODI 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 this program solves a semi-discretized form of the Burgers equation,
   9 c
  10 c     u  = -(u*u/2)  + eta * u
  11 c      t           x          xx
  12 c
  13 c for a = -1 .le. x .le. 1 = b, t .ge. 0.
  14 c Here eta = 0.05.
  15 c Boundary conditions: u(-1,t) = u(1,t) = 0.
  16 c Initial profile: square wave
  17 c     u(0,x) = 0    for 1/2 .lt. abs(x) .le. 1
  18 c     u(0,x) = 1/2  for abs(x) = 1/2
  19 c     u(0,x) = 1    for 0 .le. abs(x) .lt. 1/2
  20 c
  21 c An ODE system is generated by a simplified Galerkin treatment
  22 c of the spatial variable x.
  23 c
  24 c Reference:
  25 c R. C. Y. Chin, G. W. Hedstrom, and K. E. Karlsson,
  26 c A Simplified Galerkin Method for Hyperbolic Equations,
  27 c Math. Comp., vol. 33, no. 146 (April 1979), pp. 647-658.
  28 c
  29 c The problem is run with the DLSODI package with a 10-point mesh
  30 c and a 100-point mesh.  In each case, it is run with two tolerances
  31 c and for various appropriate values of the method flag mf.
  32 c Output is on unit lout, set to 6 in a data statement below.
  33 c-----------------------------------------------------------------------
  34       external res, addabd, addafl, jacbd, jacfl
  35       integer i, io, istate, itol, iwork, j,
  36      1   lout, liw, lrw, meth, miter, mf, ml, mu,
  37      2   n, nout, npts, nerr,
  38      3   nptsm1, n14, n34, n14m1, n14p1, n34m1, n34p1
  39       integer nm1
  40       double precision a, b, eta, delta,
  41      1   zero, fourth, half, one, hun,
  42      2   t, tout, tlast, tinit, errfac,
  43      3   atol, rtol, rwork, y, ydoti, elkup
  44       double precision eodsq, r4d
  45       dimension n(1)
  46       dimension y(99), ydoti(99), tout(4), atol(2), rtol(2)
  47       dimension rwork(2002), iwork(125)
  48 c Pass problem parameters in the Common block test1.
  49       common /test1/ r4d, eodsq, nm1
  50 c
  51 c Set problem parameters and run parameters
  52       data eta/0.05d0/, a/-1.0d0/, b/1.0d0/
  53       data zero/0.0d0/, fourth/0.25d0/, half/.5d0/, one/1.0d0/,
  54      1   hun/100.0d0/
  55       data tinit/0.0d0/, tlast/0.4d0/
  56       data tout/.10d0,.20d0,.30d0,.40d0/
  57       data ml/1/, mu/1/, lout/6/
  58       data nout/4/, lrw/2002/, liw/125/
  59       data itol/1/, rtol/1.0d-3, 1.0d-6/, atol/1.0d-3, 1.0d-6/
  60 c
  61       iwork(1) = ml
  62       iwork(2) = mu
  63       nerr = 0
  64 c
  65 c Loop over two values of npts.
  66       do 300  npts = 10, 100, 90
  67 c
  68 c Compute the mesh width delta and other parameters.
  69       delta = (b - a)/npts
  70       r4d = fourth/delta
  71       eodsq = eta/delta**2
  72       nptsm1 = npts - 1
  73       n14 = npts/4
  74       n34 = 3 * n14
  75       n14m1 = n14 - 1
  76       n14p1 = n14m1 + 2
  77       n34m1 = n34 - 1
  78       n34p1 = n34m1 + 2
  79       n(1) = nptsm1
  80       nm1 = n(1) - 1
  81 c
  82 c Set the initial profile (for output purposes only).
  83 c
  84       do 10 i = 1,n14m1
  85    10   y(i) = zero
  86       y(n14) = half
  87       do 20 i = n14p1,n34m1
  88    20   y(i) = one
  89       y(n34) = half
  90       do 30 i = n34p1,nptsm1
  91    30   y(i) = zero
  92 c
  93       if (npts .gt. 10) write (lout,1010)
  94       write (lout,1000)
  95       write (lout,1100) eta,a,b,tinit,tlast,ml,mu,n(1)
  96       write (lout,1200) zero, (y(i), i=1,n(1)), zero
  97 c
  98 c The j loop is over error tolerances.
  99 c
 100       do 200 j = 1,2
 101 c
 102 c Loop over method flag loop (for demonstration).
 103 c
 104       do 100 meth = 1,2
 105        do 100 miter = 1,5
 106         if (miter .eq. 3)  go to 100
 107         if (miter .le. 2 .and. npts .gt. 10)  go to 100
 108         if (miter .eq. 5 .and. npts .lt. 100)  go to 100
 109         mf = 10*meth + miter
 110 c
 111 c Set the initial profile.
 112 c
 113         do 40 i = 1,n14m1
 114    40     y(i) = zero
 115         y(n14) = half
 116         do 50 i = n14p1,n34m1
 117    50     y(i) = one
 118         y(n34) = half
 119         do 60 i = n34p1,nptsm1
 120    60     y(i) = zero
 121 c
 122         t = tinit
 123         istate = 0
 124 c
 125         write (lout,1500) rtol(j), atol(j), mf, npts
 126 c
 127 c  Output loop for each case
 128 c
 129         do 80 io = 1,nout
 130 c
 131 c         call DLSODI
 132           if (miter .le. 2) call dlsodi (res, addafl, jacfl, n, y,
 133      1                      ydoti, t, tout(io), itol, rtol(j), atol(j),
 134      2                      1, istate, 0, rwork, lrw, iwork, liw, mf)
 135           if (miter .ge. 4) call dlsodi (res, addabd, jacbd, n, y,
 136      1                      ydoti, t, tout(io), itol, rtol(j), atol(j),
 137      2                      1, istate, 0, rwork, lrw, iwork, liw, mf)
 138           write (lout,2000) t, rwork(11), iwork(14),(y(i), i=1,n(1))
 139 c
 140 c If istate is not 2 on return, print message and loop.
 141           if (istate .ne. 2) then
 142             write (lout,4000) mf, t, istate
 143             nerr = nerr + 1
 144             go to 100
 145           endif
 146 c
 147    80     continue
 148 c
 149         write (lout,3000) mf, iwork(11), iwork(12), iwork(13),
 150      1                iwork(17), iwork(18)
 151 c
 152 c Estimate final error and print result.
 153         itemp = n(1)
 154         errfac = elkup( itemp, y, rwork(21), itol, rtol(j), atol(j) )
 155         if (errfac .gt. hun) then
 156           write (lout,5001)  errfac
 157           nerr = nerr + 1
 158         else
 159           write (lout,5000)  errfac
 160         endif
 161   100   continue
 162   200 continue
 163   300 continue
 164 c
 165       write (lout,6000) nerr
 166 c      stop
 167 c
 168  1000 format(20x,' Demonstration Problem for DLSODI')
 169  1010 format(///80('*')///)
 170  1100 format(/10x,' Simplified Galerkin Solution of Burgers Equation'//
 171      1       13x,'Diffusion coefficient is eta =',d10.2/
 172      2       13x,'Uniform mesh on interval',d12.3,' to ',d12.3/
 173      3       13x,'Zero boundary conditions'/
 174      4       13x,'Time limits: t0 = ',d12.5,'   tlast = ',d12.5/
 175      5       13x,'Half-bandwidths ml = ',i2,'   mu = ',i2/
 176      6       13x,'System size neq = ',i3/)
 177 c
 178  1200 format('Initial profile:'/17(6d12.4/))
 179 c
 180  1500 format(///80('-')///'Run with rtol =',d12.2,'  atol =',d12.2,
 181      1       '   mf =',i3,'   npts =',i4,':'//)
 182 c
 183  2000 format('Output for time t = ',d12.5,'   current h =',
 184      1       d12.5,'   current order =',i2,':'/17(6d12.4/))
 185 c
 186  3000 format(//'Final statistics for mf = ',i2,':'/
 187      1       i4,' steps,',i5,' res,',i4,' Jacobians,',
 188      2       '   rwork size =',i6,',   iwork size =',i6)
 189 c
 190  4000 format(///80('*')//20x,'Final time reached for mf = ',i2,
 191      1       ' was t = ',d12.5/25x,'at which istate = ',i2////80('*'))
 192  5000 format('  Final output is correct to within ',d8.1,
 193      1       '  times local error tolerance')
 194  5001 format('  Final output is wrong by ',d8.1,
 195      1       '  times local error tolerance')
 196  6000 format(//80('*')//
 197      1       'Run completed.  Number of errors encountered =',i3)
 198 c
 199 c end of main program for the DLSODI demonstration problem.
 200       end
 201
 202       subroutine gfun (n, t, y, g)
 203 c This subroutine computes the right-hand side function g(y,t).
 204 c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1
 205 c from the Common block test1.
 206 c
 207       integer i, n, nm1
 208       double precision t, y, g, r4d, eodsq, two
 209       dimension g(n), y(n)
 210       common /test1/ r4d, eodsq, nm1
 211       data two/2.0d0/
 212 c
 213       g(1) = -r4d*y(2)**2 + eodsq*(y(2) - two*y(1))
 214 c
 215       do 20 i = 2,nm1
 216         g(i) = r4d*(y(i-1)**2 - y(i+1)**2)
 217      1        + eodsq*(y(i+1) - two*y(i) + y(i-1))
 218    20   continue
 219 c
 220       g(n) = r4d*y(nm1)**2 + eodsq*(y(nm1) - two*y(n))
 221 c
 222       return
 223 c end of subroutine gfun for the DLSODI demonstration problem.
 224       end
 225
 226       subroutine res (n, t, y, v, r, ires)
 227 c This subroutine computes the residual vector
 228 c   r = g(t,y) - A(t,y)*v .
 229 c It uses nm1 = n - 1 from Common.
 230 c If ires = -1, only g(t,y) is returned in r, since A(t,y) does
 231 c not depend on y.
 232 c
 233       integer i, ires, n, nm1
 234       double precision t, y, v, r, r4d, eodsq, one, four, six,
 235      1   fact1, fact4
 236       dimension y(*), v(*), r(*)
 237       dimension n(*)
 238       integer itemp
 239       common /test1/ r4d, eodsq, nm1
 240       data one /1.0d0/, four /4.0d0/, six /6.0d0/
 241 c
 242       itemp = n(1)
 243       call gfun (itemp, t, y, r)
 244       if (ires .eq. -1) return
 245 c
 246       fact1 = one/six
 247       fact4 = four/six
 248       r(1) = r(1) - (fact4*v(1) + fact1*v(2))
 249       do 10 i = 2, nm1
 250   10   r(i) = r(i) - (fact1*v(i-1) + fact4*v(i) + fact1*v(i+1))
 251       r(n(1)) = r(n(1)) - (fact1*v(nm1) + fact4*v(n(1)))
 252       return
 253 c end of subroutine res for the DLSODI demonstration problem.
 254       end
 255
 256       subroutine addabd (n, t, y, ml, mu, pa, m0)
 257 c This subroutine computes the matrix A in band form, adds it to pa,
 258 c and returns the sum in pa.   The matrix A is tridiagonal, of order n,
 259 c with nonzero elements (reading across) of  1/6, 4/6, 1/6.
 260 c
 261       integer i, n, m0, ml, mu, mup1, mup2
 262       double precision t, y, pa, fact1, fact4, one, four, six
 263       dimension y(*), pa(m0,*)
 264       dimension n(*)
 265       data one/1.0d0/, four/4.0d0/, six/6.0d0/
 266 c
 267 c Set the pointers.
 268       mup1 = mu + 1
 269       mup2 = mu + 2
 270 c Compute the elements of A.
 271       fact1 = one/six
 272       fact4 = four/six
 273 c Add the matrix A to the matrix pa (banded).
 274       do 10 i = 1,n(1)
 275         pa(mu,i) = pa(mu,i) + fact1
 276         pa(mup1,i) = pa(mup1,i) + fact4
 277         pa(mup2,i) = pa(mup2,i) + fact1
 278    10   continue
 279       return
 280 c end of subroutine addabd for the DLSODI demonstration problem.
 281       end
 282
 283       subroutine addafl (n, t, y, ml, mu, pa, m0)
 284 c This subroutine computes the matrix A in full form, adds it to
 285 c pa, and returns the sum in pa.
 286 c It uses nm1 = n - 1 from Common.
 287 c The matrix A is tridiagonal, of order n, with nonzero elements
 288 c (reading across) of  1/6, 4/6, 1/6.
 289 c
 290       integer i, n, m0, ml, mu, nm1
 291       double precision t, y, pa, r4d, eodsq, one, four, six,
 292      1   fact1, fact4
 293       dimension y(*), pa(m0,8)
 294       dimension n(1)
 295       common /test1/ r4d, eodsq, nm1
 296       data one/1.0d0/, four/4.0d0/, six/6.0d0/
 297 c
 298 c Compute the elements of A.
 299       fact1 = one/six
 300       fact4 = four/six
 301 c
 302 c Add the matrix A to the matrix pa (full).
 303 c
 304       do 110  i = 2, nm1
 305          pa(i,i+1) = pa(i,i+1) + fact1
 306          pa(i,i) = pa(i,i) + fact4
 307          pa(i,i-1) = pa(i,i-1) + fact1
 308   110    continue
 309       pa(1,2) = pa(1,2) + fact1
 310       pa(1,1) = pa(1,1) + fact4
 311       pa(n(1),n(1)) = pa(n(1),n(1)) + fact4
 312       pa(n(1),nm1) = pa(n(1),nm1) + fact1
 313       return
 314 c end of subroutine addafl for the DLSODI demonstration problem.
 315       end
 316
 317       subroutine jacbd (n, t, y, s, ml, mu, pa, m0)
 318 c This subroutine computes the Jacobian dg/dy = d(g-a*s)/dy
 319 c and stores elements
 320 c   i   j
 321 c dg /dy   in  pa(i-j+mu+1,j)  in band matrix format.
 322 c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1
 323 c from the Common block test1.
 324 c
 325       integer i, n, m0, ml, mu, mup1, mup2, nm1
 326       double precision t, y, s, pa, diag, r4d, eodsq, two, r2d
 327       dimension y(*), s(*), pa(m0,*)
 328       dimension n(*)
 329       common /test1/ r4d, eodsq, nm1
 330       data two/2.0d0/
 331 c
 332       mup1 = mu + 1
 333       mup2 = mu + 2
 334       diag = -two*eodsq
 335       r2d = two*r4d
 336 c                     1   1
 337 c Compute and store dg /dy
 338       pa(mup1,1) = diag
 339 c
 340 c                     1   2
 341 c Compute and store dg /dy
 342       pa(mu,2) = -r2d*y(2) + eodsq
 343 c
 344       do 20 i = 2,nm1
 345 c
 346 c                     i   i-1
 347 c Compute and store dg /dy
 348         pa(mup2,i-1) = r2d*y(i-1) + eodsq
 349 c
 350 c                     i   i
 351 c Compute and store dg /dy
 352       pa(mup1,i) = diag
 353 c
 354 c                     i   i+1
 355 c Compute and store dg /dy
 356         pa(mu,i+1) = -r2d*y(i+1) + eodsq
 357    20   continue
 358 c
 359 c                     n   n-1
 360 c Compute and store dg /dy
 361       pa(mup2,nm1) = r2d*y(nm1) + eodsq
 362 c
 363 c                     n   n
 364 c Compute and store dg /dy
 365       pa(mup1,n(1)) = diag
 366 c
 367       return
 368 c end of subroutine jacbd for the DLSODI demonstration problem.
 369       end
 370
 371       subroutine jacfl (n, t, y, s, ml, mu, pa, m0)
 372 c This subroutine computes the Jacobian dg/dy = d(g-a*s)/dy
 373 c and stores elements
 374 c   i   j
 375 c dg /dy   in  pa(i,j) in full matrix format.
 376 c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1
 377 c from the Common block test1.
 378 c
 379       integer i, n, m0, ml, mu, nm1
 380       double precision t, y, s, pa, diag, r4d, eodsq, two, r2d
 381       dimension y(*), s(8), pa(m0,*)
 382       dimension n(*)
 383       common /test1/ r4d, eodsq, nm1
 384       data two/2.0d0/
 385 c
 386       diag = -two*eodsq
 387       r2d = two*r4d
 388 c
 389 c                     1   1
 390 c Compute and store dg /dy
 391       pa(1,1) = diag
 392 c
 393 c                     1   2
 394 c Compute and store dg /dy
 395       pa(1,2) = -r2d*y(2) + eodsq
 396 c
 397       do 120  i = 2,nm1
 398 c
 399 c                     i   i-1
 400 c Compute and store dg /dy
 401         pa(i,i-1) = r2d*y(i-1) + eodsq
 402 c
 403 c                     i   i
 404 c Compute and store dg /dy
 405       pa(i,i) = diag
 406 c
 407 c                     i   i+1
 408 c Compute and store dg /dy
 409         pa(i,i+1) = -r2d*y(i+1) + eodsq
 410   120   continue
 411 c
 412 c                     n   n-1
 413 c Compute and store dg /dy
 414       pa(n(1),nm1) = r2d*y(nm1) + eodsq
 415 c
 416 c                     n   n
 417 c Compute and store dg /dy
 418       pa(n(1),n(1)) = diag
 419 c
 420       return
 421 c end of subroutine jacfl for the DLSODI demonstration problem.
 422       end
 423
 424       double precision function elkup (n, y, ewt, itol, rtol, atol)
 425 c This routine looks up approximately correct values of y at t = 0.4,
 426 c ytrue = y9 or y99 depending on whether n = 9 or 99.  These were
 427 c obtained by running DLSODI with very tight tolerances.
 428 c The returned value is
 429 c elkup  =  norm of  ( y - ytrue ) / ( rtol*abs(ytrue) + atol ).
 430 c
 431       integer n, itol, i
 432       double precision y, ewt, rtol, atol, y9, y99
 433 c      double precision y99a, y99b, y99c,
 434 c     1   y99d, y99e, y99f, y99g, dvnorm
 435       dimension y(n), ewt(n), y9(9), y99(99)
 436       dimension rtol(*), atol(*)
 437 c      dimension y99a(16), y99b(16), y99c(16), y99d(16), y99e(16),
 438 c     1          y99f(16), y99g(3)
 439 c      equivalence (y99a(1),y99(1)), (y99b(1),y99(17)),
 440 c     1      (y99c(1),y99(33)), (y99d(1),y99(49)), (y99e(1),y99(65)),
 441 c     1      (y99f(1),y99(81)), (y99g(1),y99(97))
 442       data y9 /
 443      1 1.07001457d-01, 2.77432492d-01, 5.02444616d-01, 7.21037157d-01,
 444      1 9.01670441d-01, 8.88832048d-01, 4.96572850d-01, 9.46924362d-02,
 445      1-6.90855199d-03 /
 446       data y99 /
 447      1 2.05114384d-03, 4.19527452d-03, 6.52533872d-03, 9.13412751d-03,
 448      1 1.21140191d-02, 1.55565301d-02, 1.95516488d-02, 2.41869487d-02,
 449      1 2.95465081d-02, 3.57096839d-02, 4.27498067d-02, 5.07328729d-02,
 450      1 5.97163151d-02, 6.97479236d-02, 8.08649804d-02, 9.30936515d-02,
 451      1 1.06448659d-01, 1.20933239d-01, 1.36539367d-01, 1.53248227d-01,
 452      1 1.71030869d-01, 1.89849031d-01, 2.09656044d-01, 2.30397804d-01,
 453      1 2.52013749d-01, 2.74437805d-01, 2.97599285d-01, 3.21423708d-01,
 454      1 3.45833531d-01, 3.70748792d-01, 3.96087655d-01, 4.21766871d-01,
 455      1 4.47702161d-01, 4.73808532d-01, 5.00000546d-01, 5.26192549d-01,
 456      1 5.52298887d-01, 5.78234121d-01, 6.03913258d-01, 6.29252015d-01,
 457      1 6.54167141d-01, 6.78576790d-01, 7.02400987d-01, 7.25562165d-01,
 458      1 7.47985803d-01, 7.69601151d-01, 7.90342031d-01, 8.10147715d-01,
 459      1 8.28963844d-01, 8.46743353d-01, 8.63447369d-01, 8.79046021d-01,
 460      1 8.93519106d-01, 9.06856541d-01, 9.19058529d-01, 9.30135374d-01,
 461      1 9.40106872d-01, 9.49001208d-01, 9.56853318d-01, 9.63702661d-01,
 462      1 9.69590361d-01, 9.74555682d-01, 9.78631814d-01, 9.81840924d-01,
 463      1 9.84188430d-01, 9.85656465d-01, 9.86196496d-01, 9.85721098d-01,
 464      1 9.84094964d-01, 9.81125395d-01, 9.76552747d-01, 9.70041743d-01,
 465      1 9.61175143d-01, 9.49452051d-01, 9.34294085d-01, 9.15063568d-01,
 466      1 8.91098383d-01, 8.61767660d-01, 8.26550038d-01, 7.85131249d-01,
 467      1 7.37510044d-01, 6.84092540d-01, 6.25748369d-01, 5.63802368d-01,
 468      1 4.99946558d-01, 4.36077986d-01, 3.74091566d-01, 3.15672765d-01,
 469      1 2.62134958d-01, 2.14330497d-01, 1.72640946d-01, 1.37031155d-01,
 470      1 1.07140815d-01, 8.23867920d-02, 6.20562432d-02, 4.53794321d-02,
 471      1 3.15789227d-02, 1.98968820d-02, 9.60472135d-03 /
 472 c
 473       if (n .eq. 99) go to 99
 474 c
 475 c Compute local error tolerance using correct y (n = 9).
 476 c
 477       call dewset( n, itol, rtol, atol, y9, ewt )
 478 c
 479 c Invert ewt and replace y by the error, y - ytrue.
 480 c
 481       do 20  i = 1, 9
 482         ewt(i) = 1.0d0/ewt(i)
 483  20     y(i) = y(i) - y9(i)
 484       go to 200
 485 c
 486 c Compute local error tolerance using correct y (n = 99).
 487 c
 488  99   call dewset( n, itol, rtol, atol, y99, ewt )
 489 c
 490 c Invert ewt and replace y by the error, y - ytrue.
 491 c
 492       do 120  i = 1, 99
 493         ewt(i) = 1.0d0/ewt(i)
 494  120    y(i) = y(i) - y99(i)
 495 c
 496 c Find weighted norm of the error and return.
 497 c
 498  200  elkup = dvnorm (n, y, ewt)
 499       return
 500 c end of function elkup for the DLSODI demonstration program.
 501       end