Discussion

This code is really just a wrapper.

The variables are similar to the ones used in the previous Runge Kutta subroutine:

y

the equation that is being solved is of the form

$$\frac{\textrm{d} \boldsymbol{y}(x)}{\textrm{d} x} = \textrm{derivs}(x, \boldsymbol{y})$$

so this variable represents a vector of y_k values, where k runs through the equations. On entry y stands for $\boldsymbol{y}\left(x_n\right)$. On exit it will stand for $\boldsymbol{y}\left(x_n + \Delta x\right)$.

dydx

is the initial value of derivs(x, y), that is $\textrm{derivs} \left(x_n, \boldsymbol{y}\left( x_n \right)\right)$.

x

is x_n

htry

is the initial guess for a good $\Delta x$ that may be changed if a requested accuracy of integration is not met.

eps

The error associated with the whole system of equations is going to be evaluated in the following manner. The Runge-Kutta subroutine rkck to be called by rkqs will return a new value $\boldsymbol{y}\left(x_n + \Delta x\right)$ plus an absolute error estimate $\Delta\boldsymbol{y}$. That error is then going to be scaled by dividing by a user supplied array yscal. This new scaled error is then going to be compared to $\epsilon$, which is what this parameter eps stands for.

yscal

is the scaling array. If you're happy with unscaled values of y simply set it to 1.

hdid

is the value of $\Delta x$ that has really been used, after all the mucking up, to make the step.

hnext

is the suggested next value of $\Delta x$.

derivs

is the subroutine used to evaluate the right hand side of

$$\frac{\textrm{d}\boldsymbol{y}(x)} {\textrm{d} x} = \boldsymbol{f}(x, \boldsymbol{y})$$

As most other codes discussed in this lecture notes the action here begins by checking that the dimensions of variables passed to the subroutine are correct. The dimension of y that is extracted is discarded:

  ndum=assert_eq(SIZE(y), SIZE(dydx), SIZE(yscal), 'rkqs')

Then the step size $\Delta x$ is set to the suggested step size of htry:

  h=htry

and we enter the DO loop within which we

1.

call rkck to evaluate the next y for the suggested value of $\Delta x$, but also to give us the truncation error in yerr:

     CALL rkck(y, dydx, x, h, ytemp, yerr, derivs)

2.

scale the returned error and compare it to $\epsilon$:

     errmax=MAXVAL(ABS(yerr(:)/yscal(:)))/eps
     IF (errmax <= 1.0) EXIT

observe that since we are working with the system of equations here, and each equation is going to have its own different value of error $\Delta y_k$, we pick up the largest error and make that stand for the error for the whole system of equations. If the error is within the prescribed limits we accept the result and exit the DO loop.

3.

If the error is too large, then we have to repeat the step with a shorter $\Delta x$. So we shrink the step, but by no more than a factor of 10:

     htemp=safety*h*(errmax**pshrnk)
     h = SIGN(MAX(ABS(htemp), 0.1_sp*ABS(h)), h)

4.

then check if the step size hasn't shrunk too much, in which case the subroutine aborts with an error message:

     xnew = x+h
     IF (xnew == x) CALL nrerror('stepsize underflow in rkqs')

5.

and return to the top of the loop when the next trial Runge Kutta step is made.

After we have finished with the looping and have the new values of y as well as the new value of error, we suggest stretching $\Delta x$, but by no more than 5 times:

  IF (errmax > errcon) THEN
     hnext=safety*h*(errmax**pgrow)
  ELSE
     hnext=5.0_sp*h
  END IF

The last three lines move the used $\Delta x$ into hdid, the new value of y goes into y, and x itself becomes updated to $x + \Delta x$:

  hdid=h
  x=x+h
  y(:)=ytemp(:)

The subroutine rkck is going to be very similar to subroutine rk4, but instead of making just one Runge-Kutta step, it has to make

1.

one full size step

2.

two half size steps

then compare the results and return the result of the two half size steps in ytemp and the observed error in yerr.

I leave it to you to develop your own version of rkck. The easiest way to do that is to rewrite subroutine rk4, which we have discussed in section 4.4.1.