c----------------------------------------------------------------------- c Demonstration program for the DLSODE package. c This is the version of 14 June 2001. c c This version is in double precision. c c The package is used to solve two simple problems, c one with a full Jacobian, the other with a banded Jacobian, c with all 8 of the appropriate values of mf in each case. c If the errors are too large, or other difficulty occurs, c a warning message is printed. All output is on unit lout = 6. c----------------------------------------------------------------------- external f1, jac1, f2, jac2 integer i, iopar, iopt, iout, istate, itask, itol, iwork, 1 leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter, 2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t, 1 tout, tout1, y dimension y(25), rwork(697), iwork(45) data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/ c nerr = 0 itol = 1 rtol = 0.0d0 atol = 1.0d-6 lrw = 697 liw = 45 iopt = 0 c c First problem c neq = 2 nout = 4 write (lout,110) neq,itol,rtol,atol 110 format(/' Demonstration program for DLSODE package'/// 1 ' Problem 1: Van der Pol oscillator:'/ 2 ' xdotdot - 3*(1 - x**2)*xdot + x = 0, ', 3 ' x(0) = 2, xdot(0) = 0'/ 4 ' neq =',i2/ 5 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) c do 195 meth = 1,2 do 190 miter = 0,3 mf = 10*meth + miter write (lout,120) mf 120 format(///' Solution with mf =',i3// 1 5x,'t x xdot nq h'//) t = 0.0d0 y(1) = 2.0d0 y(2) = 0.0d0 itask = 1 istate = 1 tout = tout1 ero = 0.0d0 do 170 iout = 1,nout call dlsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac1,mf) hu = rwork(11) nqu = iwork(14) write (lout,140) t,y(1),y(2),nqu,hu 140 format(d15.5,d16.5,d14.3,i5,d14.3) if (istate .lt. 0) go to 175 iopar = iout - 2*(iout/2) if (iopar .ne. 0) go to 170 er = abs(y(1))/atol ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,150) 150 format(//' Warning: error exceeds 1000 * tolerance'//) nerr = nerr + 1 endif 170 tout = tout + dtout 175 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (miter .eq. 2) nfea = nfe - neq*nje if (miter .eq. 3) nfea = nfe - nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 180 format(//' Final statistics for this run:'/ 1 ' rwork size =',i4,' iwork size =',i4/ 2 ' number of steps =',i5/ 3 ' number of f-s =',i5/ 4 ' (excluding J-s) =',i5/ 5 ' number of J-s =',i5/ 6 ' error overrun =',d10.2) 190 continue 195 continue c c Second problem c neq = 25 ml = 5 mu = 0 iwork(1) = ml iwork(2) = mu mband = ml + mu + 1 nout = 5 write (lout,210) neq,ml,mu,itol,rtol,atol 210 format(///70('-')/// 1 ' Problem 2: ydot = A * y , where', 2 ' A is a banded lower triangular matrix'/ 3 12x, 'derived from 2-D advection PDE'/ 4 ' neq =',i3,' ml =',i2,' mu =',i2/ 5 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) do 295 meth = 1,2 do 290 miter = 0,5 if (miter .eq. 1 .or. miter .eq. 2) go to 290 mf = 10*meth + miter write (lout,220) mf 220 format(///' Solution with mf =',i3// 1 5x,'t max.err. nq h'//) t = 0.0d0 do 230 i = 2,neq 230 y(i) = 0.0d0 y(1) = 1.0d0 itask = 1 istate = 1 tout = 0.01d0 ero = 0.0d0 do 270 iout = 1,nout call dlsode(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,mf) call edit2(y,t,erm) hu = rwork(11) nqu = iwork(14) write (lout,240) t,erm,nqu,hu 240 format(d15.5,d14.3,i5,d14.3) if (istate .lt. 0) go to 275 er = erm/atol ero = max(ero,er) if (er .gt. 1000.0d0) then write (lout,150) nerr = nerr + 1 endif 270 tout = tout*10.0d0 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (miter .eq. 5) nfea = nfe - mband*nje if (miter .eq. 3) nfea = nfe - nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 290 continue 295 continue write (lout,300) nerr 300 format(////' Number of errors encountered =',i3) stop end subroutine f1 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(neq), ydot(neq) ydot(1) = y(2) ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd double precision t, y, pd dimension y(neq), pd(nrowpd,neq) pd(1,1) = 0.0d0 pd(1,2) = 1.0d0 pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1)) return end subroutine f2 (neq, t, y, ydot) integer neq, i, j, k, ng double precision t, y, ydot, alph1, alph2, d dimension y(neq), ydot(neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ do 10 j = 1,ng do 10 i = 1,ng k = i + (j - 1)*ng d = -2.0d0*y(k) if (i .ne. 1) d = d + y(k-1)*alph1 if (j .ne. 1) d = d + y(k-ng)*alph2 10 ydot(k) = d return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng double precision t, y, pd, alph1, alph2 dimension y(neq), pd(nrowpd,neq) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ mband = ml + mu + 1 mu1 = mu + 1 mu2 = mu + 2 do 10 j = 1,neq pd(mu1,j) = -2.0d0 pd(mu2,j) = alph1 10 pd(mband,j) = alph2 do 20 j = ng,neq,ng 20 pd(mu2,j) = 0.0d0 return end subroutine edit2 (y, t, erm) integer i, j, k, ng double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt dimension y(25) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ erm = 0.0d0 if (t .eq. 0.0d0) return ex = 0.0d0 if (t .le. 30.0d0) ex = exp(-2.0d0*t) a2 = 1.0d0 do 60 j = 1,ng a1 = 1.0d0 do 50 i = 1,ng k = i + (j - 1)*ng yt = t**(i+j-2)*ex*a1*a2 er = abs(y(k)-yt) erm = max(erm,er) a1 = a1*alph1/i 50 continue a2 = a2*alph2/j 60 continue return end Demonstration program for DLSODE package Problem 1: Van der Pol oscillator: xdotdot - 3*(1 - x**2)*xdot + x = 0, x(0) = 2, xdot(0) = 0 neq = 2 itol = 1 rtol = 0.0E+00 atol = 0.1E-05 Solution with mf = 10 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 3 0.123E+00 0.36076E+01 -0.77986E-04 -0.317E+01 5 0.217E-01 0.58224E+01 -0.16801E+01 0.291E+00 3 0.475E-01 0.80372E+01 0.11669E-03 0.317E+01 5 0.234E-01 Final statistics for this run: rwork size = 52 iwork size = 20 number of steps = 297 number of f-s = 352 (excluding J-s) = 352 number of J-s = 0 error overrun = 0.12E+03 Solution with mf = 11 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00 0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01 0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01 0.80372E+01 0.25894E-04 0.317E+01 5 0.190E-01 Final statistics for this run: rwork size = 58 iwork size = 22 number of steps = 203 number of f-s = 281 (excluding J-s) = 281 number of J-s = 29 error overrun = 0.26E+02 Solution with mf = 12 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.121E+00 0.36076E+01 -0.17732E-04 -0.317E+01 5 0.187E-01 0.58224E+01 -0.16801E+01 0.291E+00 6 0.963E-01 0.80372E+01 0.25894E-04 0.317E+01 5 0.190E-01 Final statistics for this run: rwork size = 58 iwork size = 22 number of steps = 203 number of f-s = 339 (excluding J-s) = 281 number of J-s = 29 error overrun = 0.26E+02 Solution with mf = 13 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.739E-01 0.36076E+01 0.34401E-04 -0.317E+01 6 0.260E-01 0.58224E+01 -0.16801E+01 0.291E+00 4 0.133E+00 0.80372E+01 -0.59053E-04 0.317E+01 5 0.205E-01 Final statistics for this run: rwork size = 56 iwork size = 20 number of steps = 198 number of f-s = 315 (excluding J-s) = 289 number of J-s = 26 error overrun = 0.59E+02 Solution with mf = 20 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.549E-01 0.36076E+01 -0.56579E-04 -0.317E+01 5 0.143E-01 0.58224E+01 -0.16801E+01 0.291E+00 4 0.583E-01 0.80372E+01 0.10387E-03 0.317E+01 5 0.149E-01 Final statistics for this run: rwork size = 38 iwork size = 20 number of steps = 289 number of f-s = 321 (excluding J-s) = 321 number of J-s = 0 error overrun = 0.10E+03 Solution with mf = 21 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01 0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01 0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00 0.80372E+01 0.96867E-04 0.317E+01 5 0.142E-01 Final statistics for this run: rwork size = 44 iwork size = 22 number of steps = 262 number of f-s = 345 (excluding J-s) = 345 number of J-s = 30 error overrun = 0.97E+02 Solution with mf = 22 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.676E-01 0.36076E+01 -0.48977E-04 -0.317E+01 5 0.141E-01 0.58224E+01 -0.16801E+01 0.291E+00 5 0.126E+00 0.80372E+01 0.96867E-04 0.317E+01 5 0.142E-01 Final statistics for this run: rwork size = 44 iwork size = 22 number of steps = 262 number of f-s = 405 (excluding J-s) = 345 number of J-s = 30 error overrun = 0.97E+02 Solution with mf = 23 t x xdot nq h 0.13928E+01 0.16801E+01 -0.291E+00 5 0.709E-01 0.36076E+01 -0.46705E-04 -0.317E+01 5 0.139E-01 0.58224E+01 -0.16801E+01 0.291E+00 3 0.719E-01 0.80372E+01 0.54700E-04 0.317E+01 5 0.154E-01 Final statistics for this run: rwork size = 42 iwork size = 20 number of steps = 271 number of f-s = 414 (excluding J-s) = 383 number of J-s = 31 error overrun = 0.55E+02 ---------------------------------------------------------------------- Problem 2: ydot = A * y , where A is a banded lower triangular matrix derived from 2-D advection PDE neq = 25 ml = 5 mu = 0 itol = 1 rtol = 0.0E+00 atol = 0.1E-05 Solution with mf = 10 t max.err. nq h 0.10000E-01 0.556E-06 2 0.766E-02 0.10000E+00 0.655E-05 3 0.249E-01 0.10000E+01 0.274E-05 4 0.520E-01 0.10000E+02 0.114E-05 3 0.117E+00 0.10000E+03 0.221E-05 2 0.262E+00 Final statistics for this run: rwork size = 420 iwork size = 20 number of steps = 524 number of f-s = 552 (excluding J-s) = 552 number of J-s = 0 error overrun = 0.65E+01 Solution with mf = 13 t max.err. nq h 0.10000E-01 0.839E-06 2 0.949E-02 0.10000E+00 0.208E-05 3 0.250E-01 0.10000E+01 0.127E-03 3 0.168E-01 0.10000E+02 0.113E-04 3 0.385E+00 0.10000E+03 0.145E-05 2 0.149E+02 Final statistics for this run: rwork size = 447 iwork size = 20 number of steps = 129 number of f-s = 235 (excluding J-s) = 201 number of J-s = 34 error overrun = 0.13E+03 Solution with mf = 14 t max.err. nq h 0.10000E-01 0.877E-06 2 0.965E-02 0.10000E+00 0.206E-05 3 0.250E-01 0.10000E+01 0.126E-05 5 0.935E-01 0.10000E+02 0.311E-06 6 0.442E+00 0.10000E+03 0.159E-07 2 0.291E+02 Final statistics for this run: rwork size = 697 iwork size = 45 number of steps = 92 number of f-s = 113 (excluding J-s) = 113 number of J-s = 18 error overrun = 0.21E+01 Solution with mf = 15 t max.err. nq h 0.10000E-01 0.877E-06 2 0.965E-02 0.10000E+00 0.206E-05 3 0.250E-01 0.10000E+01 0.126E-05 5 0.935E-01 0.10000E+02 0.311E-06 6 0.442E+00 0.10000E+03 0.160E-07 2 0.291E+02 Final statistics for this run: rwork size = 697 iwork size = 45 number of steps = 92 number of f-s = 221 (excluding J-s) = 113 number of J-s = 18 error overrun = 0.21E+01 Solution with mf = 20 t max.err. nq h 0.10000E-01 0.465E-06 2 0.483E-02 0.10000E+00 0.131E-05 3 0.148E-01 0.10000E+01 0.427E-05 5 0.635E-01 0.10000E+02 0.192E-05 4 0.351E+00 0.10000E+03 0.929E-07 1 0.455E+00 Final statistics for this run: rwork size = 245 iwork size = 20 number of steps = 330 number of f-s = 530 (excluding J-s) = 530 number of J-s = 0 error overrun = 0.43E+01 Solution with mf = 23 t max.err. nq h 0.10000E-01 0.101E-05 2 0.598E-02 0.10000E+00 0.446E-06 3 0.146E-01 0.10000E+01 0.153E-05 5 0.738E-01 0.10000E+02 0.578E-06 4 0.324E+00 0.10000E+03 0.908E-08 1 0.992E+02 Final statistics for this run: rwork size = 272 iwork size = 20 number of steps = 180 number of f-s = 325 (excluding J-s) = 274 number of J-s = 51 error overrun = 0.15E+01 Solution with mf = 24 t max.err. nq h 0.10000E-01 0.104E-05 2 0.608E-02 0.10000E+00 0.463E-06 3 0.146E-01 0.10000E+01 0.247E-05 5 0.666E-01 0.10000E+02 0.828E-06 5 0.391E+00 0.10000E+03 0.384E-09 1 0.108E+03 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 118 number of f-s = 136 (excluding J-s) = 136 number of J-s = 18 error overrun = 0.25E+01 Solution with mf = 25 t max.err. nq h 0.10000E-01 0.104E-05 2 0.608E-02 0.10000E+00 0.463E-06 3 0.146E-01 0.10000E+01 0.247E-05 5 0.666E-01 0.10000E+02 0.828E-06 5 0.391E+00 0.10000E+03 0.384E-09 1 0.108E+03 Final statistics for this run: rwork size = 522 iwork size = 45 number of steps = 118 number of f-s = 244 (excluding J-s) = 136 number of J-s = 18 error overrun = 0.25E+01 Number of errors encountered = 0