See, for example, J. L. Synge, Geometrical Optics—An Introduction to Hamilton’s method, Vol. 37 of Cambridge Tracts in Mathematics and Mathematical Physics, F. Smithies, J. A. Todd, eds. (Cambridge U. Press, Cambridge, 1937;reprinted in 1962), Chap. 5, Sec. 20.

Note that this method is also referred to as a Runge–Kutta–Nystrom method. For a description of this method, see, for example, P. Henrici, Discrete Variable Methods in Ordinary Differential Equations (Wiley, New York, 1962), Sec. 4.2–4.

See, for example, F. B. Hildebrand, Introduction to Numerical Analysis, 2nd ed. (Dover, New York, 1987), Sees. 6.1–6.6 and 6.12.

Although the trajectory that results from the use of Sharma’s method is shown bending away from the true trajectory in Fig. 2, it is also possible for Sharma’s cubic to bend toward the true trajectory. When this happens, there is an apparent improvement that results from the use of Sharma’s method. This apparent improvement, however, is not reliable: there is no way to predict whether a given ray is helped or hindered by the added error. It is, in effect, like adding a random number to the transfer error.

The method described here is not unique. More generally, the equation for G [Eq. (15)] can be written asG=ΔtD[Rn+Δt(c0Tn+c1A+c2B+c3C)].Analogs of Eqs. (15)–(18) are then found by performing a Taylor-series expansion onTn+b1A+b2B+b3C+b4Gand matching the result to the Taylor expansion for T(tint; n) up to order tint5 by adjusting the parameters c0, c1c2,c3, b1, b2, b3, and b4. In this form, the problem is underspecified: there are more free parameters than equations of constraint. In the method presented here, this freedom was used to set c1and c2 to zero. The rest of the parameters are then determined uniquely, with the exception of c0, whose form is given by the solution of a quadratic. Notice that the expression in Eq. (16) for c0may be replaced by the second solution of the quadratic:c0=34+14(9+16ττ−3)1/2,where, again, τ= tint/Δt. In this case, c0is a monotonically decreasing function of τ, with c0(0) = 3/2 and c0(1) = 1.If this replacement is made, Eqs. (18) may still be used; however, both the numerators and denominators of the expressions for b2and b3go to zero (in such a way that the quotient remains finite) as τ approaches 1. This is a problem whenever tintapproaches Δt(i.e., τ approaches 1) when this alternative procedure is implemented on a computer. However, equivalent expressions for b2and b3can be written for which the denominators remain nonzero for 0 < τ< 1:b2=2τ2[c0(6−4τ)+τ(3τ−4)]3(2c0−1),b3=τ3(2c0−1)(6τ3−32τ2+37τ−12)6[2c0(3−4τ)−(6τ2−16τ+9)].In general, the results obtained by using one of these forms for c0are not significantly different from those obtained by using the other. The particular form used for this paper was chosen because the gradient is sampled in the region between the two Runge–Kutta points to be connected. With the alternative form of c0, the gradient is sampled beyond Rl+1.

The maximum intersection error of Sharma’s method can be seen in Fig. 3 to be approximately an order of magnitude larger than the error obtained by using the optimal trajectory. Since this is a fourth-order scheme, reduction of the step parameter by a factor of α has the effect of reducing the transfer error by a factor of α4. Consequently, reducing the step size by a factor of 101/4is equivalent, in this example, to eliminating the intersection error.