c----------------------------------------------------------------------- c demonstration program for the lsodar package. c this is the january 27, 1982 version. c this version is in double precision. c c for computer systems requiring a program card, the following (with c the c in column 1 removed) may be used.. c program lardem(larout,tape6=larout) c c the lsodar package is used to solve two simple problems, c one nonstiff and one intermittently stiff. 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, gr1, f2, jac2, gr2 integer iopt, iout, istate, itask, itol, iwork, jroot, jt, 1 kroot, leniw, lenrw, liw, lrw, lout, neq, nerr, ng, 2 nfe, nfea, nge, nje, nst double precision atol, er, ero, errt, rtol, rwork, 1 t, tout, tzero, y, yt dimension y(2), atol(2), rwork(57), iwork(22), jroot(2) data lout/6/ c nerr = 0 c----------------------------------------------------------------------- c first problem. c the initial value problem is.. c dy/dt = ((2*log(y) + 8)/t - 5)*y, y(1) = 1, 1 .le. t .le. 6 c the solution is y(t) = exp(-t**2 + 5*t - 4) c the two root functions are.. c g1 = ((2*log(y)+8)/t - 5)*y (= dy/dt) (with root at t = 2.5), c g2 = log(y) - 2.2491 (with roots at t = 2.47 and 2.53) c----------------------------------------------------------------------- c set all input parameters and print heading. neq = 1 y(1) = 1.0d0 t = 1.0d0 tout = 2.0d0 itol = 1 rtol = 0.0d0 atol(1) = 1.0d-6 itask = 1 istate = 1 iopt = 0 lrw = 44 liw = 21 jt = 2 ng = 2 write (lout,110) itol,rtol,atol(1),jt 110 format(1h1/1x,41h demonstration program for lsodar package///// 1 1x,14h first problem/// 2 1x,54h problem is dy/dt = ((2*log(y)+8)/t - 5)*y, y(1) = 1// 3 1x,41h solution is y(t) = exp(-t**2 + 5*t - 4)// 4 1x,21h root functions are../ 5 10x,30h g1 = dy/dt (root at t = 2.5)/ 6 10x,55h g2 = log(y) - 6.2491 (roots at t = 2.47 and t = 2.53)// 7 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,e10.1/ 8 1x,5h jt =,i3/////) c c call lsodar in loop over tout values 2,3,4,5,6. ero = 0.0d0 do 180 iout = 1,5 120 continue call lsodar(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jdum,jt,gr1,ng,jroot) c c print y and error in y, and print warning if error too large. yt = dexp(-t*t + 5.0d0*t - 4.0d0) er = y(1) - yt write (lout,130) t,y(1),er 130 format(1x,7h at t =,e15.7,5x,3hy =,e15.7,5x,7herror =,e12.4) if (istate .lt. 0) go to 185 er = dabs(er)/atol(1) ero = dmax1(ero,er) if (er .lt. 1000.0d0) go to 140 write (lout,135) 135 format(//1x,41h warning.. error exceeds 1000 * tolerance//) nerr = nerr + 1 140 continue if (istate .ne. 3) go to 175 c c if a root was found, write results and check root location. c then reset istate to 2 and return to lsodar call. write (lout,150) t,jroot(1),jroot(2) 150 format(/1x,18h root found at t =,e15.7,5x,7hjroot =,2i5) if (jroot(1) .eq. 1) errt = t - 2.5d0 if (jroot(2) .eq. 1 .and. t .le. 2.5d0) errt = t - 2.47d0 if (jroot(2) .eq. 1 .and. t .gt. 2.5d0) errt = t - 2.53d0 write (lout,160) errt 160 format(1x,31h error in t location of root is,e12.4//) if (dabs(errt) .lt. 1.0d-3) go to 170 write (lout,165) 165 format(//1x,36h warning.. root error exceeds 1.0e-3//) nerr = nerr + 1 170 istate = 2 go to 120 c c if no root found, increment tout and loop back. 175 tout = tout + 1.0d0 180 continue c c problem complete. print final statistics. 185 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) nge = iwork(10) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,190) lenrw,leniw,nst,nfe,nfea,nje,nge,ero 190 format(//1x,32h final statistics for this run../ 1 1x,13h rwork size =,i4,15h iwork size =,i4/ 2 1x,18h number of steps =,i5/ 3 1x,18h number of f-s =,i5/ 4 1x,18h (excluding j-s) =,i5/ 5 1x,18h number of j-s =,i5/ 6 1x,18h number of g-s =,i5/ 7 1x,16h error overrun =,e10.2) c c----------------------------------------------------------------------- c second problem (van der pol oscillator). c the initial value problem is.. c dy1/dt = y2, dy2/dt = 100*(1 - y1**2)*y2 - y1, c y1(0) = 2, y2(0) = 0, 0 .le. t .le. 200 c the root function is g = y1. c an analytic solution is not known, but the zeros of y1 are known c to 15 figures for purposes of checking the accuracy. c----------------------------------------------------------------------- c set tolerance parameters and print heading. itol = 2 rtol = 1.0d-6 atol(1) = 1.0d-6 atol(2) = 1.0d-4 write (lout,200) itol,rtol,atol(1),atol(2) 200 format(1h1/1x,40h second problem (van der pol oscillator)/// 1 1x,56h problem is dy1/dt = y2, dy2/dt = 100*(1-y1**2)*y2 - y1/ 2 1x,33h y1(0) = 2, y2(0) = 0// 3 1x,25h root function is g = y1// 4 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,2e10.1///) c c loop over jt = 1, 2. set remaining parameters and print jt. do 290 jt = 1,2 neq = 2 y(1) = 2.0d0 y(2) = 0.0d0 t = 0.0d0 tout = 20.0d0 itask = 1 istate = 1 iopt = 0 lrw = 57 liw = 22 ng = 1 write (lout,210) jt 210 format(////////1x,19h solution with jt =,i2////) c c call lsodar in loop over tout values 20,40,...,200. do 270 iout = 1,10 220 continue call lsodar(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,jt,gr2,ng,jroot) c c print y1 and y2. write (lout,230) t,y(1),y(2) 230 format(1x,7h at t =,e15.7,5x,4hy1 =,e15.7,5x,4hy2 =,e15.7) if (istate .lt. 0) go to 275 if (istate .ne. 3) go to 265 c c if a root was found, write results and check root location. c then reset istate to 2 and return to lsodar call. write (lout,240) t 240 format(/1x,18h root found at t =,e15.7) kroot = ifix(sngl(t/81.2d0 + 0.5d0)) tzero = 81.17237787055d0 + dfloat(kroot-1)*81.41853556212d0 errt = t - tzero write (lout,250) errt 250 format(1x,31h error in t location of root is,e12.4//) if (errt .lt. 1.0d0) go to 260 write (lout,255) 255 format(//1x,33h warning.. root error exceeds 1.0//) nerr = nerr + 1 260 istate = 2 go to 220 c c if no root found, increment tout and loop back. 265 tout = tout + 20.0d0 270 continue c c problem complete. print final statistics. 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) nge = iwork(10) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,280) lenrw,leniw,nst,nfe,nfea,nje,nge 280 format(//1x,32h final statistics for this run../ 1 1x,13h rwork size =,i4,15h iwork size =,i4/ 2 1x,18h number of steps =,i5/ 3 1x,18h number of f-s =,i5/ 4 1x,18h (excluding j-s) =,i5/ 5 1x,18h number of j-s =,i5/ 6 1x,18h number of g-s =,i5) 290 continue c c write (lout,300) nerr 300 format(////1x,31h number of errors encountered =,i3) stop end subroutine f1 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(1), ydot(1) ydot(1) = ((2.0d0*dlog(y(1)) + 8.0d0)/t - 5.0d0)*y(1) return end subroutine gr1 (neq, t, y, ng, groot) integer neq, ng double precision t, y, groot dimension y(1), groot(2) groot(1) = ((2.0d0*dlog(y(1)) + 8.0d0)/t - 5.0d0)*y(1) groot(2) = dlog(y(1)) - 2.2491d0 return end subroutine f2 (neq, t, y, ydot) integer neq double precision t, y, ydot dimension y(2), ydot(2) ydot(1) = y(2) ydot(2) = 100.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd double precision t, y, pd dimension y(2), pd(nrowpd,2) pd(1,1) = 0.0d0 pd(1,2) = 1.0d0 pd(2,1) = -200.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 100.0d0*(1.0d0 - y(1)*y(1)) return end subroutine gr2 (neq, t, y, ng, groot) integer neq, ng double precision t, y, groot dimension y(2), groot(1) groot(1) = y(1) return end