subroutine ode(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
implicit real*8(a-h,o-z)
c
c double precision subroutine ode integrates a system of neqn
c first order ordinary differential equations of the form:
c dy(i)/dt = f(t,y(1),y(2),...,y(neqn))
c y(i) given at t .
c the subroutine integrates from t to tout . on return the
c parameters in the call list are set for continuing the integration.
c the user has only to define a new value tout and call ode again.
c
c the differential equations are actually solved by a suite of codes
c de , step , and intrp . ode allocates virtual storage in the
c arrays work and iwork and calls de . de is a supervisor which
c directs the solution. it calls on the routines step and intrp
c to advance the integration and to interpolate at output points.
c step uses a modified divided difference form of the adams pece
c formulas and local extrapolation. it adjusts the order and step
c size to control the local error per unit step in a generalized
c sense. normally each call to step advances the solution one step
c in the direction of tout . for reasons of efficiency de
c integrates beyond tout internally, though never beyond
c t+10*(tout-t), and calls intrp to interpolate the solution at
c tout . an option is provided to stop the integration at tout but
c it should be used only if it is impossible to continue the
c integration beyond tout .
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 the parameters represent:
c f -- double precision subroutine f(t,y,yp) to evaluate
c derivatives yp(i)=dy(i)/dt
c neqn -- number of equations to be integrated (integer*4)
c y(*) -- solution vector at t (real*8)
c t -- independent variable (real*8)
c tout -- point at which solution is desired (real*8)
c relerr,abserr -- relative and absolute error tolerances for local
c error test (real*8). at each step the code requires
c dabs(local error) .le. dabs(y)*relerr + abserr
c for each component of the local error and solution vectors
c iflag -- indicates status of integration (integer*4)
c work(*) (real*8) -- arrays to hold information internal to
c iwork(*) (integer*4) which is necessary for subsequent calls
c
c first call to ode --
c
c the user must provide storage in his calling program for the arrays
c in the call list,
c y(neqn), work(100+21*neqn), iwork(5),
c declare f in an external statement, supply the double precision
c subroutine f(t,y,yp) to evaluate
c dy(i)/dt = yp(i) = f(t,y(1),y(2),...,y(neqn))
c and 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 -- point at which solution is desired
c relerr,abserr -- relative and absolute local error tolerances
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 all parameters except f , neqn and tout may be altered by the
c code on output so must be variables in the calling program.
c
c output from ode --
c
c neqn -- unchanged
c y(*) -- solution at t
c t -- last point reached in integration. normal return has
c t = tout .
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 500 steps needed
c = 5 -- integration did not reach tout because equations
c appear to be stiff
c = 6 -- invalid input parameters (fatal error)
c the value of iflag is returned negative when the input
c value is negative and the integration does not reach tout ,
c i.e., -3, -4, -5.
c work(*),iwork(*) -- information generally of no interest to the
c user but necessary for subsequent calls.
c
c subsequent calls to ode --
c
c subroutine ode returns with all information needed to continue
c the integration. if the integration reached tout , the user need
c only define a new tout and call again. if the integration did not
c reach tout and the user wants to continue, he just calls again.
c the output value of iflag is the appropriate input value for
c subsequent calls. the only situation in which it should be altered
c is to stop the integration internally at the new tout , i.e.,
c change output iflag=2 to input iflag=-2 . error tolerances may
c be changed by the user before continuing. all other parameters must
c remain unchanged.
c
c***********************************************************************
c* subroutines de and step contain machine dependent constants. *
c* be sure they are set before using ode . *
c***********************************************************************
c
logical start,phase1,nornd
dimension y(neqn),work(1),iwork(5)
external f
data ialpha,ibeta,isig,iv,iw,ig,iphase,ipsi,ix,ih,ihold,istart,
1 itold,idelsn/1,13,25,38,50,62,75,76,88,89,90,91,92,93/
iyy = 100
iwt = iyy + neqn
ip = iwt + neqn
iyp = ip + neqn
iypout = iyp + neqn
iphi = iypout + neqn
if(iabs(iflag) .eq. 1) go to 1
start = work(istart) .gt. 0.0d0
phase1 = work(iphase) .gt. 0.0d0
nornd = iwork(2) .ne. -1
1 call de(f,neqn,y,t,tout,relerr,abserr,iflag,work(iyy),
1 work(iwt),work(ip),work(iyp),work(iypout),work(iphi),
2 work(ialpha),work(ibeta),work(isig),work(iv),work(iw),work(ig),
3 phase1,work(ipsi),work(ix),work(ih),work(ihold),start,
4 work(itold),work(idelsn),iwork(1),nornd,iwork(3),iwork(4),
5 iwork(5))
work(istart) = -1.0d0
if(start) work(istart) = 1.0d0
work(iphase) = -1.0d0
if(phase1) work(iphase) = 1.0d0
iwork(2) = -1
if(nornd) iwork(2) = 1
return
end
subroutine de(f,neqn,y,t,tout,relerr,abserr,iflag,
1 yy,wt,p,yp,ypout,phi,alpha,beta,sig,v,w,g,phase1,psi,x,h,hold,
2 start,told,delsgn,ns,nornd,k,kold,isnold)
implicit real*8(a-h,o-z)
c
c ode merely allocates storage for de to relieve the user of the
c inconvenience of a long call list. consequently de is used as
c described in the comments for ode .
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
logical stiff,crash,start,phase1,nornd
dimension y(neqn),yy(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),
1 ypout(neqn),psi(12),alpha(12),beta(12),sig(13),v(12),w(12),g(13)
external f
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 data statement before using de . the routine *
c* machin calculates u . fouru and twou=2.0*u must also be *
c* inserted in subroutine step before calling de . *
c data fouru/.888d-15/ ***
c***********************************************************************
c
c the constant maxnum is the maximum number of steps allowed in one
c call to de . the user may change this limit by altering the
c following statement
data maxnum/500/
c
c *** *** ***
c test for improper parameters
c
fouru = 4.0 * d1mach(4) ***
if(neqn .lt. 1) go to 10
if(t .eq. tout) go to 10
if(relerr .lt. 0.0d0 .or. abserr .lt. 0.0d0) go to 10
eps = dmax1(relerr,abserr)
if(eps .le. 0.0d0) 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 20
10 iflag = 6
return
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 = dabs(del)
tend = t + 10.0d0*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.0d0) go to 50
c
c on start and restart also set work variables x and yy(*), store the
c direction of integration and initialize the step size
c
30 start = .true.
x = t
do 40 l = 1,neqn
40 yy(l) = y(l)
delsgn = dsign(1.0d0,del)
h = dsign(dmax1(dabs(tout-x),fouru*dabs(x)),tout-x)
c
c if already past output point, interpolate and return
c
50 if(dabs(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. dabs(tout-x) .ge. fouru*dabs(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 many steps
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 = dsign(dmin1(dabs(h),dabs(tend-x)),h)
do 110 l = 1,neqn
110 wt(l) = releps*dabs(yy(l)) + abseps
call step(x,yy,f,neqn,h,eps,wt,start,
1 hold,k,kold,crash,phi,p,yp,psi,
2 alpha,beta,sig,v,w,g,phase1,ns,nornd)
c
c test for tolerances too small
c
if(.not.crash) go to 130
iflag = isn*3
relerr = eps*releps
abserr = eps*abseps
do 120 l = 1,neqn
120 y(l) = yy(l)
t = x
told = t
isnold = 1
return
c
c augment counter on number of steps and test for stiffness
c
130 nostep = nostep + 1
kle4 = kle4 + 1
if(kold .gt. 4) kle4 = 0
if(kle4 .ge. 50) stiff = .true.
go to 50
end
subroutine step(x,y,f,neqn,h,eps,wt,start,
c 1 hold,k,kold,crash,phi,p,yp,psi)
1 hold,k,kold,crash,phi,p,yp,psi,
2 alpha,beta,sig,v,w,g,phase1,ns,nornd)
implicit real*8(a-h,o-z)
c
c double precision subroutine step
c 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 (real*8)
c y(*) -- solution vector at x (real*8)
c yp(*) -- derivative of solution vector at x after successful
c step (real*8)
c neqn -- number of equations to be integrated (integer*4)
c h -- appropriate step size for next step. normally determined by
c code (real*8)
c eps -- local error tolerance. must be variable (real*8)
c wt(*) -- vector of weights for error criterion (real*8)
c start -- logical variable set .true. for first step, .false.
c otherwise (logical*4)
c hold -- step size used for last successful step (real*8)
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. all are real*8
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 = dabs(y(l)) error relative to the most recent value of
c the l-th component of the solution,
c = dabs(yp(l)) error relative to the most recent value of
c the l-th component of the derivative,
c = dmax1(wt(l),dabs(y(l))) error relative to the largest
c magnitude of l-th component obtained so far,
c = dabs(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 dmax1(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 start,crash,phase1,nornd
dimension y(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),psi(12)
dimension alpha(12),beta(12),sig(13),w(12),v(12),g(13),
1 gstr(13),two(13)
external f
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/.444d-15,.888d-15/ ***
c***********************************************************************
data two/2.0d0,4.0d0,8.0d0,16.0d0,32.0d0,64.0d0,128.0d0,256.0d0,
1 512.0d0,1024.0d0,2048.0d0,4096.0d0,8192.0d0/
data gstr/0.500d0,0.0833d0,0.0417d0,0.0264d0,0.0188d0,0.0143d0,
1 0.0114d0,0.00936d0,0.00789d0,0.00679d0,0.00592d0,0.00524d0,
2 0.00468d0/
c data g(1),g(2)/1.0,0.5/,sig(1)/1.0/
c
c
twou = 2.0 * d1mach(4) ***
fouru = 2.0 * twou ***
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
crash = .true.
if(dabs(h) .ge. fouru*dabs(x)) go to 5
h = dsign(fouru*dabs(x),h)
return
5 p5eps = 0.5d0*eps
c
c if error tolerance is too small, increase it to an acceptable value
c
round = 0.0d0
do 10 l = 1,neqn
10 round = round + (y(l)/wt(l))**2
round = twou*dsqrt(round)
if(p5eps .ge. round) go to 15
eps = 2.0*round*(1.0d0 + fouru)
return
15 crash = .false.
g(1)=1.0d0
g(2)=0.5d0
sig(1)=1.0d0
if(.not.start) go to 99
c
c initialize. compute appropriate step size for first step
c
call f(x,y,yp)
sum = 0.0d0
do 20 l = 1,neqn
phi(l,1) = yp(l)
phi(l,2) = 0.0d0
20 sum = sum + (yp(l)/wt(l))**2
sum = dsqrt(sum)
absh = dabs(h)
if(eps .lt. 16.0d0*sum*h*h) absh = 0.25d0*dsqrt(eps/sum)
h = dsign(dmax1(absh,fouru*dabs(x)),h)
hold = 0.0d0
k = 1
kold = 0
start = .false.
phase1 = .true.
nornd = .true.
if(p5eps .gt. 100.0d0*round) go to 99
nornd = .false.
do 25 l = 1,neqn
25 phi(l,15) = 0.0d0
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.0d0
realns = ns
alpha(ns) = 1.0d0/realns
temp1 = h*realns
sig(nsp1) = 1.0d0
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.0d0/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.0d0/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.0d0
220 p(l) = 0.0d0
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 = dabs(h)
call f(x,p,yp)
c
c estimate errors at orders k,k-1,k-2
c
erkm2 = 0.0d0
erkm1 = 0.0d0
erk = 0.0d0
do 265 l = 1,neqn
temp3 = 1.0d0/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)*dsqrt(erkm2)
275 erkm1 = absh*sig(k)*gstr(km1)*dsqrt(erkm1)
280 temp5 = absh*dsqrt(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(dmax1(erkm1,erkm2) .le. erk) knew = km1
go to 299
290 if(erkm1 .le. 0.5d0*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.0d0/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.5d0
if(ifail - 3) 335,330,325
325 if(p5eps .lt. 0.25d0*erk) temp2 = dsqrt(p5eps/erk)
330 knew = 1
335 h = temp2*h
k = knew
if(dabs(h) .ge. fouru*dabs(x)) go to 340
crash = .true.
h = dsign(fouru*dabs(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.0d0
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)*dsqrt(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.5d0*erk) go to 460
go to 450
445 if(erkm1 .le. dmin1(erk,erkp1)) go to 455
if(erkp1 .ge. erk .or. k .eq. 12) go to 460
c
c here erkp1 .lt. erk .lt. dmax1(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.0d0/temp2)
hnew = absh*dmax1(0.5d0,dmin1(0.9d0,r))
hnew = dsign(dmax1(hnew,fouru*dabs(x)),h)
465 h = hnew
return
c *** end block 4 ***
end
subroutine intrp(x,y,xout,yout,ypout,neqn,kold,phi,psi)
implicit real*8(a-h,o-z)
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 double precision
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.0d0/,rho(1)/1.0d0/
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.0d0/temp1
term = 0.0d0
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.0d0
20 yout(l) = 0.0d0
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