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