subroutine deroot(f,neqn,y,t,tout,relerr,abserr,iflag,
* g,reroot,aeroot)
external f
integer neqn, iflag
real y(neqn), t, tout, relerr, abserr, g,
* reroot, aeroot
c
c subroutine deroot integrates a system of up to 20 first order
c ordinary differential equations of the form
c dy(i)/dt = f(t,y(1),...,y(neqn))
c y(i) given at t.
c the subroutine integrates from t in the direction of tout until
c it locates the first root of the nonlinear equation
c g(t,y(1),...,y(neqn),yp(1),...,yp(neqn)) = 0.
c upon finding the root, the code returns with all parameters in the
c call list set for continuing the integration to the next root or
c the first root of a new function g . if no root is found, the
c integration proceeds to tout . again all parameters are set to
c continue.
c
c the differential equations are actually solved by a suite of codes,
c deroot , step , and intrp . deroot is a supervisor which directs
c the integration. it calls on step to advance the solution and
c intrp to interpolate the solution and its derivative. step uses
c a modified divided difference form of the adams pece formulas and
c local extrapolation. it adjusts the order and step size to control
c the local error per unit step in a generalized sense. normally each
c call to step advances the solution one step in the direction of
c tout . for reasons of efficiency deroot integrates beyond tout
c internally, though never beyond t+10*(tout-t), and calls intrp to
c interpolate the solution and derivative at tout . an option is
c provided to stop the integration at tout but it should be used
c only if it is impossible to continue the integration beyond tout .
c
c after each internal step, deroot evaluates the function g and
c checks for a change in sign in the function value from the
c preceding step. such a change indicates a root lies in the
c interval of the step just completed. deroot then calls subroutine
c root to reduce the bracketing interval until the root is
c determined to the desired accuracy. subroutine root uses a
c combination of the secant rule and bisection to do this. the
c solution and derivative values required are obtained by
c interpolation with intrp . the code locates only those roots
c for which g changes sign in (t,tout) and for which a
c bracketing interval exists. in particular, it does not locate
c a root at the initial point t .
c
c the codes step and intrp and that portion of deroot which
c directs the integration are completely explained and documented in
c the text, computer solution of ordinary differential equations,
c the initial value problem by l. f. shampine and m. k. gordon.
c subroutine root is a slightly modified version of the root-solver
c discussed in the text, numerical computing, an introduction by
c l. f. shampine and r. c. allen.
c
c the parameters for deroot are
c f -- subroutine f(t,y,yp) to evaluate derivatives yp(i)=dy(i)/dt
c neqn -- number of equations to be integrated
c y(*) -- solution vector at t
c t -- independent variable
c tout -- arbitrary point beyond the root desired
c relerr,abserr -- relative and absolute error tolerances for local
c error test. at each step the code requires
c abs(local error) .le. abs(y)*relerr + abserr
c for each component of the local error and solution vectors
c iflag -- indicates status of integration
c g - function of t, y(*), yp(*) whose root is desired.
c reroot, aeroot -- relative and absolute error tolerances for
c accepting the root. the interval containing the root is
c reduced until it satisfies
c 0.5*abs(length of interval) .le. reroot*abs(root)+aeroot
c where root is that endpoint yielding the smaller value of
c g in magnitude. pure relative error is not recommended
c if the root might be zero.
c
c first call to deroot --
c
c the user must provide storage in his calling program for the
c array in the call list,
c y(neqn)
c and declare f , g in an external statement. he must supply the
c subroutine f(t,y,yp) to evaluate
c dy(i)/dt = yp(i) = f(t,y(1),...,y(neqn))
c and the function g(t,y,yp) to evaluate
c g = g(t,y(1),...,y(neqn),yp(1),...,yp(neqn)).
c note that the array yp is an input argument and should not be
c computed in the function subprogram. finally the user must
c initialize the parameters
c neqn -- number of equations to be integrated
c y(*) -- vector of initial conditions
c t -- starting point of integration
c tout -- arbitrary point beyond the root desired
c relerr,abserr -- relative and absolute local error tolerances
c for integrating the equations
c iflag -- +1,-1. indicator to initialize the code. normal input
c is +1. the user should set iflag=-1 only if it is
c impossible to continue the integration beyond tout .
c reroot,aeroot -- relative and absolute error tolerances for
c computing the root of g
c
c all parameters except f, g, neqn, tout, reroot and aeroot may be
c altered by the code on output so must be variables in the calling
c program.
c
c output from deroot --
c
c neqn -- unchanged
c y(*) -- solution at t
c t -- last point reached in integration. normal return has
c t = tout or t = root
c tout -- unchanged
c relerr,abserr -- normal return has tolerances unchanged. iflag=3
c signals tolerances increased
c iflag = 2 -- normal return. integration reached tout
c = 3 -- integration did not reach tout because error
c tolerances too small. relerr , abserr increased
c appropriately for continuing
c = 4 -- integration did not reach tout because more than
c maxnum steps needed
c = 5 -- integration did not reach tout because equations
c appear to be stiff
c = 6 -- invalid input parameters (fatal error)
c = 7 -- normal return. a root was found which satisfied
c the error criterion or had a zero residual
c = 8 -- abnormal return. an odd order pole of g was
c found.
c = 9 -- abnormal return. too many evaluations of g were
c required (as programmed 500 are allowed.)
c the value of iflag is returned negative when the input
c value is negative and the integration does not reach
c tout , i.e., -3,-4,-5,-7,-8,-9.
c reroot,aeroot -- unchanged
c
c subsequent calls to deroot --
c
c subroutine deroot returns with all information needed to continue
c the integration. if the integration did not reach tout and the
c user wants to continue, he just calls again. if the integration
c reached tout , the user need only define a new tout and call
c again. the output value of iflag is the appropriate input value
c for subsequent calls. the only situation in which it should be
c altered is to stop the integration internally at the new tout ,
c i.e., change output iflag=2 to input iflag=-2 . only the error
c tolerances and the function g may be changed by the user before
c continuing. all other parameters must remain unchanged.
c
logical stiff,crash,start
real fouru, eps, gxold, gx, x, delsgn,
* b, c, gc, del, tend, releps, abseps, troot,
* absdel, h, hold, told, gt
real yy(20),wt(20),phi(20,16),p(20),yp(20),ypout(20),
* psi(12)
real abs, amax1, amin1, r1mach, sign ***
c
c***********************************************************************
c* the only machine dependent constant is based on the machine unit *
c* roundoff error u which is the smallest positive number such that *
c* 1.0+u .gt. 1.0 . u must be calculated and fouru=4.0*u inserted *
c* in the following statement before using deroot . the subroutine *
c* machin calculates u . fouru and twou=2.0*u must also be *
c* inserted in subroutine step before calling deroot . *
c data fouru/8.8d-16/ ***
c***********************************************************************
c
c the constant maxnum is the maximum number of steps allowed in one
c call to deroot . the user may change this limit by altering the
c following statement
data maxnum/500/
c
c this version of deroot allows a maximum of 20 equations. to
c increase this number, only the number 20 in the dimension statement
c and in the following statement need be changed
c *** *** ***
c test for improper parameters
c
c-----------------------------------------------------------------
fouru = 4.0*r1mach(4) ***
if(neqn .lt. 1 .or. neqn .gt. 20) go to 10
if(t .eq. tout) go to 10
if(relerr .lt. 0.0 .or. abserr .lt. 0.0) go to 10
eps = amax1(relerr,abserr)
if(eps .le. 0.0) go to 10
if(reroot .lt. 0.0 .or. aeroot .lt. 0.0) go to 10
if(reroot+aeroot .le. 0.0) go to 10
if(iflag .eq. 0) go to 10
isn = isign(1,iflag)
iflag = iabs(iflag)
if(iflag .eq. 1) go to 20
if(t .ne. told) go to 10
if(iflag .ge. 2 .and. iflag .le. 5) go to 15
if(iflag .ge. 7 .and. iflag .le. 9) go to 15
10 iflag = 6
return
c
c for a new function g, check for a root in interval of step
c just completed
c
15 gxold = gx
gx = g(x,yy,yp)
if(gx .eq. gxold) go to 20
if(iflag .gt. 2 .and. iflag .le. 5) go to 20
if(isnold .lt. 0 .or. delsgn*(tout-t) .lt. 0.0) go to 20
jflag = 1
b = x
c = t
call intrp(x,yy,c,y,ypout,neqn,kold,phi,psi)
gc = g(c,y,ypout)
if(sign(1.0,gc)*sign(1.0,gx) .lt. 0.0) go to 140
if(gc .eq. 0.0 .or. gx .eq. 0.0) go to 140
c
c on each call set interval of integration and counter for number of
c steps. adjust input error tolerances to define weight vector for
c subroutine step
c
20 del = tout - t
absdel = abs(del)
tend = t + 10.0*del
if(isn .lt. 0) tend = tout
nostep = 0
kle4 = 0
stiff = .false.
releps = relerr/eps
abseps = abserr/eps
if(iflag .eq. 1) go to 30
if(isnold .lt. 0) go to 30
if(delsgn*del .gt. 0.0) go to 50
c
c on start and restart also set work variables x and yy(*), store the
c direction of integration, initialize the step size, and evaluate g
c
30 start = .true.
x = t
troot = t
do 40 l = 1,neqn
40 yy(l) = y(l)
delsgn = sign(1.0,del)
h = sign(amax1(abs(tout-x),fouru*abs(x)),tout-x)
call f(x,yy,yp)
gx = g(x,yy,yp)
c
c if already past output point, interpolate and return
c
50 gxold = gx
if(abs(x-t) .lt. absdel) go to 60
call intrp(x,yy,tout,y,ypout,neqn,kold,phi,psi)
iflag = 2
t = tout
told = t
isnold = isn
return
c
c if cannot go past output point and sufficiently close,
c extrapolate and return
c
60 if(isn .gt. 0 .or. abs(tout-x) .ge. fouru*abs(x)) go to 80
h = tout - x
call f(x,yy,yp)
do 70 l = 1,neqn
70 y(l) = yy(l) + h*yp(l)
iflag = 2
t = tout
told = t
isnold = isn
return
c
c test for too much work
c
80 if(nostep .lt. maxnum) go to 100
iflag = isn*4
if(stiff) iflag = isn*5
do 90 l = 1,neqn
90 y(l) = yy(l)
t = x
told = t
isnold = 1
return
c
c limit step size, set weight vector and take a step
c
100 h = sign(amin1(abs(h),abs(tend-x)),h)
do 110 l = 1,neqn
110 wt(l) = releps*abs(yy(l)) + abseps
call step(x,yy,f,neqn,h,eps,wt,start,
* hold,k,kold,crash,phi,p,yp,psi)
c
c test for tolerances too small. if so, set the derivative at x
c before returning
c
if(.not.crash) go to 130
iflag = isn*3
relerr = eps*releps
abserr = eps*abseps
do 120 l = 1,neqn
yp(l) = phi(l,1)
120 y(l) = yy(l)
t = x
told = t
isnold = 1
return
c
c augment counter on work and test for stiffness. also test for a
c root in the step just completed
c
130 nostep = nostep + 1
kle4 = kle4 + 1
if(kold .gt. 4) kle4 = 0
if(kle4 .ge. 50) stiff = .true.
gx = g(x,yy,yp)
if(sign(1.0,gxold)*sign(1.0,gx) .lt. 0.0) go to 135
if(gx .eq. 0.0 .or. gxold .eq. 0.0) go to 135
go to 50
c
c locate root of g. interpolate with intrp for solution and
c derivative values
c
135 b = x
c = x - hold
jflag = 1
140 call root(t,gt,b,c,reroot,aeroot,jflag)
if(jflag .gt. 0) go to 150
call intrp(x,yy,t,y,ypout,neqn,kold,phi,psi)
gt = g(t,y,ypout)
go to 140
150 continue
iflag = jflag + 6
if(jflag .eq. 2 .or. jflag .eq. 4) iflag = 7
if(jflag .eq. 3) iflag = 8
if(jflag .eq. 5) iflag = 9
iflag = iflag*isn
call intrp(x,yy,b,y,ypout,neqn,kold,phi,psi)
t = b
if(abs(t-troot) .le. reroot*abs(t) + aeroot) go to 50
troot = t
told = t
isnold = 1
return
end
subroutine intrp(x,y,xout,yout,ypout,neqn,kold,phi,psi)
integer neqn, kold
real x, y(neqn), xout, yout(neqn), ypout(neqn),
* phi(20,16), psi(12)
c
c the methods in subroutine step approximate the solution near x
c by a polynomial. subroutine intrp approximates the solution at
c xout by evaluating the polynomial there. information defining this
c polynomial is passed from step so intrp cannot be used alone.
c
c this code is completely explained and documented in the text,
c computer solution of ordinary differential equations the initial
c value problem by l. f. shampine and m. k. gordon.
c
c input to intrp --
c
c the user provides storage in the calling program for the arrays in
c the call list
c and defines
c xout -- point at which solution is desired.
c the remaining parameters are defined in step and passed to intrp
c from that subroutine
c
c output from intrp --
c
c yout(*) -- solution at xout
c ypout(*) -- derivative of solution at xout
c the remaining parameters are returned unaltered from their input
c values. integration with step may be continued.
c
real hi, temp1, term, psijm1, gamma, eta, temp2,
* temp3
real g(13),w(13),rho(13)
data g(1)/1.0/,rho(1)/1.0/
c
hi = xout - x
ki = kold + 1
kip1 = ki + 1
c
c initialize w(*) for computing g(*)
c
do 5 i = 1,ki
temp1 = i
5 w(i) = 1.0/temp1
term = 0.0
c
c compute g(*)
c
do 15 j = 2,ki
jm1 = j - 1
psijm1 = psi(jm1)
gamma = (hi + term)/psijm1
eta = hi/psijm1
limit1 = kip1 - j
do 10 i = 1,limit1
10 w(i) = gamma*w(i) - eta*w(i+1)
g(j) = w(1)
rho(j) = gamma*rho(jm1)
15 term = psijm1
c
c interpolate
c
do 20 l = 1,neqn
ypout(l) = 0.0
20 yout(l) = 0.0
do 30 j = 1,ki
i = kip1 - j
temp2 = g(i)
temp3 = rho(i)
do 25 l = 1,neqn
yout(l) = yout(l) + temp2*phi(l,i)
25 ypout(l) = ypout(l) + temp3*phi(l,i)
30 continue
do 35 l = 1,neqn
35 yout(l) = y(l) + hi*yout(l)
return
end
subroutine root(t,ft,b,c,relerr,abserr,iflag)
real t, ft, b, c, relerr, abserr
integer iflag
c
c root computes a root of the nonlinear equation f(x)=0
c where f(x) is a continuous real function of a single real
c variable x. the method used is a combination of bisection
c and the secant rule.
c
c normal input consists of a continuous function f and an
c interval (b,c) such that f(b)*f(c).le.0.0. each iteration
c finds new values of b and c such that the interval (b,c) is
c shrunk and f(b)*f(c).le.0.0. the stopping criterion is
c
c abs(b-c).le.2.0*(relerr*abs(b)+abserr)
c
c where relerr=relative error and abserr=absolute error are
c input quantities. set the flag, iflag, positive to initialize
c the computation. as b,c and iflag are used for both input and
c output, they must be variables in the calling program.
c
c if 0 is a possible root, one should not choose abserr=0.0.
c
c the output value of b is the better approximation to a root
c as b and c are always redefined so that abs(f(b)).le.abs(f(c)).
c
c to solve the equation, root must evaluate f(x) repeatedly. this
c is done in the calling program. when an evaluation of f is
c needed at t, root returns with iflag negative. evaluate ft=f(t)
c and call root again. do not alter iflag.
c
c when the computation is complete, root returns to the calling
c program with iflag positive:
c
c iflag=1 if f(b)*f(c).lt.0 and the stopping criterion is met.
c
c =2 if a value b is found such that the computed value
c f(b) is exactly zero. the interval (b,c) may not
c satisfy the stopping criterion.
c
c =3 if abs(f(b)) exceeds the input values abs(f(b)),
c abs(f(c)). in this case it is likely that b is close
c to a pole of f.
c
c =4 if no odd order root was found in the interval. a
c local minimum may have been obtained.
c
c =5 if too many function evaluations were made.
c (as programmed, 500 are allowed.)
c
c this code is a modification of the code zeroin which is completely
c explained and documented in the text, numerical computing: an
c introduction by l. f. shampine and r. c. allen.
c
real u, re, ae, acbs, a, fa, fb, fc, fx, cmb,
* acmb, tol, p, q
real abs, amax1, r1mach, sign
c***********************************************************************
c* the only machine dependent constant is based on the machine unit *
c* roundoff error u which is the smallest positive number such that *
c* 1.0+u .gt. 1.0 . u must be calculated and inserted in the *
c* following data statement before using root . the routine machin *
c* calculates u . *
c data u/2.2d-16/
c***********************************************************************
c
u = r1mach(4)
if(iflag.ge.0) go to 100
iflag=iabs(iflag)
go to (200,300,400), iflag
100 re=amax1(relerr,u)
ae=amax1(abserr,0.0)
ic=0
acbs=abs(b-c)
a=c
t=a
iflag=-1
return
200 fa=ft
t=b
iflag=-2
return
300 fb=ft
fc=fa
kount=2
fx=amax1(abs(fb),abs(fc))
1 if(abs(fc).ge.abs(fb))go to 2
c
c interchange b and c so that abs(f(b)).le.abs(f(c)).
c
a=b
fa=fb
b=c
fb=fc
c=a
fc=fa
2 cmb=0.5*(c-b)
acmb=abs(cmb)
tol=re*abs(b)+ae
c
c test stopping criterion and function count.
c
if(acmb.le.tol)go to 8
if(kount.ge.500)go to 12
c
c calculate new iterate implicitly as b+p/q
c where we arrange p.ge.0. the implicit
c form is used to prevent overflow.
c
p=(b-a)*fb
q=fa-fb
if(p.ge.0.0)go to 3
p=-p
q=-q
c
c update a, check if reduction in the size of bracketing
c interval is satisfactory. if not, bisect until it is.
c
3 a=b
fa=fb
ic=ic+1
if(ic.lt.4)go to 4
if(8.0*acmb.ge.acbs)go to 6
ic=0
acbs=acmb
c
c test for too small a change.
c
4 if(p.gt.abs(q)*tol)go to 5
c
c increment by tolerance.
c
b=b+sign(tol,cmb)
go to 7
c
c root ought to be between b and (c+b)/2.
c
5 if(p.ge.cmb*q)go to 6
c
c use secant rule.
c
b=b+p/q
go to 7
c
c use bisection.
c
6 b=0.5*(c+b)
c
c have completed computation for new iterate b.
c
7 t=b
iflag=-3
return
400 fb=ft
if(fb.eq.0.0)go to 9
kount=kount+1
if(sign(1.0,fb).ne.sign(1.0,fc))go to 1
c=a
fc=fa
go to 1
c
c finished. set iflag.
c
8 if(sign(1.0,fb).eq.sign(1.0,fc))go to 11
if(abs(fb).gt.fx)go to 10
iflag=1
return
9 iflag=2
return
10 iflag=3
return
11 iflag=4
return
12 iflag=5
return
end
subroutine step(x,y,f,neqn,h,eps,wt,start,
* hold,k,kold,crash,phi,p,yp,psi)
external f
integer neqn, k, kold
logical start, crash
real x, y(neqn), h, eps, wt(20), hold,
* phi(20,16), p(20), yp(20), psi(12)
c
c subroutine step integrates a system of first order ordinary
c differential equations one step, normally from x to x+h, using a
c modified divided difference form of the adams pece formulas. local
c extrapolation is used to improve absolute stability and accuracy.
c the code adjusts its order and step size to control the local error
c per unit step in a generalized sense. special devices are included
c to control roundoff error and to detect when the user is requesting
c too much accuracy.
c
c this code is completely explained and documented in the text,
c computer solution of ordinary differential equations the initial
c value problem by l. f. shampine and m. k. gordon.
c
c
c the parameters represent
c x -- independent variable
c y(*) -- solution vector at x
c yp(*) -- derivative of solution vector at x after successful
c step
c neqn -- number of equations to be integrated
c h -- appropriate step size for next step. normally determined by
c code
c eps -- local error tolerance. must be variable
c wt(*) -- vector of weights for error criterion
c start -- logical variable set .true. for first step, .false.
c otherwise
c hold -- step size used for last successful step
c k -- appropriate order for next step (determined by code)
c kold -- order used for last successful step
c crash -- logical variable set .true. when no step can be taken,
c .false. otherwise.
c the arrays phi, psi are required for the interpolation subroutine
c intrp . the array p is internal to the code.
c
c input to step
c
c first call --
c
c the user must provide storage in his driver program for all arrays
c in the call list, namely
c
c dimension y(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),psi(12)
c
c the user must also declare start and crash logical variables
c and f an external subroutine, supply the subroutine f(x,y,yp)
c to evaluate
c dy(i)/dx = yp(i) = f(x,y(1),y(2),...,y(neqn))
c and initialize only the following parameters
c x -- initial value of the independent variable
c y(*) -- vector of initial values of dependent variables
c neqn -- number of equations to be integrated
c h -- nominal step size indicating direction of integration
c and maximum size of step. must be variable
c eps -- local error tolerance per step. must be variable
c wt(*) -- vector of non-zero weights for error criterion
c start -- .true.
c
c step requires the l2 norm of the vector with components
c local error(l)/wt(l) be less than eps for a successful step. the
c array wt allows the user to specify an error test appropriate
c for his problem. for example,
c wt(l) = 1.0 specifies absolute error,
c = abs(y(l)) error relative to the most recent value of the
c l-th component of the solution,
c = abs(yp(l)) error relative to the most recent value of
c the l-th component of the derivative,
c = amax1(wt(l),abs(y(l))) error relative to the largest
c magnitude of l-th component obtained so far,
c = abs(y(l))*relerr/eps + abserr/eps specifies a mixed
c relative-absolute test where relerr is relative
c error, abserr is absolute error and eps =
c amax1(relerr,abserr) .
c
c subsequent calls --
c
c subroutine step is designed so that all information needed to
c continue the integration, including the step size h and the order
c k , is returned with each step. with the exception of the step
c size, the error tolerance, and the weights, none of the parameters
c should be altered. the array wt must be updated after each step
c to maintain relative error tests like those above. normally the
c integration is continued just beyond the desired endpoint and the
c solution interpolated there with subroutine intrp . if it is
c impossible to integrate beyond the endpoint, the step size may be
c reduced to hit the endpoint since the code will not take a step
c larger than the h input. changing the direction of integration,
c i.e., the sign of h , requires the user set start = .true. before
c calling step again. this is the only situation in which start
c should be altered.
c
c output from step
c
c successful step --
c
c the subroutine returns after each successful step with start and
c crash set .false. . x represents the independent variable
c advanced one step of length hold from its value on input and y
c the solution vector at the new value of x . all other parameters
c represent information corresponding to the new x needed to
c continue the integration.
c
c unsuccessful step --
c
c when the error tolerance is too small for the machine precision,
c the subroutine returns without taking a step and crash = .true. .
c an appropriate step size and error tolerance for continuing are
c estimated and all other information is restored as upon input
c before returning. to continue with the larger tolerance, the user
c just calls the code again. a restart is neither required nor
c desirable.
c
logical phase1,nornd
real twou, fouru, p5eps, round, sum, absh, realns,
* temp1, temp2, reali, temp3, temp4, temp5, temp6,
* tau, xold, erkm1, erkm2, erk, err, rho, erkp1,
* r, hnew
real alpha(12),beta(12),sig(13),w(12),v(12),g(13),
* gstr(13),two(13)
real abs, amax1, amin1, r1mach, sign, sqrt
c***********************************************************************
c* the only machine dependent constants are based on the machine unit *
c* roundoff error u which is the smallest positive number such that *
c* 1.0+u .gt. 1.0 . the user must calculate u and insert *
c* twou=2.0*u and fouru=4.0*u in the data statement before calling *
c* the code. the routine machin calculates u . *
c data twou,fouru/4.4d-16, 8.8d-16/ ***
c***********************************************************************
data two/2.0,4.0,8.0,16.0,32.0,64.0,128.0,256.0,
* 512.0,1024.0,2048.0,4096.0,8192.0/
data gstr/0.500,0.0833,0.0417,0.0264,0.0188,0.0143,
* 0.0144,0.00936,0.00789,0.00679,0.00592,
* 0.00524,0.00468/
data g(1),g(2)/1.0,0.5/,sig(1)/1.0/
c
c
c *** begin block 0 ***
c check if step size or error tolerance is too small for machine
c precision. if first step, initialize phi array and estimate a
c starting step size.
c ***
c
c if step size is too small, determine an acceptable one
c
twou = 2.0*r1mach(4) ***
fouru = 2.0*twou ***
crash = .true.
if(abs(h) .ge. fouru*abs(x)) go to 5
h = sign(fouru*abs(x),h)
return
5 p5eps = 0.5*eps
c
c if error tolerance is too small, increase it to an acceptable value
c
round = 0.0
do 10 l = 1,neqn
10 round = round + (y(l)/wt(l))**2
round = twou*sqrt(round)
if(p5eps .ge. round) go to 15
eps = 2.0*round*(1.0 + fouru)
return
15 crash = .false.
if(.not.start) go to 99
c
c initialize. compute appropriate step size for first step
c
call f(x,y,yp)
sum = 0.0
do 20 l = 1,neqn
phi(l,1) = yp(l)
phi(l,2) = 0.0
20 sum = sum + (yp(l)/wt(l))**2
sum = sqrt(sum)
absh = abs(h)
if(eps .lt. 16.0*sum*h*h) absh = 0.25*sqrt(eps/sum)
h = sign(amax1(absh,fouru*abs(x)),h)
hold = 0.0
k = 1
kold = 0
start = .false.
phase1 = .true.
nornd = .true.
if(p5eps .gt. 100.0*round) go to 99
nornd = .false.
do 25 l = 1,neqn
25 phi(l,15) = 0.0
99 ifail = 0
c *** end block 0 ***
c
c *** begin block 1 ***
c compute coefficients of formulas for this step. avoid computing
c those quantities not changed when step size is not changed.
c ***
c
100 kp1 = k+1
kp2 = k+2
km1 = k-1
km2 = k-2
c
c ns is the number of steps taken with size h, including the current
c one. when k.lt.ns, no coefficients change
c
if(h .ne. hold) ns = 0
if (ns .le. kold) ns=ns+1
nsp1 = ns+1
if (k .lt. ns) go to 199
c
c compute those components of alpha(*),beta(*),psi(*),sig(*) which
c are changed
c
beta(ns) = 1.0
realns = ns
alpha(ns) = 1.0/realns
temp1 = h*realns
sig(nsp1) = 1.0
if(k .lt. nsp1) go to 110
do 105 i = nsp1,k
im1 = i-1
temp2 = psi(im1)
psi(im1) = temp1
beta(i) = beta(im1)*psi(im1)/temp2
temp1 = temp2 + h
alpha(i) = h/temp1
reali = i
105 sig(i+1) = reali*alpha(i)*sig(i)
110 psi(k) = temp1
c
c compute coefficients g(*)
c
c initialize v(*) and set w(*). g(2) is set in data statement
c
if(ns .gt. 1) go to 120
do 115 iq = 1,k
temp3 = iq*(iq+1)
v(iq) = 1.0/temp3
115 w(iq) = v(iq)
go to 140
c
c if order was raised, update diagonal part of v(*)
c
120 if(k .le. kold) go to 130
temp4 = k*kp1
v(k) = 1.0/temp4
nsm2 = ns-2
if(nsm2 .lt. 1) go to 130
do 125 j = 1,nsm2
i = k-j
125 v(i) = v(i) - alpha(j+1)*v(i+1)
c
c update v(*) and set w(*)
c
130 limit1 = kp1 - ns
temp5 = alpha(ns)
do 135 iq = 1,limit1
v(iq) = v(iq) - temp5*v(iq+1)
135 w(iq) = v(iq)
g(nsp1) = w(1)
c
c compute the g(*) in the work vector w(*)
c
140 nsp2 = ns + 2
if(kp1 .lt. nsp2) go to 199
do 150 i = nsp2,kp1
limit2 = kp2 - i
temp6 = alpha(i-1)
do 145 iq = 1,limit2
145 w(iq) = w(iq) - temp6*w(iq+1)
150 g(i) = w(1)
199 continue
c *** end block 1 ***
c
c *** begin block 2 ***
c predict a solution p(*), evaluate derivatives using predicted
c solution, estimate local error at order k and errors at orders k,
c k-1, k-2 as if constant step size were used.
c ***
c
c change phi to phi star
c
if(k .lt. nsp1) go to 215
do 210 i = nsp1,k
temp1 = beta(i)
do 205 l = 1,neqn
205 phi(l,i) = temp1*phi(l,i)
210 continue
c
c predict solution and differences
c
215 do 220 l = 1,neqn
phi(l,kp2) = phi(l,kp1)
phi(l,kp1) = 0.0
220 p(l) = 0.0
do 230 j = 1,k
i = kp1 - j
ip1 = i+1
temp2 = g(i)
do 225 l = 1,neqn
p(l) = p(l) + temp2*phi(l,i)
225 phi(l,i) = phi(l,i) + phi(l,ip1)
230 continue
if(nornd) go to 240
do 235 l = 1,neqn
tau = h*p(l) - phi(l,15)
p(l) = y(l) + tau
235 phi(l,16) = (p(l) - y(l)) - tau
go to 250
240 do 245 l = 1,neqn
245 p(l) = y(l) + h*p(l)
250 xold = x
x = x + h
absh = abs(h)
call f(x,p,yp)
c
c estimate errors at orders k,k-1,k-2
c
erkm2 = 0.0
erkm1 = 0.0
erk = 0.0
do 265 l = 1,neqn
temp3 = 1.0/wt(l)
temp4 = yp(l) - phi(l,1)
if(km2)265,260,255
255 erkm2 = erkm2 + ((phi(l,km1)+temp4)*temp3)**2
260 erkm1 = erkm1 + ((phi(l,k)+temp4)*temp3)**2
265 erk = erk + (temp4*temp3)**2
if(km2)280,275,270
270 erkm2 = absh*sig(km1)*gstr(km2)*sqrt(erkm2)
275 erkm1 = absh*sig(k)*gstr(km1)*sqrt(erkm1)
280 temp5 = absh*sqrt(erk)
err = temp5*(g(k)-g(kp1))
erk = temp5*sig(kp1)*gstr(k)
knew = k
c
c test if order should be lowered
c
if(km2)299,290,285
285 if(amax1(erkm1,erkm2) .le. erk) knew = km1
go to 299
290 if(erkm1 .le. 0.5*erk) knew = km1
c
c test if step successful
c
299 if(err .le. eps) go to 400
c *** end block 2 ***
c
c *** begin block 3 ***
c the step is unsuccessful. restore x, phi(*,*), psi(*) .
c if third consecutive failure, set order to one. if step fails more
c than three times, consider an optimal step size. double error
c tolerance and return if estimated step size is too small for machine
c precision.
c ***
c
c restore x, phi(*,*) and psi(*)
c
phase1 = .false.
x = xold
do 310 i = 1,k
temp1 = 1.0/beta(i)
ip1 = i+1
do 305 l = 1,neqn
305 phi(l,i) = temp1*(phi(l,i) - phi(l,ip1))
310 continue
if(k .lt. 2) go to 320
do 315 i = 2,k
315 psi(i-1) = psi(i) - h
c
c on third failure, set order to one. thereafter, use optimal step
c size
c
320 ifail = ifail + 1
temp2 = 0.5
if(ifail - 3) 335,330,325
325 if(p5eps .lt. 0.25*erk) temp2 = sqrt(p5eps/erk)
330 knew = 1
335 h = temp2*h
k = knew
if(abs(h) .ge. fouru*abs(x)) go to 340
crash = .true.
h = sign(fouru*abs(x),h)
eps = eps + eps
return
340 go to 100
c *** end block 3 ***
c
c *** begin block 4 ***
c the step is successful. correct the predicted solution, evaluate
c the derivatives using the corrected solution and update the
c differences. determine best order and step size for next step.
c ***
400 kold = k
hold = h
c
c correct and evaluate
c
temp1 = h*g(kp1)
if(nornd) go to 410
do 405 l = 1,neqn
rho = temp1*(yp(l) - phi(l,1)) - phi(l,16)
y(l) = p(l) + rho
405 phi(l,15) = (y(l) - p(l)) - rho
go to 420
410 do 415 l = 1,neqn
415 y(l) = p(l) + temp1*(yp(l) - phi(l,1))
420 call f(x,y,yp)
c
c update differences for next step
c
do 425 l = 1,neqn
phi(l,kp1) = yp(l) - phi(l,1)
425 phi(l,kp2) = phi(l,kp1) - phi(l,kp2)
do 435 i = 1,k
do 430 l = 1,neqn
430 phi(l,i) = phi(l,i) + phi(l,kp1)
435 continue
c
c estimate error at order k+1 unless
c in first phase when always raise order,
c already decided to lower order,
c step size not constant so estimate unreliable
c
erkp1 = 0.0
if(knew .eq. km1 .or. k .eq. 12) phase1 = .false.
if(phase1) go to 450
if(knew .eq. km1) go to 455
if(kp1 .gt. ns) go to 460
do 440 l = 1,neqn
440 erkp1 = erkp1 + (phi(l,kp2)/wt(l))**2
erkp1 = absh*gstr(kp1)*sqrt(erkp1)
c
c using estimated error at order k+1, determine appropriate order
c for next step
c
if(k .gt. 1) go to 445
if(erkp1 .ge. 0.5*erk) go to 460
go to 450
445 if(erkm1 .le. amin1(erk,erkp1)) go to 455
if(erkp1 .ge. erk .or. k .eq. 12) go to 460
c
c here erkp1 .lt. erk .lt. amax1(erkm1,erkm2) else order would have
c been lowered in block 2. thus order is to be raised
c
c raise order
c
450 k = kp1
erk = erkp1
go to 460
c
c lower order
c
455 k = km1
erk = erkm1
c
c with new order determine appropriate step size for next step
c
460 hnew = h + h
if(phase1) go to 465
if(p5eps .ge. erk*two(k+1)) go to 465
hnew = h
if(p5eps .ge. erk) go to 465
temp2 = k+1
r = (p5eps/erk)**(1.0/temp2)
hnew = absh*amax1(0.5,amin1(0.9,r))
hnew = sign(amax1(hnew,fouru*abs(x)),h)
465 h = hnew
return
c *** end block 4 ***
end
subroutine intrp(x,y,xout,yout,ypout,neqn,kold,phi,psi)
c
c the methods in subroutine step approximate the solution near x
c by a polynomial. subroutine intrp approximates the solution at
c xout by evaluating the polynomial there. information defining this
c polynomial is passed from step so intrp cannot be used alone.
c
c this code is completely explained and documented in the text,
c computer solution of ordinary differential equations: the initial
c value problem by l. f. shampine and m. k. gordon.
c
c input to intrp --
c
c all floating point variables are real
c the user provides storage in the calling program for the arrays in
c the call list
dimension y(neqn),yout(neqn),ypout(neqn),phi(neqn,16),psi(12)
c and defines
c xout -- point at which solution is desired.
c the remaining parameters are defined in step and passed to intrp
c from that subroutine
c
c output from intrp --
c
c yout(*) -- solution at xout
c ypout(*) -- derivative of solution at xout
c the remaining parameters are returned unaltered from their input
c values. integration with step may be continued.
c
dimension g(13),w(13),rho(13)
data g(1)/1.0/,rho(1)/1.0/
c
hi = xout - x
ki = kold + 1
kip1 = ki + 1
c
c initialize w(*) for computing g(*)
c
do 5 i = 1,ki
temp1 = i
5 w(i) = 1.0/temp1
term = 0.0
c
c compute g(*)
c
do 15 j = 2,ki
jm1 = j - 1
psijm1 = psi(jm1)
gamma = (hi + term)/psijm1
eta = hi/psijm1
limit1 = kip1 - j
do 10 i = 1,limit1
10 w(i) = gamma*w(i) - eta*w(i+1)
g(j) = w(1)
rho(j) = gamma*rho(jm1)
15 term = psijm1
c
c interpolate
c
do 20 l = 1,neqn
ypout(l) = 0.0
20 yout(l) = 0.0
do 30 j = 1,ki
i = kip1 - j
temp2 = g(i)
temp3 = rho(i)
do 25 l = 1,neqn
yout(l) = yout(l) + temp2*phi(l,i)
25 ypout(l) = ypout(l) + temp3*phi(l,i)
30 continue
do 35 l = 1,neqn
35 yout(l) = y(l) + hi*yout(l)
return
end