c----------------------------------------------------------------------- c Demonstration program for the SLSODE package. c This is the version of 12 June 2001. c c This version is in single 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 real 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.39283880203e0/, dtout/2.214773875e0/ c nerr = 0 itol = 1 rtol = 0.0e0 atol = 1.0e-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 SLSODE 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 =',e10.1,' atol =',e10.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.0e0 y(1) = 2.0e0 y(2) = 0.0e0 itask = 1 istate = 1 tout = tout1 ero = 0.0e0 do 170 iout = 1,nout call slsode(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(e15.5,e16.5,e14.3,i5,e14.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.0e0) 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 =',e10.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 =',e10.1,' atol =',e10.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.0e0 do 230 i = 2,neq 230 y(i) = 0.0e0 y(1) = 1.0e0 itask = 1 istate = 1 tout = 0.01e0 ero = 0.0e0 do 270 iout = 1,nout call slsode(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(e15.5,e14.3,i5,e14.3) if (istate .lt. 0) go to 275 er = erm/atol ero = max(ero,er) if (er .gt. 1000.0e0) then write (lout,150) nerr = nerr + 1 endif 270 tout = tout*10.0e0 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 real t, y, ydot dimension y(neq), ydot(neq) ydot(1) = y(2) ydot(2) = 3.0e0*(1.0e0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd real t, y, pd dimension y(neq), pd(nrowpd,neq) pd(1,1) = 0.0e0 pd(1,2) = 1.0e0 pd(2,1) = -6.0e0*y(1)*y(2) - 1.0e0 pd(2,2) = 3.0e0*(1.0e0 - y(1)*y(1)) return end subroutine f2 (neq, t, y, ydot) integer neq, i, j, k, ng real t, y, ydot, alph1, alph2, d dimension y(neq), ydot(neq) data alph1/1.0e0/, alph2/1.0e0/, ng/5/ do 10 j = 1,ng do 10 i = 1,ng k = i + (j - 1)*ng d = -2.0e0*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 real t, y, pd, alph1, alph2 dimension y(neq), pd(nrowpd,neq) data alph1/1.0e0/, alph2/1.0e0/, ng/5/ mband = ml + mu + 1 mu1 = mu + 1 mu2 = mu + 2 do 10 j = 1,neq pd(mu1,j) = -2.0e0 pd(mu2,j) = alph1 10 pd(mband,j) = alph2 do 20 j = ng,neq,ng 20 pd(mu2,j) = 0.0e0 return end subroutine edit2 (y, t, erm) integer i, j, k, ng real y, t, erm, alph1, alph2, a1, a2, er, ex, yt dimension y(25) data alph1/1.0e0/, alph2/1.0e0/, ng/5/ erm = 0.0e0 if (t .eq. 0.0e0) return ex = 0.0e0 if (t .le. 30.0e0) ex = exp(-2.0e0*t) a2 = 1.0e0 do 60 j = 1,ng a1 = 1.0e0 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