*DECK DHSTRT SUBROUTINE DHSTRT (DF, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL, + BIG, SPY, PV, YP, SF, RPAR, IPAR, H) C***BEGIN PROLOGUE DHSTRT C***SUBSIDIARY C***PURPOSE Subsidiary to DDEABM, DDEBDF and DDERKF C***LIBRARY SLATEC C***TYPE DOUBLE PRECISION (HSTART-S, DHSTRT-D) C***AUTHOR Watts, H. A., (SNLA) C***DESCRIPTION C C DHSTRT computes a starting step size to be used in solving initial C value problems in ordinary differential equations. C C ********************************************************************** C ABSTRACT C C Subroutine DHSTRT computes a starting step size to be used by an C initial value method in solving ordinary differential equations. C It is based on an estimate of the local Lipschitz constant for the C differential equation (lower bound on a norm of the Jacobian) , C a bound on the differential equation (first derivative) , and C a bound on the partial derivative of the equation with respect to C the independent variable. C (all approximated near the initial point A) C C Subroutine DHSTRT uses a function subprogram DHVNRM for computing C a vector norm. The maximum norm is presently utilized though it C can easily be replaced by any other vector norm. It is presumed C that any replacement norm routine would be carefully coded to C prevent unnecessary underflows or overflows from occurring, and C also, would not alter the vector or number of components. C C ********************************************************************** C On input you must provide the following C C DF -- This is a subroutine of the form C DF(X,U,UPRIME,RPAR,IPAR) C which defines the system of first order differential C equations to be solved. For the given values of X and the C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must C evaluate the NEQ components of the system of differential C equations DU/DX=DF(X,U) and store the derivatives in the C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for C equations I=1,...,NEQ. C C Subroutine DF must not alter X or U(*). You must declare C the name DF in an external statement in your program that C calls DHSTRT. You must dimension U and UPRIME in DF. C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C program and subroutine DF. They are not used or altered by C DHSTRT. If you do not need RPAR or IPAR, ignore these C parameters by treating them as dummy arguments. If you do C choose to use them, dimension them in your program and in C DF as arrays of appropriate length. C C NEQ -- This is the number of (first order) differential equations C to be integrated. C C A -- This is the initial point of integration. C C B -- This is a value of the independent variable used to define C the direction of integration. A reasonable choice is to C set B to the first point at which a solution is desired. C You can also use B, if necessary, to restrict the length C of the first integration step because the algorithm will C not compute a starting step length which is bigger than C ABS(B-A), unless B has been chosen too close to A. C (it is presumed that DHSTRT has been called with B C different from A on the machine being used. Also see the C discussion about the parameter SMALL.) C C Y(*) -- This is the vector of initial values of the NEQ solution C components at the initial point A. C C YPRIME(*) -- This is the vector of derivatives of the NEQ C solution components at the initial point A. C (defined by the differential equations in subroutine DF) C C ETOL -- This is the vector of error tolerances corresponding to C the NEQ solution components. It is assumed that all C elements are positive. Following the first integration C step, the tolerances are expected to be used by the C integrator in an error test which roughly requires that C ABS(LOCAL ERROR) .LE. ETOL C for each vector component. C C MORDER -- This is the order of the formula which will be used by C the initial value method for taking the first integration C step. C C SMALL -- This is a small positive machine dependent constant C which is used for protecting against computations with C numbers which are too small relative to the precision of C floating point arithmetic. SMALL should be set to C (approximately) the smallest positive DOUBLE PRECISION C number such that (1.+SMALL) .GT. 1. on the machine being C used. The quantity SMALL**(3/8) is used in computing C increments of variables for approximating derivatives by C differences. Also the algorithm will not compute a C starting step length which is smaller than C 100*SMALL*ABS(A). C C BIG -- This is a large positive machine dependent constant which C is used for preventing machine overflows. A reasonable C choice is to set big to (approximately) the square root of C the largest DOUBLE PRECISION number which can be held in C the machine. C C SPY(*),PV(*),YP(*),SF(*) -- These are DOUBLE PRECISION work C arrays of length NEQ which provide the routine with needed C storage space. C C RPAR,IPAR -- These are parameter arrays, of DOUBLE PRECISION and C INTEGER type, respectively, which can be used for C communication between your program and the DF subroutine. C They are not used or altered by DHSTRT. C C ********************************************************************** C On Output (after the return from DHSTRT), C C H -- is an appropriate starting step size to be attempted by the C differential equation method. C C All parameters in the call list remain unchanged except for C the working arrays SPY(*),PV(*),YP(*), and SF(*). C C ********************************************************************** C C***SEE ALSO DDEABM, DDEBDF, DDERKF C***ROUTINES CALLED DHVNRM C***REVISION HISTORY (YYMMDD) C 820301 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 891024 Changed references from DVNORM to DHVNRM. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) C 910722 Updated AUTHOR section. (ALS) C***END PROLOGUE DHSTRT