1c----------------------------------------------------------------------- 2c Demonstration program for the DLSODAR package. 3c This is the version of 14 June 2001. 4c 5c This version is in double precision. 6c 7c The DLSODAR package is used to solve two simple problems, 8c one nonstiff and one intermittently stiff. 9c If the errors are too large, or other difficulty occurs, 10c a warning message is printed. All output is on unit lout = 6. 11c----------------------------------------------------------------------- 12 external f1, gr1, f2, jac2, gr2 13 integer iopt, iout, istate, itask, itol, iwork, jroot, jt, 14 1 kroot, leniw, lenrw, liw, lrw, lout, neq, nerr, ng, 15 2 nfe, nfea, nge, nje, nst 16 double precision atol, er, ero, errt, rtol, rwork, 17 1 t, tout, tzero, y, yt 18 dimension y(2), atol(2), rwork(57), iwork(22), jroot(2) 19 data lout/6/ 20c 21 nerr = 0 22c----------------------------------------------------------------------- 23c First problem. 24c The initial value problem is: 25c dy/dt = ((2*log(y) + 8)/t - 5)*y, y(1) = 1, 1 .le. t .le. 6 26c The solution is y(t) = exp(-t**2 + 5*t - 4) 27c The two root functions are: 28c g1 = ((2*log(y)+8)/t - 5)*y (= dy/dt) (with root at t = 2.5), 29c g2 = log(y) - 2.2491 (with roots at t = 2.47 and 2.53) 30c----------------------------------------------------------------------- 31c Set all input parameters and print heading. 32 neq = 1 33 y(1) = 1.0d0 34 t = 1.0d0 35 tout = 2.0d0 36 itol = 1 37 rtol = 1.0d-6 38 atol(1) = 1.0d-6 39 itask = 1 40 istate = 1 41 iopt = 0 42 lrw = 44 43 liw = 21 44 jt = 2 45 ng = 2 46 write (lout,110) itol,rtol,atol(1),jt 47 110 format(/' Demonstration program for DLSODAR package'//// 48 1 ' First problem'/// 49 2 ' Problem is dy/dt = ((2*log(y)+8)/t - 5)*y, y(1) = 1'// 50 3 ' Solution is y(t) = exp(-t**2 + 5*t - 4)'// 51 4 ' Root functions are:'/ 52 5 10x,' g1 = dy/dt (root at t = 2.5)'/ 53 6 10x,' g2 = log(y) - 2.2491 (roots at t = 2.47 and t = 2.53)'// 54 7 ' itol =',i3,' rtol =',d10.1,' atol =',d10.1// 55 8 ' jt =',i3///) 56c 57c Call DLSODAR in loop over tout values 2,3,4,5,6. 58 ero = 0.0d0 59 do 180 iout = 1,5 60 120 continue 61 call dlsodar(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 62 1 iopt,rwork,lrw,iwork,liw,jdum,jt,gr1,ng,jroot) 63c 64c Print y and error in y, and print warning if error too large. 65 yt = exp(-t*t + 5.0d0*t - 4.0d0) 66 er = y(1) - yt 67 write (lout,130) t,y(1),er 68 130 format(' At t =',d15.7,5x,'y =',d15.7,5x,'error =',d12.4) 69 if (istate .lt. 0) go to 185 70 er = abs(er)/(rtol*abs(y(1)) + atol(1)) 71 ero = max(ero,er) 72 if (er .gt. 1000.0d0) then 73 write (lout,140) 74 140 format(//' Warning: error exceeds 1000 * tolerance'//) 75 nerr = nerr + 1 76 endif 77 if (istate .ne. 3) go to 175 78c 79c If a root was found, write results and check root location. 80c Then reset istate to 2 and return to DLSODAR call. 81 write (lout,150) t,jroot(1),jroot(2) 82 150 format(/' Root found at t =',d15.7,5x,'jroot =',2i5) 83 if (jroot(1) .eq. 1) errt = t - 2.5d0 84 if (jroot(2) .eq. 1 .and. t .le. 2.5d0) errt = t - 2.47d0 85 if (jroot(2) .eq. 1 .and. t .gt. 2.5d0) errt = t - 2.53d0 86 write (lout,160) errt 87 160 format(' Error in t location of root is',d12.4/) 88 if (abs(errt) .gt. 1.0d-3) then 89 write (lout,170) 90 170 format(//' Warning: root error exceeds 1.0d-3'//) 91 nerr = nerr + 1 92 endif 93 istate = 2 94 go to 120 95c 96c If no root found, increment tout and loop back. 97 175 tout = tout + 1.0d0 98 180 continue 99c 100c Problem complete. Print final statistics. 101 185 continue 102 if (istate .lt. 0) nerr = nerr + 1 103 nst = iwork(11) 104 nfe = iwork(12) 105 nje = iwork(13) 106 nge = iwork(10) 107 lenrw = iwork(17) 108 leniw = iwork(18) 109 nfea = nfe 110 if (jt .eq. 2) nfea = nfe - neq*nje 111 write (lout,190) lenrw,leniw,nst,nfe,nfea,nje,nge,ero 112 190 format(//' Final statistics for this run:'/ 113 1 ' rwork size =',i4,' iwork size =',i4/ 114 2 ' number of steps =',i5/ 115 3 ' number of f-s =',i5/ 116 4 ' (excluding j-s) =',i5/ 117 5 ' number of j-s =',i5/ 118 6 ' number of g-s =',i5/ 119 7 ' error overrun =',d10.2) 120c 121c----------------------------------------------------------------------- 122c Second problem (Van der Pol oscillator). 123c The initial value problem is (after reduction of 2nd order ODE): 124c dy1/dt = y2, dy2/dt = 100*(1 - y1**2)*y2 - y1, 125c y1(0) = 2, y2(0) = 0, 0 .le. t .le. 200 126c The root function is g = y1. 127c An analytic solution is not known, but the zeros of y1 are known 128c to 15 figures for purposes of checking the accuracy. 129c----------------------------------------------------------------------- 130c Set tolerance parameters and print heading. 131 itol = 2 132 rtol = 1.0d-6 133 atol(1) = 1.0d-6 134 atol(2) = 1.0d-4 135 write (lout,200) itol,rtol,atol(1),atol(2) 136 200 format(////80('*')//' Second problem (Van der Pol oscillator)'// 137 1 ' Problem is dy1/dt = y2, dy2/dt = 100*(1-y1**2)*y2 - y1'/ 138 2 ' y1(0) = 2, y2(0) = 0'// 139 3 ' Root function is g = y1'// 140 4 ' itol =',i3,' rtol =',d10.1,' atol =',2d10.1) 141c 142c Loop over jt = 1, 2. Set remaining parameters and print jt. 143 do 290 jt = 1,2 144 neq = 2 145 y(1) = 2.0d0 146 y(2) = 0.0d0 147 t = 0.0d0 148 tout = 20.0d0 149 itask = 1 150 istate = 1 151 iopt = 0 152 lrw = 57 153 liw = 22 154 ng = 1 155 write (lout,210) jt 156 210 format(///' Solution with jt =',i2//) 157c 158c Call DLSODAR in loop over tout values 20,40,...,200. 159 do 270 iout = 1,10 160 220 continue 161 call dlsodar(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 162 1 iopt,rwork,lrw,iwork,liw,jac2,jt,gr2,ng,jroot) 163c 164c Print y1 and y2. 165 write (lout,230) t,y(1),y(2) 166 230 format(' At t =',d15.7,5x,'y1 =',d15.7,5x,'y2 =',d15.7) 167 if (istate .lt. 0) go to 275 168 if (istate .ne. 3) go to 265 169c 170c If a root was found, write results and check root location. 171c Then reset istate to 2 and return to DLSODAR call. 172 write (lout,240) t 173 240 format(/' Root found at t =',d15.7) 174 kroot = int(t/81.2d0 + 0.5d0) 175 tzero = 81.17237787055d0 + (kroot-1)*81.41853556212d0 176 errt = t - tzero 177 write (lout,250) errt 178 250 format(' Error in t location of root is',d12.4//) 179 if (abs(errt) .gt. 1.0d-1) then 180 write (lout,260) 181 260 format(//' Warning: root error exceeds 1.0d-1'//) 182 nerr = nerr + 1 183 endif 184 istate = 2 185 go to 220 186c 187c If no root found, increment tout and loop back. 188 265 tout = tout + 20.0d0 189 270 continue 190c 191c Problem complete. Print final statistics. 192 275 continue 193 if (istate .lt. 0) nerr = nerr + 1 194 nst = iwork(11) 195 nfe = iwork(12) 196 nje = iwork(13) 197 nge = iwork(10) 198 lenrw = iwork(17) 199 leniw = iwork(18) 200 nfea = nfe 201 if (jt .eq. 2) nfea = nfe - neq*nje 202 write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,nge 203 280 format(//' Final statistics for this run:'/ 204 1 ' rwork size =',i4,' iwork size =',i4/ 205 2 ' number of steps =',i5/ 206 3 ' number of f-s =',i5/ 207 4 ' (excluding j-s) =',i5/ 208 5 ' number of j-s =',i5/ 209 6 ' number of g-s =',i5) 210 290 continue 211c 212 write (lout,300) nerr 213 300 format(///' Total number of errors encountered =',i3) 214 stop 215 end 216 217 subroutine f1 (neq, t, y, ydot) 218 integer neq 219 double precision t, y, ydot 220 dimension y(1), ydot(1) 221 ydot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1) 222 return 223 end 224 225 subroutine gr1 (neq, t, y, ng, groot) 226 integer neq, ng 227 double precision t, y, groot 228 dimension y(1), groot(2) 229 groot(1) = ((2.0d0*log(y(1)) + 8.0d0)/t - 5.0d0)*y(1) 230 groot(2) = log(y(1)) - 2.2491d0 231 return 232 end 233 234 subroutine f2 (neq, t, y, ydot) 235 integer neq 236 double precision t, y, ydot 237 dimension y(2), ydot(2) 238 ydot(1) = y(2) 239 ydot(2) = 100.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) 240 return 241 end 242 243 subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) 244 integer neq, ml, mu, nrowpd 245 double precision t, y, pd 246 dimension y(2), pd(nrowpd,2) 247 pd(1,1) = 0.0d0 248 pd(1,2) = 1.0d0 249 pd(2,1) = -200.0d0*y(1)*y(2) - 1.0d0 250 pd(2,2) = 100.0d0*(1.0d0 - y(1)*y(1)) 251 return 252 end 253 254 subroutine gr2 (neq, t, y, ng, groot) 255 integer neq, ng 256 double precision t, y, groot 257 dimension y(2), groot(1) 258 groot(1) = y(1) 259 return 260 end 261