1c----------------------------------------------------------------------- 2c Demonstration program for the DLSODE package. 3c This is the version of 14 June 2001. 4c 5c This version is in double precision. 6c 7c The package is used to solve two simple problems, 8c one with a full Jacobian, the other with a banded Jacobian, 9c with all 8 of the appropriate values of mf in each case. 10c If the errors are too large, or other difficulty occurs, 11c a warning message is printed. All output is on unit lout = 6. 12c----------------------------------------------------------------------- 13 external f1, jac1 14 integer i, iopar, iopt, iout, istate, itask, itol, iwork, 15 1 leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter, 16 2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst 17 double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t, 18 1 tout, tout1, y 19 dimension y(25), rwork(697), iwork(45) 20 data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/ 21c 22 nerr = 0 23 itol = 1 24 rtol = 0.0d0 25 atol = 1.0d-6 26 lrw = 697 27 liw = 45 28 iopt = 0 29c 30c First problem 31c 32 neq = 2 33 nout = 4 34 write (lout,110) neq,itol,rtol,atol 35 110 format(/' Demonstration program for DLSODE package'/// 36 1 ' Problem 1: Van der Pol oscillator:'/ 37 2 ' xdotdot - 3*(1 - x**2)*xdot + x = 0, ', 38 3 ' x(0) = 2, xdot(0) = 0'/ 39 4 ' neq =',i2/ 40 5 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1//) 41c 42 do 195 meth = 1,2 43 do 190 miter = 0,3 44 mf = 10*meth + miter 45 write (lout,120) mf 46 120 format(///' Solution with mf =',i3// 47 1 5x,'t x xdot nq h'//) 48 t = 0.0d0 49 y(1) = 2.0d0 50 y(2) = 0.0d0 51 itask = 1 52 istate = 1 53 tout = tout1 54 ero = 0.0d0 55 do 170 iout = 1,nout 56 call dlsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 57 1 iopt,rwork,lrw,iwork,liw,jac1,mf) 58 hu = rwork(11) 59 nqu = iwork(14) 60 write (lout,140) t,y(1),y(2),nqu,hu 61 140 format(d15.5,d16.5,d14.3,i5,d14.3) 62 if (istate .lt. 0) go to 175 63 iopar = iout - 2*(iout/2) 64 if (iopar .ne. 0) go to 170 65 er = abs(y(1))/atol 66 ero = max(ero,er) 67 if (er .gt. 1000.0d0) then 68 write (lout,150) 69 150 format(//' Warning: error exceeds 1000 * tolerance'//) 70 nerr = nerr + 1 71 endif 72 170 tout = tout + dtout 73 175 continue 74 if (istate .lt. 0) nerr = nerr + 1 75 nst = iwork(11) 76 nfe = iwork(12) 77 nje = iwork(13) 78 lenrw = iwork(17) 79 leniw = iwork(18) 80 nfea = nfe 81 if (miter .eq. 2) nfea = nfe - neq*nje 82 if (miter .eq. 3) nfea = nfe - nje 83 write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 84 180 format(//' Final statistics for this run:'/ 85 1 ' rwork size =',i4,' iwork size =',i4/ 86 2 ' number of steps =',i5/ 87 3 ' number of f-s =',i5/ 88 4 ' (excluding J-s) =',i5/ 89 5 ' number of J-s =',i5/ 90 6 ' error overrun =',d10.2) 91 190 continue 92 195 continue 93 94 write (lout,300) nerr 95 300 format(////' Number of errors encountered =',i3) 96 stop 97 end 98 99c$$$ subroutine f1 (neq, t, y, ydot) 100c$$$ integer neq 101c$$$ double precision t, y, ydot 102c$$$ dimension y(neq), ydot(neq) 103c$$$ ydot(1) = y(2) 104c$$$ ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) 105c$$$ return 106c$$$ end 107c$$$ 108c$$$ subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd) 109c$$$ integer neq, ml, mu, nrowpd 110c$$$ double precision t, y, pd 111c$$$ dimension y(neq), pd(nrowpd,neq) 112c$$$ pd(1,1) = 0.0d0 113c$$$ pd(1,2) = 1.0d0 114c$$$ pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0 115c$$$ pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1)) 116c$$$ return 117c$$$ end 118