OSA's Digital Library

Journal of the Optical Society of America A

Journal of the Optical Society of America A


  • Vol. 7, Iss. 2 — Feb. 1, 1990
  • pp: 248–254

Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media

Bryan D. Stone and G. W. Forbes  »View Author Affiliations

JOSA A, Vol. 7, Issue 2, pp. 248-254 (1990)

View Full Text Article

Enhanced HTML    Acrobat PDF (901 KB)

Browse Journals / Lookup Meetings

Browse by Journal and Year


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools



Runge–Kutta integration schemes are well suited to the determination of ray trajectories in inhomogeneous media. There is a fundamental difference, however, between Runge–Kutta schemes and many other schemes for numerically integrating ordinary differential equations: Runge–Kutta schemes are not based on approximating the continuous trajectory by a polynomial. Consequently, these schemes do not implicitly provide a continuous trajectory; they yield only a series of points through which the ray passes, together with the ray direction at those points. A supplementary method must be devised when a continuous trajectory is required. The accuracy of a continuous trajectory for Runge–Kutta schemes is limited by the error introduced in a single iteration of the integrator. A trajectory that attains this limit is referred to here as optimal. The existing method of calculating trajectories for a widely used Runge–Kutta scheme is, in fact, not optimal. Accordingly, an efficient method of determining optimal intermediate trajectories is presented. This new technique is shown to be superior to the existing method for locating ray–surface intersections and allows accuracy doubling (a recently proposed method for accelerating the analysis of systems with inhomogeneous elements) to be used to full advantage.

© 1990 Optical Society of America

Original Manuscript: April 20, 1989
Revised Manuscript: September 7, 1989
Published: February 1, 1990

Bryan D. Stone and G. W. Forbes, "Optimal interpolants for Runge–Kutta ray tracing in inhomogeneous media," J. Opt. Soc. Am. A 7, 248-254 (1990)

Sort:  Author  |  Year  |  Journal  |  Reset  


  1. 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.
  2. A. Sharma, D. V. Kumar, A. K. Ghatak, “Tracing rays through graded-index media: a new method,” Appl. Opt. 21, 984–987 (1982). [CrossRef] [PubMed]
  3. 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.
  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.
  5. G. W. Forbes, “Accuracy doubling in the determination of final ray configurations,” J. Opt. Soc. Am. A 6, 1776–1783 (1989). [CrossRef]
  6. J. M. Fine, “Interpolants for Runge–Kutta–Nystrom methods,” Computing 39, 27–42 (1987). [CrossRef]
  7. A. Sharma, A. K. Ghatak, “Ray tracing in gradient-index lenses: computation of ray–surface intersection,” Appl. Opt. 25, 3409–3412 (1986). [CrossRef] [PubMed]
  8. J. B. Caldwell, D. T. Moore, “Design of gradient-index lens systems for disc format cameras,” Appl. Opt. 25, 3351–3355 (1986), lens GT2. [CrossRef] [PubMed]
  9. 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.
  10. 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.
  11. 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.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


Fig. 1 Fig. 2 Fig. 3

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited