*DECK HSTART SUBROUTINE HSTART (F, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL, + BIG, SPY, PV, YP, SF, RPAR, IPAR, H) C***BEGIN PROLOGUE HSTART C***SUBSIDIARY C***PURPOSE Subsidiary to DEABM, DEBDF and DERKF C***LIBRARY SLATEC C***TYPE SINGLE PRECISION (HSTART-S, DHSTRT-D) C***AUTHOR Watts, H. A., (SNLA) C***DESCRIPTION C C HSTART computes a starting step size to be used in solving initial C value problems in ordinary differential equations. C ********************************************************************** C Abstract C C Subroutine HSTART 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 HSTART uses a function subprogram HVNRM 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 F -- This is a subroutine of the form C F(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=F(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 F must not alter X or U(*). You must declare C the name F in an EXTERNAL statement in your program that C calls HSTART. You must dimension U and UPRIME in F. C C RPAR and IPAR are real and integer parameter arrays which C you can use for communication between your program and C subroutine F. They are not used or altered by HSTART. If C you do not need RPAR or IPAR, ignore these parameters by C treating them as dummy arguments. If you do choose to use C them, dimension them in your program and in F as arrays C 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 HSTART has been called with B C different from A on the machine being used. Also see C the 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 F) 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 real number such C that (1.+SMALL) .GT. 1. on the machine being used. the C quantity SMALL**(3/8) is used in computing increments of C variables for approximating derivatives by differences. C also the algorithm will not compute a starting step length C which is smaller than 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 real number which can be held in the machine. C C SPY(*),PV(*),YP(*),SF(*) -- These are real work arrays of length C NEQ which provide the routine with needed storage space. C C RPAR,IPAR -- These are parameter arrays, of real and integer C type, respectively, which can be used for communication C between your program and the F subroutine. They are not C used or altered by HSTART. C C ********************************************************************** C On Output (after the return from HSTART), 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 DEABM, DEBDF, DERKF C***ROUTINES CALLED HVNRM C***REVISION HISTORY (YYMMDD) C 800501 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891024 Changed references from VNORM to HVNRM. (WRB) C 891024 REVISION DATE from Version 3.2 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 HSTART C DIMENSION Y(*),YPRIME(*),ETOL(*),SPY(*),PV(*),YP(*),SF(*), 1 RPAR(*),IPAR(*) EXTERNAL F C C....................................................................... C C***FIRST EXECUTABLE STATEMENT HSTART DX = B - A ABSDX = ABS(DX) RELPER = SMALL**0.375 YNORM = HVNRM(Y,NEQ) C C....................................................................... C C COMPUTE A WEIGHTED APPROXIMATE BOUND (DFDXB) ON THE PARTIAL C DERIVATIVE OF THE EQUATION WITH RESPECT TO THE C INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW. ALSO C COMPUTE A WEIGHTED BOUND (FBND) ON THE FIRST DERIVATIVE LOCALLY. C DA = SIGN(MAX(MIN(RELPER*ABS(A),ABSDX),100.*SMALL*ABS(A)),DX) IF (DA .EQ. 0.) DA = RELPER*DX CALL F(A+DA,Y,SF,RPAR,IPAR) C IF (MORDER .EQ. 1) GO TO 20 POWER = 2./(MORDER+1) DO 10 J=1,NEQ WTJ = ETOL(J)**POWER SPY(J) = SF(J)/WTJ YP(J) = YPRIME(J)/WTJ 10 PV(J) = SPY(J) - YP(J) GO TO 40 C 20 DO 30 J=1,NEQ SPY(J) = SF(J)/ETOL(J) YP(J) = YPRIME(J)/ETOL(J) 30 PV(J) = SPY(J) - YP(J) C 40 DELF = HVNRM(PV,NEQ) DFDXB = BIG IF (DELF .LT. BIG*ABS(DA)) DFDXB = DELF/ABS(DA) YPNORM = HVNRM(YP,NEQ) FBND = MAX(HVNRM(SPY,NEQ),YPNORM) C C....................................................................... C C COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ CONSTANT FOR C THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS ALSO REPRESENTS AN C ESTIMATE OF THE NORM OF THE JACOBIAN LOCALLY. C THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO ESTIMATE THE C LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES. THE FIRST C PERTURBATION VECTOR IS BASED ON THE INITIAL DERIVATIVES AND C DIRECTION OF INTEGRATION. THE SECOND PERTURBATION VECTOR IS C FORMED USING ANOTHER EVALUATION OF THE DIFFERENTIAL EQUATION. C THE THIRD PERTURBATION VECTOR IS FORMED USING PERTURBATIONS BASED C ONLY ON THE INITIAL VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS C CHANGED TO NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN C INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT COMPONENTS C OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE CONSISTENT WITH C THE SLOPES OF LOCAL SOLUTION CURVES. C ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST DERIVATIVE. C NO ATTEMPT IS MADE TO KEEP THE PERTURBATION VECTOR SIZE CONSTANT. C IF (YPNORM .EQ. 0.) GO TO 60 C USE INITIAL DERIVATIVES FOR FIRST PERTURBATION ICASE = 1 DO 50 J=1,NEQ SPY(J) = YPRIME(J) 50 YP(J) = YPRIME(J) GO TO 80 C CANNOT HAVE A NULL PERTURBATION VECTOR 60 ICASE = 2 DO 70 J=1,NEQ SPY(J) = YPRIME(J) 70 YP(J) = ETOL(J) C 80 DFDUB = 0. LK = MIN(NEQ+1,3) DO 260 K=1,LK C SET YPNORM AND DELX YPNORM = HVNRM(YP,NEQ) IF (ICASE .EQ. 1 .OR. ICASE .EQ. 3) GO TO 90 DELX = SIGN(1.0,DX) GO TO 120 C TRY TO ENFORCE MEANINGFUL PERTURBATION VALUES 90 DELX = DX IF (ABS(DELX)*YPNORM .GE. RELPER*YNORM) GO TO 100 DELXB = BIG IF (RELPER*YNORM .LT. BIG*YPNORM) DELXB = RELPER*YNORM/YPNORM DELX = SIGN(DELXB,DX) 100 DO 110 J=1,NEQ IF (ABS(DELX*YP(J)) .GT. ETOL(J)) DELX=SIGN(ETOL(J)/YP(J),DX) 110 CONTINUE C DEFINE PERTURBED VECTOR OF INITIAL VALUES 120 DO 130 J=1,NEQ 130 PV(J) = Y(J) + DELX*YP(J) IF (K .EQ. 2) GO TO 150 C EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED C VECTOR AND COMPUTE CORRESPONDING DIFFERENCES CALL F(A,PV,YP,RPAR,IPAR) DO 140 J=1,NEQ 140 PV(J) = YP(J) - YPRIME(J) GO TO 170 C USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE C IN COMPUTING ONE ESTIMATE 150 CALL F(A+DA,PV,YP,RPAR,IPAR) DO 160 J=1,NEQ 160 PV(J) = YP(J) - SF(J) C CHOOSE LARGEST BOUND ON THE WEIGHTED FIRST C DERIVATIVE 170 IF (MORDER .EQ. 1) GO TO 190 DO 180 J=1,NEQ 180 YP(J) = YP(J)/ETOL(J)**POWER GO TO 210 190 DO 200 J=1,NEQ 200 YP(J) = YP(J)/ETOL(J) 210 FBND = MAX(FBND,HVNRM(YP,NEQ)) C COMPUTE BOUND ON A LOCAL LIPSCHITZ CONSTANT DELF = HVNRM(PV,NEQ) IF (DELF .EQ. 0.) GO TO 220 DELY = ABS(DELX)*YPNORM IF (DELF .GE. BIG*DELY) GO TO 270 DFDUB = MAX(DFDUB,DELF/DELY) C 220 IF (K .EQ. LK) GO TO 280 C CHOOSE NEXT PERTURBATION VECTOR DO 250 J=1,NEQ IF (K .EQ. LK-1) GO TO 230 ICASE = 3 DY = ABS(PV(J)) IF (DY .EQ. 0.) DY = MAX(DELF,ETOL(J)) GO TO 240 230 ICASE = 4 DY = MAX(RELPER*ABS(Y(J)),ETOL(J)) 240 IF (SPY(J) .EQ. 0.) SPY(J) = YP(J) IF (SPY(J) .NE. 0.) DY = SIGN(DY,SPY(J)) 250 YP(J) = DY 260 CONTINUE C C PROTECT AGAINST AN OVERFLOW 270 DFDUB = BIG C C....................................................................... C C COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE C 280 YDPB = DFDXB + DFDUB*FBND C C....................................................................... C C COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND SECOND C DERIVATIVE INFORMATION C C RESTRICT THE STEP LENGTH TO BE NOT BIGGER THAN C ABS(B-A). (UNLESS B IS TOO CLOSE TO A) H = ABSDX C IF (YDPB .NE. 0. .OR. FBND .NE. 0.) GO TO 290 C C BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND C DERIVATIVE TERM (YDPB) ARE ZERO GO TO 310 C 290 IF (YDPB .NE. 0.) GO TO 300 C C ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO IF (1.0 .LT. FBND*ABSDX) H = 1./FBND GO TO 310 C C SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO 300 SRYDPB = SQRT(0.5*YDPB) IF (1.0 .LT. SRYDPB*ABSDX) H = 1./SRYDPB C C FURTHER RESTRICT THE STEP LENGTH TO BE NOT C BIGGER THAN 1/DFDUB 310 IF (H*DFDUB .GT. 1.) H = 1./DFDUB C C FINALLY, RESTRICT THE STEP LENGTH TO BE NOT C SMALLER THAN 100*SMALL*ABS(A). HOWEVER, IF C A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO, C THE ALGORITHM RETURNS SMALL*ABS(B) FOR THE C STEP LENGTH. H = MAX(H,100.*SMALL*ABS(A)) IF (H .EQ. 0.) H = SMALL*ABS(B) C C NOW SET DIRECTION OF INTEGRATION H = SIGN(H,DX) C RETURN END