subroutine rk(ni,y,ysc,t,h,eps,hdone,ierr) C C This routine does Runge-Kutta integration of an ODE whose r.h.s. C is given by an external function func. C Parameters: C C ni [in]: number of equations. C y(ni) [in/out]: current value of the vector of solutions. C ysc(ni) [in]: error is scaled to ysc: abs(y) < eps*ysc. Has to C be positive. C t [in/out]: current values of the independent variable. C h [in/out]: on the input, the trial size of the time-step. On the C output, the new trial time-step for the next step. C eps [in]: the requested error relative to ysc. C hdone [out]: the time-step actually taken. C ierr [out]: error index: ierr=0 means success; ierr=-1 means C ni is too large, the internal parameter NE needs to C be increased; ierr>0 means a failure, in this case C ierr gives the number of the offending equation. C C This routine needs another routine with the following signature: C C subroutine func(ni,y,f,t) C C where: C C ni [in]: number of equations. C y(ni) [in]: current value of the vector of solutions. C f(ni) [out]: value of the r.h.s. C t [in]: current value of the independent variable. C C ******************************************************* C C RK used two sets of parameters. One of them has to be commented out. C C Feldberg set of parameters C parameter(a2 = 2.0/9.0) parameter(a3 = 1.0/3.0) parameter(a4 = 3.0/4.0) parameter(a5 = 1.0) parameter(a6 = 5.0/6.0) parameter(b21 = 2.0/9.0) parameter(b31 = 1.0/12.0) parameter(b32 = 1.0/4.0) parameter(b41 = 69.0/128.0) parameter(b42 = -243.0/128.0) parameter(b43 = 135.0/64.0) parameter(b51 = -17.0/12.0) parameter(b52 = 27.0/4.0) parameter(b53 = -54.0/10.0) parameter(b54 = 16.0/15.0) parameter(b61 = 65.0/432.0) parameter(b62 = -5.0/16.0) parameter(b63 = 13.0/16.0) parameter(b64 = 4.0/27.0) parameter(b65 = 5.0/144.0) parameter(c1 = 47.0/450.0) parameter(c2 = 0.0) parameter(c3 = 12.0/25.0) parameter(c4 = 32.0/225.0) parameter(c5 = 1.0/30.0) parameter(c6 = 6.0/25.0) parameter(d1 = 1.0/150.0) parameter(d2 = 0.0) parameter(d3 = -3.0/100.0) parameter(d4 = 16.0/75.0) parameter(d5 = 1.0/20.0) parameter(d6 = -6.0/25.0) C C Cash-Karp set of parameters C * parameter(a2 = 1.0/5.0) * parameter(a3 = 3.0/10.0) * parameter(a4 = 3.0/5.0) * parameter(a5 = 1.0) * parameter(a6 = 7.0/8.0) * parameter(b21 = 1.0/5.0) * parameter(b31 = 3.0/40.0) * parameter(b32 = 9.0/40.0) * parameter(b41 = 3.0/10.0) * parameter(b42 = -9.0/10.0) * parameter(b43 = 6.0/5.0) * parameter(b51 = -11.0/54.0) * parameter(b52 = 5.0/2.0) * parameter(b53 = -70.0/27.0) * parameter(b54 = 35.0/27.0) * parameter(b61 = 1631.0/55296.0) * parameter(b62 = 175.0/512.0) * parameter(b63 = 575.0/13824.0) * parameter(b64 = 44275/110592.0) * parameter(b65 = 253.0/4096.0) * parameter(c1 = 37.0/378.0) * parameter(c2 = 0.0) * parameter(c3 = 250.0/621.0) * parameter(c4 = 125.0/594.0) * parameter(c5 = 0.0) * parameter(c6 = 512/1771.0) * parameter(d1 = -277/64512.0) * parameter(d2 = 0.0) * parameter(d3 = 6925.0/370944.0) * parameter(d4 = -6925.0/202752.0) * parameter(d5 = -277.0/14336.0) * parameter(d6 = 277.0/7084.0) C C itmax is the allowed maximum number of attempts to decrease the time-step C and redo the calculation is the error tolerance is not satisfied. C parameter(itmax = 99) C C NE is the maximum number of equation accepted. C parameter(NE=1000) C C Input arrays C dimension y(*), ysc(*) C C Internal storage C dimension y0(NE), f0(NE), f(NE) dimension ak1(NE), ak2(NE), ak3(NE) dimension ak4(NE), ak5(NE), ak6(NE) C C First: set ierr to zero, check that the internal limits for the number C of equations is sufficient, and store the initial values of y C and t. Reset the iteration counter it. C ierr = 0 if(ni .gt. NE) then ierr = -1 return endif do j=1,ni y0(j) = y(j) enddo t0 = t it = 0 C C Calculate the r.h.s. at the initial values of y and t. Needs to be done C only once per call. C call func(ni,y0,f0,t0) C C Start the iteration loop. If the initial h is not sufficient, we come C back here, if it is, we do not come back. C 10 it = it + 1 C C Set k1 C do j=1,ni ak1(j) = h*f0(j) enddo C C Set k2 C t = t0 + a2*h do j=1,ni y(j) = y0(j) + b21*ak1(j) enddo call func(ni,y,f,t) do j=1,ni ak2(j) = h*f(j) enddo C C Set k3 C t = t0 + a3*h do j=1,ni y(j) = y0(j) + b31*ak1(j) + b32*ak2(j) enddo call func(ni,y,f,t) do j=1,ni ak3(j) = h*f(j) enddo C C Set k4 C t = t0 + a4*h do j=1,ni y(j) = y0(j) + b41*ak1(j) + b42*ak2(j) + b43*ak3(j) enddo call func(ni,y,f,t) do j=1,ni ak4(j) = h*f(j) enddo C C Set k5 C t = t0 + a5*h do j=1,ni y(j) = y0(j) + b51*ak1(j) + b52*ak2(j) + b53*ak3(j) + . b54*ak4(j) enddo call func(ni,y,f,t) do j=1,ni ak5(j) = h*f(j) enddo C C Set k6 C t = t0 + a6*h do j=1,ni y(j) = y0(j) + b61*ak1(j) + b62*ak2(j) + b63*ak3(j) + . b64*ak4(j) + b65*ak5(j) enddo call func(ni,y,f,t) do j=1,ni ak6(j) = h*f(j) enddo C C Find the new y and atthe same time find the maximum error errmax and C the equation number jmax where it is reached. C errmax = -1. jmax = 0 do j=1,ni y(j) = y0(j) + c1*ak1(j) + c2*ak2(j) + c3*ak3(j) + . c4*ak4(j) + c5*ak5(j) + c6*ak6(j) errabs = d1*ak1(j) + d2*ak2(j) + d3*ak3(j) + . d4*ak4(j) + d5*ak5(j) + d6*ak6(j) err = abs(errabs)/ysc(j) if(err .gt. errmax) then errmax = err jmax = j endif enddo C C Check the maximum error and the tolerance. C if(errmax .gt. eps) then C C Error is too large. Reduce the time-step, reset t, and start again C unless we have reached the iteration limit. In that case report the C failure. C t = t0 h = 0.9*h*(eps/errmax)**0.25 if(it .lt. itmax) goto 10 ierr = jmax hdone = 0.0 return else C C Error is acceptable. Update the independent variable, compute the C estimate for the next time-step, and exit. C t = t0 + h hdone = h h = 0.9*h*(eps/(0.004*eps+errmax))**0.2 endif C C We are done! C return end