c----------------------------------------------------------------------- c demonstration program for the lsoda package. c this is the january 26, 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 lsadem(laout,tape6=laout) c c the package is used to solve two simple problems, c one with a full jacobian, the other with a banded jacobian, c with the 2 appropriate values of jt in each case. 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, jac1, f2, jac2 integer i, iopar, iopt, iout, istate, itask, itol, iwork, 1 jt, leniw, lenrw, liw, lout, lrw, mband, mused, 2 ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t, 1 tout, tout1, tsw, y dimension y(25), rwork(522), iwork(45) data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/ c nerr = 0 itol = 1 rtol = 0.0d0 atol = 1.0d-6 lrw = 522 liw = 45 iopt = 0 c c first problem c neq = 2 nout = 4 write (lout,110) neq,itol,rtol,atol 110 format(1h1/1x,40h demonstration program for lsoda package//// 1 1x,39h problem 1.. van der pol oscillator../ 2 1x,38h xdotdot - 3*(1 - x**2)*xdot + x = 0, , 3 24h x(0) = 2, xdot(0) = 0/ 4 1x,6h neq =,i2/ 5 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,e10.1//) c do 190 jt = 1,2 write (lout,120) jt 120 format(////1x,5h jt =,i3/// 1 6x,1ht,15x,1hx,15x,4hxdot,7x,4hmeth,3x,2hnq,6x,1hh,12x,3htsw//) t = 0.0d0 y(1) = 2.0d0 y(2) = 0.0d0 itask = 1 istate = 1 tout = tout1 ero = 0.0d0 do 170 iout = 1,nout call lsoda(f1,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac1,jt) hu = rwork(11) tsw = rwork(15) nqu = iwork(14) mused = iwork(19) write (lout,140) t,y(1),y(2),mused,nqu,hu,tsw 140 format(1x,e15.5,e16.5,e14.3,2i6,2e14.3) if (istate .lt. 0) go to 175 iopar = iout - 2*(iout/2) if (iopar .ne. 0) go to 170 er = dabs(y(1))/atol ero = dmax1(ero,er) if (er .lt. 1000.0d0) go to 170 write (lout,150) 150 format(//1x,41h warning.. error exceeds 1000 * tolerance//) nerr = nerr + 1 170 tout = tout + dtout 175 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 2) nfea = nfe - neq*nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 180 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,16h error overrun =,e10.2) 190 continue c c second problem c neq = 25 ml = 5 mu = 0 iwork(1) = ml iwork(2) = mu mband = ml + mu + 1 nout = 5 write (lout,210) neq,ml,mu,itol,rtol,atol 210 format(1h1/1x,40h demonstration program for lsoda package//// 1 1x,33h problem 2.. ydot = a * y , where, 2 39h a is a banded lower triangular matrix/ 2 1x,32h derived from 2-d advection pde/ 3 1x,6h neq =,i3,7h ml =,i2,7h mu =,i2/ 4 1x,7h itol =,i3,9h rtol =,e10.1,9h atol =,e10.1//) do 290 jt = 4,5 write (lout,220) jt 220 format(////1x,5h jt =,i3/// 1 6x,1ht,13x,8hmax.err.,5x,4hmeth,3x,2hnq,6x,1hh,12x,3htsw//) t = 0.0d0 do 230 i = 2,neq 230 y(i) = 0.0d0 y(1) = 1.0d0 itask = 1 istate = 1 tout = 0.01d0 ero = 0.0d0 do 270 iout = 1,nout call lsoda(f2,neq,y,t,tout,itol,rtol,atol,itask,istate, 1 iopt,rwork,lrw,iwork,liw,jac2,jt) call edit2(y,t,erm) hu = rwork(11) tsw = rwork(15) nqu = iwork(14) mused = iwork(19) write (lout,240) t,erm,mused,nqu,hu,tsw 240 format(1x,e15.5,e14.3,2i6,2e14.3) if (istate .lt. 0) go to 275 er = erm/atol ero = dmax1(ero,er) if (er .le. 1000.0d0) go to 270 write (lout,150) nerr = nerr + 1 270 tout = tout*10.0d0 275 continue if (istate .lt. 0) nerr = nerr + 1 nst = iwork(11) nfe = iwork(12) nje = iwork(13) lenrw = iwork(17) leniw = iwork(18) nfea = nfe if (jt .eq. 5) nfea = nfe - mband*nje write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero 290 continue 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(2), ydot(2) ydot(1) = y(2) ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1) return end subroutine jac1 (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) = -6.0d0*y(1)*y(2) - 1.0d0 pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1)) return end subroutine f2 (neq, t, y, ydot) integer neq, i, j, k, ng double precision t, y, ydot, alph1, alph2, d dimension y(1), ydot(1) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ do 10 j = 1,ng do 10 i = 1,ng k = i + (j - 1)*ng d = -2.0d0*y(k) if (i .ne. 1) d = d + y(k-1)*alph1 if (j .ne. 1) d = d + y(k-ng)*alph2 10 ydot(k) = d return end subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd) integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng double precision t, y, pd, alph1, alph2 dimension y(1), pd(nrowpd,1) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ mband = ml + mu + 1 mu1 = mu + 1 mu2 = mu + 2 do 10 j = 1,neq pd(mu1,j) = -2.0d0 pd(mu2,j) = alph1 10 pd(mband,j) = alph2 do 20 j = ng,neq,ng 20 pd(mu2,j) = 0.0d0 return end subroutine edit2 (y, t, erm) integer i, j, k, ng double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt dimension y(25) data alph1/1.0d0/, alph2/1.0d0/, ng/5/ erm = 0.0d0 if (t .eq. 0.0d0) return ex = 0.0d0 if (t .le. 30.0d0) ex = dexp(-2.0d0*t) a2 = 1.0d0 do 60 j = 1,ng a1 = 1.0d0 do 50 i = 1,ng k = i + (j - 1)*ng yt = t**(i+j-2)*ex*a1*a2 er = dabs(y(k)-yt) erm = dmax1(erm,er) a1 = a1*alph1/dfloat(i) 50 continue a2 = a2*alph2/dfloat(j) 60 continue return end