12 subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs,ierror)
15 double precision eps,h1,hmin,x1,x2,ystart(nvar)
16 double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
17 double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
21 integer,
intent(out) :: ierror
25 COMMON /path/ kmax,kount,dxsav,xp,yp
38 if (kmax.gt.0) xsav=x-2.d0*dxsav
43 yscal(i)=abs(y(i))+abs(h*dydx(i))+
tiny
47 if(abs(x-xsav).gt.abs(dxsav))
then
48 if(kount.lt.kmax-1)
then
59 if((x+h-x2)*(x+h-x1).gt.0.d0) h=x2-x
61 if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. dydx(1:nvar)) .or. x .ne. x)
then
62 write(*,*)
"ODEINT BEFORE CALL RKQS: NaN in x, dydx, or y!"
65 write(*,*)
"dydx",dydx
66 write(*,*)
"exiting..."
71 call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
73 if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. dydx(1:nvar)) .or. x .ne. x)
then
74 write(*,*)
"ODEINT AFTER CALL RKQS: NaN in x, dydx, or y!"
77 write(*,*)
"dydx",dydx
78 write(*,*)
"exiting..."
89 if((x-x2)*(x2-x1).ge.0.d0)
then
105 if(abs(hnext).lt.hmin)
then
115 SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
118 double precision eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
123 double precision errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,PSHRNK,ERRCON
124 parameter(safety=0.9d0,pgrow=-.2d0,pshrnk=-.25d0,errcon=1.89d-4)
128 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x .ne. x)
then
129 write(*,*)
"RKQS BEFORE CALL RKCK: NaN in x, dydx, or y!"
132 write(*,*)
"dydx",dydx
133 write(*,*)
"exiting..."
137 1
call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
138 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x .ne. x)
then
139 write(*,*)
"RKQS AFTER CALL RKCK: NaN in x, dydx, or y!"
142 write(*,*)
"dydx",dydx
143 write(*,*)
"exiting..."
149 errmax=max(errmax,abs(yerr(i)/yscal(i)))
152 if(errmax.gt.1.d0)
then
153 h=safety*h*(errmax**pshrnk)
158 if(xnew.eq.x)
call stop(
'stepsize underflow in rkqs')
161 if(errmax.gt.errcon)
then
162 hnext=safety*h*(errmax**pgrow)
176 SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs)
179 double precision h,x,dydx(n),y(n),yerr(n),yout(n)
184 double precision ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX),ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51, &
186 B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6
187 parameter(a2=.2d0,a3=.3d0,a4=.6d0,a5=1.d0, &
188 a6=.875d0,b21=.2d0,b31=3.d0/40.d0,&
189 b32=9.d0/40.d0,b41=.3d0,b42=-.9d0,b43=1.2d0,&
190 b51=-11.d0/54.d0,b52=2.5d0,&
191 b53=-70.d0/27.d0,b54=35.d0/27.d0,b61=1631.d0/55296.d0, &
193 b63=575.d0/13824.d0,b64=44275.d0/110592.d0, &
194 b65=253.d0/4096.d0,c1=37.d0/378.d0,&
195 c3=250.d0/621.d0,c4=125.d0/594.d0,c6=512.d0/1771.d0,&
196 dc1=c1-2825.d0/27648.d0,&
197 dc3=c3-18575.d0/48384.d0,dc4=c4-13525.d0/55296.d0,&
198 dc5=-277.d0/14336.d0,&
200 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)))
then
201 write(*,*)
"NaNs IN RKCK, STEP 0!"
202 write(*,*)
"y0",y(1:n)
203 write(*,*)
"derivs",dydx(1:n)
204 write(*,*)
"ABORTING..."
208 ytemp(i)=y(i)+b21*h*dydx(i)
210 call derivs(x+a2*h,ytemp,ak2)
211 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak2(1:n) .ne. ak2(1:n)))
then
212 write(*,*)
"NaNs IN RKCK, STEP 1!"
213 write(*,*)
"y1",ytemp(1:n)
214 write(*,*)
"derivs",ak2(1:n)
215 write(*,*)
"ABORTING..."
219 ytemp(i)=y(i)+h*(b31*dydx(i)+b32*ak2(i))
221 call derivs(x+a3*h,ytemp,ak3)
222 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak3(1:n) .ne. ak3(1:n)))
then
223 write(*,*)
"NaNs IN RKCK, STEP 2!"
224 write(*,*)
"y2",ytemp(1:n)
225 write(*,*)
"derivs",ak3(1:n)
226 write(*,*)
"ABORTING..."
230 ytemp(i)=y(i)+h*(b41*dydx(i)+b42*ak2(i)+b43*ak3(i))
232 call derivs(x+a4*h,ytemp,ak4)
233 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak4(1:n) .ne. ak4(1:n)))
then
234 write(*,*)
"NaNs IN RKCK, STEP 3!"
235 write(*,*)
"y3",ytemp(1:n)
236 write(*,*)
"derivs",ak4(1:n)
237 write(*,*)
"ABORTING..."
241 ytemp(i)=y(i)+h*(b51*dydx(i)+b52*ak2(i)+b53*ak3(i)+b54*ak4(i))
243 call derivs(x+a5*h,ytemp,ak5)
244 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak5(1:n) .ne. ak5(1:n)))
then
245 write(*,*)
"NaNs IN RKCK, STEP 4!"
246 write(*,*)
"y4",ytemp(1:n)
247 write(*,*)
"derivs",ak5(1:n)
248 write(*,*)
"ABORTING..."
252 ytemp(i)=y(i)+h*(b61*dydx(i)+b62*ak2(i)+b63*ak3(i)+b64*ak4(i)+b65*ak5(i))
254 call derivs(x+a6*h,ytemp,ak6)
255 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak6(1:n) .ne. ak6(1:n)))
then
256 write(*,*)
"NaNs IN RKCK, STEP 1!"
257 write(*,*)
"y15",ytemp(1:n)
258 write(*,*)
"derivs",ak6(1:n)
259 write(*,*)
"ABORTING..."
263 yout(i)=y(i)+h*(c1*dydx(i)+c3*ak3(i)+c4*ak4(i)+c6*ak6(i))
266 yerr(i)=h*(dc1*dydx(i)+dc3*ak3(i)+dc4*ak4(i)+dc5*ak5(i)+dc6*ak6(i))