c===========================================================
c tlsoda: Demo program which uses ODEPACK routine LSODA
c to solve the second-order ODE
c
c u''(t) = -u(t), 0 <= t <= tmax
c
c (' = d/dt), with initial conditions
c
c u(0) = u0, u'(0) = du0
c
c The exact solution is
c
c u_xct(t) = du0 * sin(t) + u0 * 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 = 0
c t_ntout = tmax
c ntout = 2**olevel + 1
c
c Output to standard error is the RMS error of the
c approximate solution.
c===========================================================
program tlsoda
implicit none
integer iargc, i4arg
real*8 r8arg
c-----------------------------------------------------------
c Command-line arguments:
c
c tmax: Final integration time
c u0: Initial value: u(0)
c du0: Initial value: u'(0)
c tol: Error tolerance (this program uses LSODA's
c pure absolute error control)
c olevel: Output level:
c dtout = tmax/(2**olevel+1)
c-----------------------------------------------------------
real*8 tmax, u0, du0, 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; if and when needed, LSODA will compute
c a finite-difference approximation to the Jacobian.
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; we use jt = 2 here which, as mentioned
c above, instructs LSODA to compute a finite-difference
c approximation to the Jacobian if and when needed.
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. 5 ) go to 900
tmax = r8arg(1,r8_never)
u0 = r8arg(2,r8_never)
du0 = r8arg(3,r8_never)
tol = r8arg(4,r8_never)
olevel = i4arg(5,-1)
if( tmax .eq. r8_never .or. u0 .eq. r8_never .or.
& du0 .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 / (ntout - 1)
t = 0.0d0
y(1) = u0
y(2) = du0
c-----------------------------------------------------------
c Output initial solution and error and initialize
c rms error.
c-----------------------------------------------------------
err = du0 * sin(t) + u0 * 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 replaces 't' with the value
c of 'tout' if the 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 ',1p,e11.3,0p,
& ' .. ',1p,e11.3,0p/)
go to 500
end if
c-----------------------------------------------------------
c Output the solution and error, and update RMS error
c accumulator.
c-----------------------------------------------------------
err = du0 * sin(t) + u0 * cos(t) - y(1)
write(*,*) t, y(1), err
rmserr = rmserr + err**2
end do
c-----------------------------------------------------------
c Output the RMS error to standard error.
c-----------------------------------------------------------
rmserr = sqrt(rmserr / ntout)
write(0,*) 'rmserr: ', rmserr
500 continue
stop
900 continue
write(0,*) 'usage: tlsoda <tmax> <u0> <du0> '//
& '<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