subroutine gaussj(n, alpha, beta, b, t, w) c c this routine computes the nodes t(j) and weights c w(j) for gauss-jacobi quadrature formulas. c these are used when one wishes to approximate c c integral (from a to b) f(x) w(x) dx c c n c by sum w f(t ) c j=1 j j c c (w(x) and w(j) have no connection with each other) c where w(x) is the weight function c c w(x) = (1-x)**alpha * (1+x)**beta c c on (-1, 1), alpha, beta .gt. -1. c c input: c c n the number of points used for the quadrature rule c alpha see above c beta see above c b real scratch array of length n c c output parameters (both real arrays of length n) c c t will contain the desired nodes. c w will contain the desired weights w(j). c c subroutines required: solve, class, imtql2, gamma c c reference: c c the routine has been adapted from the more general c routine gaussq by golub and welsch. see c golub, g. h., and welsch, j. h., "calculation of gaussian c quadrature rules," mathematics of computation 23 (april, c 1969), pp. 221-230. c real b(n), t(n), w(n), muzero c call class (n, alpha, beta, b, t, muzero) w(1) = 1. do 105 i = 2, n 105 w(i) = 0.0 call imtql2 (n, t, b, w, ierr) do 110 i = 1, n 110 w(i) = muzero * w(i) * w(i) return end subroutine class(n, alpha, beta, b, a, muzero) c c this procedure supplies the coefficients a(j), b(j) of the c recurrence relation c c b p (x) = (x - a ) p (x) - b p (x) c j j j j-1 j-1 j-2 c c for the various classical (normalized) orthogonal polynomials, c and the zero-th moment c c muzero = integral w(x) dx c c of the given polynomial's weight function w(x). since the c polynomials are orthonormalized, the tridiagonal matrix is c guaranteed to be symmetric. c real a(n), b(n), muzero c nm1 = n - 1 c 50 ab = alpha + beta abi = 2. + ab muzero = 2. ** (ab + 1.) * gamma(alpha + 1.) * gamma( x beta + 1.) / gamma(abi) a(1) = (beta - alpha)/abi b(1) = sqrt(4.*(1. + alpha)*(1. + beta)/((abi + 1.)*abi*abi)) a2b2 = beta*beta - alpha*alpha do 51 i = 2, nm1 abi = 2.*i + ab a(i) = a2b2/((abi - 2.)*abi) 51 b(i) = sqrt (4.*i*(i + alpha)*(i + beta)*(i + ab)/ 1 ((abi*abi - 1)*abi*abi)) abi = 2.*n + ab a(n) = a2b2/((abi - 2.)*abi) return end subroutine imtql2(n, d, e, z, ierr) c c this is a modified version of the eispack routine imtql2. c it finds the eigenvalues and first components of the c eigenvectors of a symmetric tridiagonal matrix by the implicit ql c method. c real d(n), e(n), z(n), machep c c :::::::::: machep is a machine dependent parameter specifying c the relative precision of floating point arithmetic. data machep/1.e-7/ c ierr = 0 if (n .eq. 1) go to 1001 c e(n) = 0. do 240 l = 1, n j = 0 c :::::::::: look for small sub-diagonal element :::::::::: 105 do 110 m = l, n if (m .eq. n) go to 120 if (abs(e(m)) .le. machep * (abs(d(m)) + abs(d(m+1)))) x go to 120 110 continue c 120 p = d(l) if (m .eq. l) go to 240 if (j .eq. 30) go to 1000 j = j + 1 c :::::::::: form shift :::::::::: g = (d(l+1) - p) / (2. * e(l)) r = sqrt(g*g+1.) g = d(m) - p + e(l) / (g + sign(r, g)) s = 1. c = 1. p = 0. mml = m - l c c :::::::::: for i=m-1 step -1 until l do -- :::::::::: do 200 ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) if (abs(f) .lt. abs(g)) go to 150 c = g / f r = sqrt(c*c+1.) e(i+1) = f * r s = 1. / r c = c * s go to 160 150 s = f / g r = sqrt(s*s+1.) e(i+1) = g * r c = 1. / r s = s * c 160 g = d(i+1) - p r = (d(i) - g) * s + 2. * c * b p = s * r d(i+1) = g + p g = c * r - b c :::::::::: form first component of vector :::::::::: f = z(i+1) z(i+1) = s * z(i) + c * f 200 z(i) = c * z(i) - s * f c d(l) = d(l) - p e(l) = g e(m) = 0. go to 105 240 continue c c :::::::::: order eigenvalues and eigenvectors :::::::::: do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p p = z(i) z(i) = z(k) z(k) = p 300 continue c go to 1001 c :::::::::: set error -- no convergence to an c eigenvalue after 30 iterations :::::::::: 1000 ierr = l 1001 return end function gamma (x) c c computes the gamma function. needed by gaussj. c this function is adapted from an imsl routine but should c be replaced by any locally available alternative. c logical mflag dimension p(5),q(4),p4(3) c coefficients for minimax c approximation to gamma(x), c 2.0 .le. x .le. 3.0 data p(1)/-51.49952/,p(2)/80.05398/, * p(3)/-201.4659/,p(4)/-1.889439/, * p(5)/9.895546/ data q(1)/130.5263/,q(2)/-303.5898/, * q(3)/26.84174/,q(4)/-19.52375/ c approximation to ln(gamma(x)), c 12.0 .le. x data p4(1)/.9189385/,p4(2)/.8333332e-01/, * p4(3)/-.2770927e-02/ data iend/5/,iend1/4/,iend2/3/ data xinf/1.7014e+38/ data pi/3.141593/ c gamma(xmin) .approx. xinf c gamma(big1) .approx. xinf data xmin/5.8775e-39/ data big1/34.844/ c first executable statement ier = 0 mflag = .false. t = x if (abs(t).gt.xmin) go to 5 ier = 130 gamma = xinf if (t.le.0.0) gamma = -xinf go to 9000 5 if (abs(t).lt.big1) go to 10 ier = 129 gamma = xinf go to 9000 10 if (t.gt.0.0) go to 25 c argument is negative mflag = .true. t = -t r = aint(t) sign = 1.0 if (amod(r,2.0) .eq. 0.0) sign = -1. r = t-r if (r.ne.0.0) go to 20 ier = 130 gamma = xinf if (sign.eq.-1.0) gamma = -xinf go to 9000 c argument is not a negative integer 20 r = pi/sin(r*pi)*sign t = t+1.0 c evaluate approximation for gamma(t) c t .gt. xmin 25 if (t.gt.12.0) go to 60 i = t a = 1.0 if (i.gt.2) go to 40 i = i+1 go to (30,35,50),i c 0.0 .lt. t .lt. 1.0 30 a = a/(t*(t+1.0)) t = t+2.0 go to 50 c 1.0 .le. t .lt. 2.0 35 a = a/t t = t+1.0 go to 50 c 3.0 .le. t .le. 12.0 40 do 45 j=3,i t = t-1.0 a = a*t 45 continue c 2.0 .le. t .le. 3.0 50 top = p(iend1)*t+p(iend) den = t+q(iend1) do 55 j=1,iend2 top = top*t+p(j) den = den*t+q(j) 55 continue y = (top/den)*a if (mflag) y = r/y gamma = y go to 9005 c t .gt. 12.0 60 top = alog(t) top = t*(top-1.0)-.5*top t = 1.0/t b = t*t y = (p4(3)*b+p4(2))*t+p4(1)+top y = exp(y) if (mflag) y = r/y gamma = y go to 9005 9000 continue write (6,201) 201 format (' *** error in gaussj gamma function routine') 9005 return end subroutine ns01a (n,x,f,ajinv,dstep,dmax,acc,maxfun,iprint,w, 2 calfun) c c this subroutine solves a system of nonlinear algebraic equations of c the form f(x) = 0, where f and x are vectors of length n. the c function f should have continuous first derivatives, but there is c no need to calculate these, because the algorithm automatically c computes finite difference approximations. the method used in c seeking a solution is a compromise between newton's method and c steepest descent. the iterative nature of the algorithm requires c that an initial estimate of the solution be supplied, but this can c be fairly poor and the process still will usually converge. this c is because the algorithm tends to take steepest descent steps when c it is far from the solution, switching to the newton direction for c rapid convergence as it nears the solution. c c program author: m. j. d. powell c c documentation : technical report aere-r.5947 c harwell, england c november, 1968 c c the variables in the calling sequence are defined as follows: c c n number of variables (equations). must have n>1. c c x vector of length n containing the variables of the c equations. on entry this should contain the initial c estimate of the solution. on return it will contain the c best estimate of the solution. c c f vector of length n used as working space. on return it c contains the values of the functions at the solution x. c c ajinv two dimensional array (20 by 20) used as working space. c on return it contains an approximation to the inverse of c the jacobian matrix of the system at the solution x. c c dstep step length for use in approximating first derivatives of c the functions by differences between function values. c c dmax generous estimate of the euclidean distance of the c solution from the initial estimate. used to limit maximum c size of a single correction step. must have dmax > dstep c c acc accuracy requirement. the iterative process is assumed to c have converged and a normal return takes place when an x c is found for which the euclidean norm of f(x) is less than c acc. c c maxfun maximum number of function evaluations (calls to calfun) c to be allowed. (generally 10*n is a sufficient number of c evaluations, although for some problems more may be c needed.) c c iprint flag which controls printing: c .lt.0 - no printout except error messages c .ge.0 - sum-of-squares error at each iteration c .gt.0 - x and f at each iteration c c w vector of length n*(2*n+5) which is used as working space c c calfun the subroutine defining the system of nonlinear equations c it must have the form c c subroutine calfun(n,x,f) c c where n, x, and f have the same meaning as given above. c calfun should set the components of f to be the function c values at x. must be declared external in the calling c program. c c notes: c c there are four possible error exits, each of which is accompanied c by an appropriate diagnostic message. these occur when c c 1. the gradient of f at x becomes so small that it is c predicted that steps much larger than dmax are needed. c c 2. the number of calls to calfun exceeds maxfun. c c 3. f(x) fails to decrease for n+4 consecutive iterations when c a decrease is predicted, and x is within dstep of the most c successful vector of variables. c c 4. f(x) fails to decrease even though a new (therefore, c presumably reliable) jacobian matrix has just been c obtained by finite differences, and x is within distance c dstep from the point where the jacobian was calculated. c c these error conditions can be caused by a number of factors, such c as too small values for maxfun or dmax, or very poor values for c dstep or the initial guess for x. in addition, the problem itself c may be so badly conditioned that rounding errors in the c computations make finding a solution impossible, or the equations c may be inconsistent so that no solution exists. c c another important fact to note is that only one step size, dstep, c is used for all the variables. this implies that all the variable c are expected to be of roughly the same order of magnitude. scaling c may be required to accomplish this, but this enhances the numerical c stability of the problem anyway and is definitely advisable. c c note that if the order of the system is greater than 20, then ajin c (as well as several other arrays not appearing in the calling c sequence) must have their dimensions increased. c c matrix inversion is performed by the linpack routines c sgefa followed by sgedi. these require in addition the c linpack basic linear algebra subroutines saxpy, sscal, c isamax, and sswap. c external calfun real x(1),f(1),lawork(20),ajinv(20,20),w(1),det(2) integer ipvt(20) c set various parameters maxc=0 c 'maxc' counts the number of calls of calfun nt=n+4 ntest=nt c 'nt' and 'ntest' cause an error return if f(x) does not decrease dtest=float(n+n)-0.5 c 'dtest' is used to maintain linear independence nx=n*n nf=nx+n nw=nf+n mw=nw+n ndc=mw+n nd=ndc+n c these parameters separate the working space array w fmin=0. c usually 'fmin' is the least calculated value of f(x), c and the best x is in w(nx+1) to w(nx+n) dd=0. c usually dd is the square of the current step length dss=dstep*dstep dm=dmax*dmax dmm=4.*dm is=5 c 'is' controls a 'go to' statement following a call of calfun tinc=1. c 'tinc' is used in the criterion to increase the step length c start a new page for printing if (iprint) 1,1,85 85 print 86 86 format ('1') c call the subroutine calfun 1 maxc=maxc+1 call calfun (n,x,f) c test for convergence fsq=0. do 2 i=1,n fsq=fsq+f(i)*f(i) 2 continue if (fsq .gt. acc**2) go to 4 c provide printing of final solution if requested 3 if (iprint) 5,999,6 6 print 7,maxc 7 format (/5x,'the final ns01a solution required', 1i5,' calls of calfun, and is') print 8,(i,x(i),f(i),i=1,n) 8 format (4x,'i',8x,'x(i)',13x,'f(i)'/(i5,2e17.8)) 999 print 9,maxc,fsq 9 format (' the sum-of-squares error at step',i4,' is',e11.4) 5 return c test for error return because f(x) does not decrease 4 go to (10,11,11,10,11),is 10 if (fsq-fmin) 15,20,20 20 if (dd-dss) 12,12,11 12 ntest=ntest-1 if (ntest) 13,14,11 14 print 16,nt 16 format (//5x,'error return from ns01a because',i5, 1' calls of calfun failed to improve the residuals') 17 do 18 i=1,n x(i)=w(nx+i) f(i)=w(nf+i) 18 continue fsq=fmin go to 3 c error return because a new jacobian is unsuccessful 13 print 19 19 format (//5x,'error return from ns01a because f(x) ', 1'failed to decrease using a new jacobian') go to 17 15 ntest=nt c test whether there have been maxfun calls of calfun 11 if (maxfun-maxc) 21,21,22 21 print 23,maxc 23 format (//5x,'error return from ns01a because there have been', 1i5,' calls of calfun') if (fsq-fmin) 3,17,17 c provide printing if requested 22 if (iprint) 24,998,25 25 print 26,maxc 26 format (/5x,'at the',i5,' th call of calfun we have') print 8,(i,x(i),f(i),i=1,n) 998 print 9,maxc,fsq 24 go to (27,28,29,87,30),is c store the result of the initial call of calfun 30 fmin=fsq do 31 i=1,n w(nx+i)=x(i) w(nf+i)=f(i) 31 continue c calculate a new jacobian approximation 32 ic=0 is=3 33 ic=ic+1 x(ic)=x(ic)+dstep go to 1 29 k=ic do 34 i=1,n w(k)=(f(i)-w(nf+i))/dstep k=k+n 34 continue x(ic)=w(nx+ic) if (ic-n) 33,35,35 c calculate the inverse of the jacobian and set the direction matrix 35 k=0 do 36 i=1,n do 37 j=1,n k=k+1 ajinv(i,j)=w(k) w(nd+k)=0. 37 continue w(ndc+k+i)=1. w(ndc+i)=1.+float(n-i) 36 continue c c invert matrix with linpack routines: lda = 20 info = 0 call sgefa(ajinv,lda,n,ipvt,info) if (info.ne.0) goto 44 ijob = 1 call sgedi(ajinv,lda,n,ipvt,det,lawork,ijob) c c start iteration by predicting the descent and newton minima 38 ds=0. dn=0. sp=0. do 39 i=1,n x(i)=0. f(i)=0. k=i do 40 j=1,n x(i)=x(i)-w(k)*w(nf+j) f(i)=f(i)-ajinv(i,j)*w(nf+j) k=k+n 40 continue ds=ds+x(i)*x(i) dn=dn+f(i)*f(i) sp=sp+x(i)*f(i) 39 continue c test whether a nearby stationary point is predicted if (fmin*fmin-dmm*ds) 41,41,42 c if so then return or revise jacobian 42 go to (43,43,44),is 44 print 45 45 format (//5x,'error return from ns01a because a nearby ', 1'stationary point of f(x) is predicted') go to 17 43 ntest=0 do 46 i=1,n x(i)=w(nx+i) 46 continue go to 32 c test whether to apply the full newton correction 41 is=2 if (dn-dd) 47,47,48 47 dd=max(dn,dss) ds=0.25*dn tinc=1. if (dn-dss) 49,58,58 49 is=4 go to 80 c calculate the length of the steepest descent step 48 k=0 dmult=0. do 51 i=1,n dw=0. do 52 j=1,n k=k+1 dw=dw+w(k)*x(j) 52 continue dmult=dmult+dw*dw 51 continue dmult=ds/dmult ds=ds*dmult*dmult c test whether to use the steepest descent direction if (ds-dd) 53,54,54 c test whether the initial value of dd has been set 54 if (dd) 55,55,56 55 dd=max(dss,min(dm,ds)) ds=ds/(dmult*dmult) go to 41 c set the multiplier of the steepest descent direction 56 anmult=0. dmult=dmult*sqrt(dd/ds) go to 98 c interpolate between the steepest descent and the newton direction 53 sp=sp*dmult anmult=(dd-ds)/((sp-ds)+sqrt((sp-dd)**2+(dn-dd)*(dd-ds))) dmult=dmult*(1.-anmult) c calculate the change in x and its angle with the first direction 98 dn=0. sp=0. do 57 i=1,n f(i)=dmult*x(i)+anmult*f(i) dn=dn+f(i)*f(i) sp=sp+f(i)*w(nd+i) 57 continue ds=0.25*dn c test whether an extra step is needed for independence if (w(ndc+1)-dtest) 58,58,59 59 if (sp*sp-ds) 60,58,58 c take the extra step and update the direction matrix 50 is=2 60 do 61 i=1,n x(i)=w(nx+i)+dstep*w(nd+i) w(ndc+i)=w(ndc+i+1)+1. 61 continue w(nd)=1. do 62 i=1,n k=nd+i sp=w(k) do 63 j=2,n w(k)=w(k+n) k=k+n 63 continue w(k)=sp 62 continue go to 1 c express the new direction in terms of those of the direction c matrix, and update the counts in w(ndc+1) etc. 58 sp=0. k=nd do 64 i=1,n x(i)=dw dw=0. do 65 j=1,n k=k+1 dw=dw+f(j)*w(k) 65 continue go to (68,66),is 66 w(ndc+i)=w(ndc+i)+1. sp=sp+dw*dw if (sp-ds) 64,64,67 67 is=1 kk=i x(1)=dw go to 69 68 x(i)=dw 69 w(ndc+i)=w(ndc+i+1)+1. 64 continue w(nd)=1. c reorder the directions so that kk is first if (kk-1) 70,70,71 71 ks=ndc+kk*n do 72 i=1,n k=ks+i sp=w(k) do 73 j=2,kk w(k)=w(k-n) k=k-n 73 continue w(k)=sp 72 continue c generate the new orthogonal direction matrix 70 do 74 i=1,n w(nw+i)=0. 74 continue sp=x(1)*x(1) k=nd do 75 i=2,n ds=sqrt(sp*(sp+x(i)*x(i))) dw=sp/ds ds=x(i)/ds sp=sp+x(i)*x(i) do 76 j=1,n k=k+1 w(nw+j)=w(nw+j)+x(i-1)*w(k) w(k)=dw*w(k+n)-ds*w(nw+j) 76 continue 75 continue sp=1./sqrt(dn) do 77 i=1,n k=k+1 w(k)=sp*f(i) 77 continue c calculate the next vector x, and predict the right hand sides 80 fnp=0. k=0 do 78 i=1,n x(i)=w(nx+i)+f(i) w(nw+i)=w(nf+i) do 79 j=1,n k=k+1 w(nw+i)=w(nw+i)+w(k)*f(j) 79 continue fnp=fnp+w(nw+i)**2 78 continue c call calfun using the new vector of variables go to 1 c update the step size 27 dmult=0.9*fmin+0.1*fnp-fsq if (dmult) 82,81,81 82 dd=max(dss,0.25*dd) tinc=1. if (fsq-fmin) 83,28,28 c try the test to decide whether to increase the step length 81 sp=0. ss=0. do 84 i=1,n sp=sp+abs(f(i)*(f(i)-w(nw+i))) ss=ss+(f(i)-w(nw+i))**2 84 continue pj=1.+dmult/(sp+sqrt(sp*sp+dmult*ss)) sp=min(4.,tinc,pj) tinc=pj/sp dd=min(dm,sp*dd) go to 83 c if f(x) improves store the new value of x 87 if (fsq-fmin) 83,50,50 83 fmin=fsq do 88 i=1,n sp=x(i) x(i)=w(nx+i) w(nx+i)=sp sp=f(i) f(i)=w(nf+i) w(nf+i)=sp w(nw+i)=-w(nw+i) 88 continue if (is-1) 28,28,50 c calculate the changes in f and in x 28 do 89 i=1,n x(i)=x(i)-w(nx+i) f(i)=f(i)-w(nf+i) 89 continue c update the approximations to j and to ajinv k=0 do 90 i=1,n w(mw+i)=x(i) w(nw+i)=f(i) do 91 j=1,n w(mw+i)=w(mw+i)-ajinv(i,j)*f(j) k=k+1 w(nw+i)=w(nw+i)-w(k)*x(j) 91 continue 90 continue sp=0. ss=0. do 92 i=1,n ds=0. do 93 j=1,n ds=ds+ajinv(j,i)*x(j) 93 continue sp=sp+ds*f(i) ss=ss+x(i)*x(i) f(i)=ds 92 continue dmult=1. if (abs(sp)-0.1*ss) 94,95,95 94 dmult=0.8 95 pj=dmult/ss pa=dmult/(dmult*sp+(1.-dmult)*ss) k=0 do 96 i=1,n sp=pj*w(nw+i) ss=pa*w(mw+i) do 97 j=1,n k=k+1 w(k)=w(k)+sp*x(j) ajinv(i,j)=ajinv(i,j)+ss*f(j) 97 continue 96 continue go to 38 end subroutine sgefa(a,lda,n,ipvt,info) integer ipvt(1) real a(lda,1) c c linpack routine. c sgefa factors a single precision matrix by gaussian elimination. c requires blas saxpy,sscal,isamax c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1./a(k,k) call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.) info = n return end subroutine sgedi(a,lda,n,ipvt,det,work,job) integer ipvt(1) real a(lda,1),det(2),work(1) c c linpack routine. c sgedi computes the determinant and inverse of a matrix c using the factors computed by sgeco or sgefa. c c requires blas saxpy,sscal,sswap c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1. det(2) = 0. ten = 10. do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (det(1) .eq. 0.) go to 60 10 if (abs(det(1)) .ge. 1.) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1. go to 10 20 continue 30 if (abs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1. go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n a(k,k) = 1./a(k,k) t = -a(k,k) call sscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0. call saxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = 0. 110 continue do 120 j = kp1, n t = work(j) call saxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call sswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end subroutine saxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c jack dongarra, linpack, 3/11/78. c real dx(1),dy(1) if(n.le.0)return if (da .eq. 0.) return if(incx.eq.1.and.incy.eq.1)go to 20 ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end subroutine sscal(n,da,dx,incx) c c scales a vector by a constant. c jack dongarra, linpack, 3/11/78. c real dx(1) if(n.le.0)return if(incx.eq.1)go to 20 nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end integer function isamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c real dx(1) isamax = 0 if( n .lt. 1 ) return isamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 ix = 1 dmax = abs(dx(1)) ix = ix + incx do 10 i = 2,n if(abs(dx(ix)).le.dmax) go to 5 isamax = i dmax = abs(dx(ix)) 5 ix = ix + incx 10 continue return 20 dmax = abs(dx(1)) do 30 i = 2,n if(abs(dx(i)).le.dmax) go to 30 isamax = i dmax = abs(dx(i)) 30 continue return end subroutine sswap (n,dx,incx,dy,incy) c c interchanges two vectors. c jack dongarra, linpack, 3/11/78. c real dx(1),dy(1) if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end subroutine ode(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork) c c 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 -- subroutine f(t,y,yp) to evaluate c derivatives yp(i)=dy(i)/dt c neqn -- number of equations to be integrated (integer) c y(*) -- solution vector at t (real) c t -- independent variable (real) c tout -- point at which solution is desired (real) c relerr,abserr -- relative and absolute error tolerances for local c error test (real). 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 (integer) c work(*) (real) -- arrays to hold information internal to c iwork(*) (integer) 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 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= 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(abs(iflag) .eq. 1) go to 1 start = work(istart) .gt. 0. phase1 = work(iphase) .gt. 0. 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. if(start) work(istart) = 1. work(iphase) = -1. if(phase1) work(iphase) = 1. 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) 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 . * data fouru/.3e-6/ 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 if(neqn .lt. 1) go to 10 if(t .eq. tout) go to 10 if(relerr .lt. 0. .or. abserr .lt. 0.) go to 10 eps = max(relerr,abserr) if(eps .le. 0.) go to 10 if(iflag .eq. 0) go to 10 isn = sign(1,iflag) iflag = abs(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 = abs(del) tend = t + 10.*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.) 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 = sign(1.,del) h = sign(max(abs(tout-x),fouru*abs(x)),tout-x) c c if already past output point, interpolate and return c 50 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 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 = sign(min(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, 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) c c 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) c y(*) -- solution vector at x (real) c yp(*) -- derivative of solution vector at x after successful c step (real) c neqn -- number of equations to be integrated (integer) c h -- appropriate step size for next step. normally determined by c code (real) c eps -- local error tolerance. must be variable (real) c wt(*) -- vector of weights for error criterion (real) c start -- logical variable set .true. for first step, .false. c otherwise (logical) c hold -- step size used for last successful step (real) 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 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 c the 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 = max(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 max(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 . * data twou,fouru/.15e-6,.3e-6/ c*********************************************************************** data two/2.,4.,8.,16.,32.,64.,128.,256., 1 512.,1024.,2048.,4096.,8192./ data gstr/0.50,0.0833,0.0417,0.0264,0.0188,0.0143, 1 0.0114,0.00936,0.00789,0.00679,0.00592,0.00524, 2 0.00468/ c 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 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. 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. + fouru) return 15 crash = .false. g(1)=1. g(2)=0.5 sig(1)=1. if(.not.start) go to 99 c c initialize. compute appropriate step size for first step c call f(x,y,yp) sum = 0. do 20 l = 1,neqn phi(l,1) = yp(l) phi(l,2) = 0. 20 sum = sum + (yp(l)/wt(l))**2 sum = sqrt(sum) absh = abs(h) if(eps .lt. 16.*sum*h*h) absh = 0.25*sqrt(eps/sum) h = sign(max(absh,fouru*abs(x)),h) hold = 0. k = 1 kold = 0 start = .false. phase1 = .true. nornd = .true. if(p5eps .gt. 100.*round) go to 99 nornd = .false. do 25 l = 1,neqn 25 phi(l,15) = 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. realns = ns alpha(ns) = 1./realns temp1 = h*realns sig(nsp1) = 1. 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./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./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. 220 p(l) = 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. erkm1 = 0. erk = 0. do 265 l = 1,neqn temp3 = 1./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(max(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./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. 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. min(erk,erkp1)) go to 455 if(erkp1 .ge. erk .or. k .eq. 12) go to 460 c c here erkp1 .lt. erk .lt. max(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./temp2) hnew = absh*max(0.5,min(0.9,r)) hnew = sign(max(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 single 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./,rho(1)/1./ 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./temp1 term = 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. 20 yout(l) = 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