STEPS(3)

SYNOPSIS
SUBROUTINE STEPS (F, NEQN, Y, X, H, EPS, WT, START, HOLD, K, KOLD, + CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G, PHASE1, NS, + NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV, KGI, GI, + RPAR, IPAR)

NAME
steps - Integrate a system of first order ordinary differential equations one step.

LIBRARY
SLATEC (DEPAC)

CATEGORY
I1A1B

TYPE
SINGLE PRECISION (STEPS-S, DSTEPS-D)

KEYWORDS
ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE, ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR

AUTHOR
Shampine, L. F., (SNLA) Gordon, M. K., (SNLA) MODIFIED BY H.A. WATTS

DESCRIPTION
Written by L. F. Shampine and M. K. Gordon Abstract Subroutine STEPS is normally used indirectly through subroutine DEABM . Because DEABM suffices for most problems and is much easier to use, using it should be considered before using STEPS alone. Subroutine STEPS integrates a system of NEQN first order ordinary differential equations one step, normally from X to X+H, using a modified divided difference form of the Adams Pece formulas. Local extrapolation is used to improve absolute stability and accuracy. The code adjusts its order and step size to control the local error per unit step in a generalized sense. Special devices are included to control roundoff error and to detect when the user is requesting too much accuracy. This code is completely explained and documented in the text, Computer Solution of Ordinary Differential Equations, The Initial Value Problem by L. F. Shampine and M. K. Gordon. Further details on use of this code are available in "Solving Ordinary Differential Equations with ODE, STEP, and INTRP", by L. F. Shampine and M. K. Gordon, SLA-73-1060. The parameters represent -- F -- subroutine to evaluate derivatives NEQN -- number of equations to be integrated Y(*) -- solution vector at X X -- independent variable H -- appropriate step size for next step. Normally determined by code EPS -- local error tolerance WT(*) -- vector of weights for error criterion START -- logical variable set .TRUE. for first step, .FALSE. otherwise HOLD -- step size used for last successful step K -- appropriate order for next step (determined by code) KOLD -- order used for last successful step CRASH -- logical variable set .TRUE. when no step can be taken, YP(*) -- derivative of solution vector at X after successful step KSTEPS -- counter on attempted steps TWOU -- 2.*U where U is machine unit roundoff quantity FOURU -- 4.*U where U is machine unit roundoff quantity RPAR,IPAR -- parameter arrays which you may choose to use for communication between your program and subroutine F. They are not altered or used by STEPS. The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G, W,P,IV and GI are required for the interpolation subroutine SINTRP. The remaining variables and arrays are included in the call list only to eliminate local retention of variables between calls. Input to STEPS First call -- The user must provide storage in his calling program for all arrays in the call list, namely DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12), 1 ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10), 2 RPAR(*),IPAR(*) **Note** The user must also declare START , CRASH , PHASE1 and NORND logical variables and F an EXTERNAL subroutine, supply the subroutine F(X,Y,YP) to evaluate DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN)) and initialize only the following parameters. NEQN -- number of equations to be integrated Y(*) -- vector of initial values of dependent variables X -- initial value of the independent variable H -- nominal step size indicating direction of integration and maximum size of step. Must be variable EPS -- local error tolerance per step. Must be variable WT(*) -- vector of non-zero weights for error criterion START -- .TRUE. YP(*) -- vector of initial derivative values KSTEPS -- set KSTEPS to zero TWOU -- 2.*U where U is machine unit roundoff quantity FOURU -- 4.*U where U is machine unit roundoff quantity Define U to be the machine unit roundoff quantity by calling the function routine R1MACH, U = R1MACH(4), or by computing U so that U is the smallest positive number such that 1.0+U .GT. 1.0. STEPS requires that the L2 norm of the vector with components LOCAL ERROR(L)/WT(L) be less than EPS for a successful step. The array WT allows the user to specify an error test appropriate for his problem. For example, WT(L) = 1.0 specifies absolute error, = ABS(Y(L)) error relative to the most recent value of the L-th component of the solution, = ABS(YP(L)) error relative to the most recent value of the L-th component of the derivative, = MAX(WT(L),ABS(Y(L))) error relative to the largest magnitude of L-th component obtained so far, = ABS(Y(L))*RELERR/EPS + ABSERR/EPS specifies a mixed relative-absolute test where RELERR is relative error, ABSERR is absolute error and EPS = MAX(RELERR,ABSERR) . Subsequent calls -- Subroutine STEPS is designed so that all information needed to continue the integration, including the step size H and the order K , is returned with each step. With the exception of the step size, the error tolerance, and the weights, none of the parameters should be altered. The array WT must be updated after each step to maintain relative error tests like those above. Normally the integration is continued just beyond the desired endpoint and the solution interpolated there with subroutine SINTRP . If it is impossible to integrate beyond the endpoint, the step size may be reduced to hit the endpoint since the code will not take a step larger than the H input. Changing the direction of integration, i.e., the sign of H , requires the user set START = .TRUE. before calling STEPS again. This is the only situation in which START should be altered. Output from STEPS Successful Step -- The subroutine returns after each successful step with START and CRASH set .FALSE. . X represents the independent variable advanced one step of length HOLD from its value on input and Y the solution vector at the new value of X . All other parameters represent information corresponding to the new X needed to continue the integration. Unsuccessful Step -- When the error tolerance is too small for the machine precision, the subroutine returns without taking a step and CRASH = .TRUE. . An appropriate step size and error tolerance for continuing are estimated and all other information is restored as upon input before returning. To continue with the larger tolerance, the user just calls the code again. A restart is neither required nor desirable.

REFERENCES
L. F. Shampine and M. K. Gordon, Solving ordinary differential equations with ODE, STEP, and INTRP, Report SLA-73-1060, Sandia Laboratories, 1973.

ROUTINES CALLED
HSTART, R1MACH

REVISION HISTORY (YYMMDD)
740101 DATE WRITTEN
890531 Changed all specific intrinsics to generic. (WRB)
890831 Modified array declarations. (WRB)
890831 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
920501 Reformatted the REFERENCES section. (WRB)