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----------------------------------------------------------- c usage: tlsoda c----------------------------------------------------------- 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: dtout = tmax/2**olevel 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 c Set istate to 1 to indicate initial call, istate c should be set to 2 for subsequent calls, but lsoda c will automatically do this so long as the initial c call is successful. c----------------------------------------------------------- istate = 1 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 Note that LSODA replaces 't' with the value c of 'tout' if the integration is successful. c----------------------------------------------------------- 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 '// & ' ' 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