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