Source file: tlsoda.f

c===========================================================
c     tlsoda: Demo program which uses ODEPACK routine LSODA
c     to solve the second-order ODE
c
c        u''(t) = -u(t),  tmin <= t <= tmax
c
c     (' = d/dt), with initial conditions
c
c        u(tmin) = utmin,   u'(tmin) = dutmin
c
c     The exact solution is
c
c        u_xct(t) = dutmin * sin(t) + utmin * cos(t)
c
c     Output to standard out is
c
c        t_it    u_it   [u_xct - u]_it
c
c     it = 1, 2, ... nout, where
c
c        t_1     = tmin
c        t_ntout = tmax
c        ntout   = 2**olevel + 1
c
c     Output to standard error is rms error.
c===========================================================
      program         tlsoda
      implicit        none
      integer         iargc,        i4arg
      real*8          r8arg
c-----------------------------------------------------------
c     Command-line arguments:
c
c     tmin:    Initial integration time
c     tmax:    Final integration time
c     utmin:   Initial value: u(tmin)
c     dutmin:  Initial value: u'(tmin)
c     tol:     Error tolerance (this program uses LSODA's
c              pure absolute error control)
c     olevel:  Output level:
c              dtout = (tmax-tmin)/(2**olevel+1)
c-----------------------------------------------------------
      real*8         tmin,       tmax,       utmin,
     &               dutmin,     tol
      integer        olevel
      real*8         r8_never
      parameter    ( r8_never = -1.0d-60 )
c===========================================================
c     Start of LSODA declarations
c===========================================================
c-----------------------------------------------------------
c     Note that 'fcn' and 'jac' are user supplied subroutines
c     (not functions) which evaluate the RHSs of the ODEs and
c     the Jacobian of the system.  Under normal operation,
c     (as in this case), the Jacobian evaluator can be a
c     'dummy' routine.
c-----------------------------------------------------------
      external       fcn,        jac
c-----------------------------------------------------------
c     Number of ODEs (when written in canonical first order
c     form).
c-----------------------------------------------------------
      integer         neq
      parameter     ( neq = 2 )
c-----------------------------------------------------------
c     y(neq): Storage for approximate solution
c     t:      Initial time for LSODA integration sub-interval
c     tout:   Final time for LSODA integration sub-interval
c-----------------------------------------------------------
      real*8         y(neq)
      real*8         t,          tout
c-----------------------------------------------------------
c     Tolerance parameters:
c

                                                                 c     The following comment block is extracted from the
                                                                 c     LSODA documentation.
                                                                 c-----------------------------------------------------------
                                                                 c rtol   = relative tolerance parameter (scalar).
                                                                 c atol   = absolute tolerance parameter (scalar or array).
                                                                 c   the estimated local error in y(i) will be controlled so
                                                                 c   as to be less than
                                                                 c        ewt(i) = rtol*abs(y(i)) + atol     if itol = 1, or
                                                                 c        ewt(i) = rtol*abs(y(i)) + atol(i)  if itol = 2.
                                                                 c   thus the local error test passes if, in each component,
                                                                 c   either the absolute error is less than atol (or atol(i)),
                                                                 c   or the relative error is less than rtol.
                                                                 c   use rtol = 0.0 for pure absolute error control, and
                                                                 c   use atol = 0.0 (or atol(i) = 0.0) for pure relative error
                                                                 c   control.  CAUTION.. actual (global) errors may exceed
                                                                 c   these local tolerances, so choose them CONSERVATIVELY.
                                                                 c-----------------------------------------------------------
                                                                       real*8         rtol,       atol
                                                                       integer        itol
                                                                 c-----------------------------------------------------------
                                                                 c     Control parameters and return code (see below).
                                                                 c-----------------------------------------------------------
                                                                       integer        itask,      istate,     iopt
                                                                 c-----------------------------------------------------------
                                                                 c     Work arrays.
                                                                 c-----------------------------------------------------------
                                                                       integer        lrw
                                                                       parameter    ( lrw = 22 + neq * 16 )
                                                                       real*8         rwork(lrw)
                                                                       integer        liw
                                                                       parameter    ( liw = 20 + neq )
                                                                       integer        iwork(liw)
                                                                 c-----------------------------------------------------------
                                                                 c     'jt' defines which type of Jacobian is supplied or
                                                                 c     computed.
                                                                 c-----------------------------------------------------------
                                                                       integer        jt
                                                                 c===========================================================
                                                                 c     End of LSODA declarations
                                                                 c===========================================================
                                                                 c-----------------------------------------------------------
                                                                 c     Miscellaneous variables
                                                                 c-----------------------------------------------------------
                                                                       real*8          dtout,    err,     rmserr
                                                                       integer         it,       ntout
                                                                 c-----------------------------------------------------------
                                                                 c     Argument parsing.
                                                                 c-----------------------------------------------------------
                                                                       if( iargc() .ne. 6 )  go to 900
                                                                       tmin   = r8arg(1,r8_never)
                                                                       tmax   = r8arg(2,r8_never)
                                                                       utmin  = r8arg(3,r8_never)
                                                                       dutmin = r8arg(4,r8_never)
                                                                       tol    = r8arg(5,r8_never)
                                                                       olevel = i4arg(6,-1)
                                                                       if( tmin  .eq. r8_never .or. tmax   .eq. r8_never .or.
                                                                      &    utmin .eq. r8_never .or. dutmin .eq. r8_never .or.
                                                                      &    tol   .eq. r8_never .or. olevel .lt. 0 )
                                                                      &  go to 900
                                                                 c-----------------------------------------------------------
                                                                 c     Set LSODA parameters ... see LSODA documentation
                                                                 c     for fuller description.
                                                                 c-----------------------------------------------------------
                                                                       itol   = 1         ! Indicates that 'atol' is scalar
                                                                       rtol   = 0.0d0     ! Use pure absolute tolerance
                                                                       atol   = tol       ! Absolute tolerance
                                                                       itask  = 1         ! Normal computation
                                                                       iopt   = 0         ! Indicates no optional inputs
                                                                       jt     = 2         ! Jacobian type
                                                                 c-----------------------------------------------------------
                                                                 c     Compute number of output times and output interval,
                                                                 c     and intialize sub-interval start time and solution

c     estimate.
c-----------------------------------------------------------
      ntout  = 2**olevel + 1
      dtout  = (tmax - tmin) / (ntout - 1)
      t      = tmin
      y(1)   = utmin
      y(2)   = dutmin
c-----------------------------------------------------------
c     Output initial solution and error and initialize
c     rms error.
c-----------------------------------------------------------
      err = dutmin * sin(t) + utmin * cos(t) - y(1)
      write(*,*) t, y(1), err
      rmserr = err**2
c-----------------------------------------------------------
c     Loop over requested output times ...
c-----------------------------------------------------------
      do it = 2, ntout
c-----------------------------------------------------------
c        Set final integration time for current interval ...
c-----------------------------------------------------------
         tout = t + dtout
c-----------------------------------------------------------
c        Call lsoda to integrate system on [t ... tout]
c
c        'istate' should always be set to 1 prior to
c        invocation; LSODA also uses it as a return code.
c
c        Also note that LSODA overwrites 't' with 'tout'
c        if integration is successful.
c-----------------------------------------------------------
         istate = 1
         call lsoda(fcn,neq,y,t,tout,
     &              itol,rtol,atol,itask,
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)
c-----------------------------------------------------------
c        Check return code and exit with error message if
c        there was trouble.
c-----------------------------------------------------------
         if( istate .lt. 0 ) then
            write(0,1000) istate, it, ntout, t, t + dtout
1000        format(/' sode: Error return ',i2,
     &              ' from integrator LSODA.'/
     &              ' sode: At output time ',i5,' of ',i5/
     &              ' sode: Interval ',1pe11.3,' .. ',1pe11.3/)
            go to 500
         end if
c-----------------------------------------------------------
c        Output the solution and error, and RMS error
c        accumulator.
c-----------------------------------------------------------
         err = dutmin * sin(t) + utmin * cos(t) - y(1)
         write(*,*) t, y(1), err
         rmserr = rmserr + err**2
      end do
      rmserr = sqrt(rmserr / ntout)
      write(0,*) 'rmserr: ', rmserr
 500  continue
      stop
 900  continue
         write(0,*) 'usage: tlsoda <tmin> <tmax> <u_tmin> '//
     &              '<du_tmin> <tol> <olevel>'
      stop
      end
c===========================================================
c     Implements differential equations:
c
c     u'' = -u
c
c     y(1) := u
c     y(2) := u'
c
c     y(1)' :=  y(2)

                                                                 c     y(2)' := -y(1)
                                                                 c
                                                                 c     Called by ODEPACK routine LSODA.
                                                                 c===========================================================
                                                                       subroutine fcn(neq,t,y,yprime)
                                                                          implicit    none
                                                                          integer     neq
                                                                          real*8      t,     y(neq),    yprime(neq)
                                                                          yprime(1) =  y(2)
                                                                          yprime(2) = -y(1)
                                                                          return
                                                                       end
                                                                 c===========================================================
                                                                 c     Implements Jacobian (optional).  Dummy routine in
                                                                 c     this case.
                                                                 c===========================================================
                                                                       subroutine jac
                                                                          implicit    none
                                                                          return
                                                                       end

 Source file: lsodafrag.f

      subroutine lsoda (f, neq, y, t, tout, itol, rtol, atol, itask,
     1            istate, iopt, rwork, lrw, iwork, liw, jac, jt)
      external f, jac
      integer neq, itol, itask, istate, iopt, lrw, iwork, liw, jt
      double precision y, t, tout, rtol, atol, rwork
      dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw)
c-----------------------------------------------------------------------
c this is the march 30, 1987 version of
c lsoda.. livermore solver for ordinary differential equations, with
c         automatic method switching for stiff and nonstiff problems.
c
c this version is in double precision.
c
c lsoda solves the initial value problem for stiff or nonstiff
c systems of first order ode-s,
c     dy/dt = f(t,y) ,  or, in component form,
c     dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq).
c
c this a variant version of the lsode package.
c it switches automatically between stiff and nonstiff methods.
c this means that the user does not have to determine whether the
c problem is stiff or not, and the solver will automatically choose the
c appropriate method.  it always starts with the nonstiff method.
c
c authors..
c                linda r. petzold  and  alan c. hindmarsh,
c                computing and mathematics research division, l-316
c                lawrence livermore national laboratory
c                livermore, ca 94550.
c
c references..
c 1.  alan c. hindmarsh,  odepack, a systematized collection of ode
c     solvers, in scientific computing, r. s. stepleman et al. (eds.),
c     north-holland, amsterdam, 1983, pp. 55-64.
c 2.  linda r. petzold, automatic selection of methods for solving
c     stiff and nonstiff systems of ordinary differential equations,
c     siam j. sci. stat. comput. 4 (1983), pp. 136-148.
c-----------------------------------------------------------------------
c summary of usage.
c
                               .
                               .
                               .

Source file: Makefile

                                                                 .IGNORE:
                                                                 F77_COMPILE  = $(F77) $(F77FLAGS) $(F77CFLAGS)
                                                                 F77_LOAD     = $(F77) $(F77FLAGS) $(F77LFLAGS)
                                                                 .f.o:
                                                                    $(F77_COMPILE) $*.f
                                                                 EXECUTABLES = tlsoda
                                                                 all: $(EXECUTABLES)
                                                                 # Note that usage of 'odepack' library also requires linking to
                                                                 # 'linpack' library (LINPACK is an antecedent of LAPACK)
                                                                 tlsoda: tlsoda.o dvfrom.o
                                                                    $(F77_LOAD) tlsoda.o dvfrom.o -lp410f \
                                                                                -lodepack -llinpack $(LIBBLAS) -o tlsoda
                                                                 clean:
                                                                    rm *.o
                                                                    rm $(EXECUTABLES)

 Source file: sgi-output

#############################################################
# Building 'tlsoda' and sample output on sgi1.
#############################################################
sgi1% pwd; ls
/usr/people/phys410/ode/ex1
Makefile   dvfrom.f   tlsoda.f
#############################################################
#  NOTE: 'lsoda' uses routines from a predecessor of LAPACK
#  called LINPACK as well as BLAS routines.  All three of the
#  libraries: 'odepack', 'linpack', and 'blas' must be
#  specified at load time.
#############################################################
sgi1% make
f77 -g -64 -c tlsoda.f
f77 -g -64 -c dvfrom.f
f77 -g -64 -L/usr/local/lib tlsoda.o dvfrom.o -lp410f \
               -lodepack -llinpack -lblas -o tlsoda
sgi1% tlsoda
 usage: tlsoda <y_i> <dy/dx_i>
#############################################################
# Generate a uniform-mesh of x output-values using 'dvmesh'
# and pipe it into tlsoda with initial conditions
# u(0) = 0, u'(0) = 1  --> u(x) = sin(x)
#############################################################
sgi1% dvmesh 0.0 1.0 11 | tlsoda 0.0 1.0
   0.0000000000000000E+00  0.0000000000000000E+00
   0.1000000000000000      9.9833904829367021E-02
   0.2000000000000000      0.1986703509344916
   0.3000000000000000      0.2955222777237400
   0.4000000000000000      0.3894213005717221
   0.5000000000000000      0.4794291455964967
   0.6000000000000000      0.5646464181480721
   0.7000000000000000      0.6442215920888248
   0.7999999999999999      0.7173595216370443
   0.9000000000000000      0.7833295096871006
    1.000000000000000      0.8414729721292930
#############################################################
# Use 'nf' to compute error in solution.  Note that, at least
# at the initial time, the error in the solution is roughly
# within the specified tolerance of 1.0d-6.
#############################################################
sgi1% dvmesh 0.0 1.0 11 | tlsoda 0.0 1.0 | nf _1 'sin(_1) - _2'
0.0000000000000000E+00   0
0.1000000000000000   -4.88182538865845e-07
0.2000000000000000   -1.02013943037438e-06
0.3000000000000000   -2.07106240046606e-06
0.4000000000000000   -2.95826307156188e-06
0.5000000000000000   -3.606992293681e-06
0.6000000000000000   -3.94475303677755e-06
0.7000000000000000   -3.90485113377359e-06
0.7999999999999999   -3.43073752162582e-06
0.9000000000000000   -2.60005961716381e-06
1.000000000000000   -1.98732139644608e-06