Abstract

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

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  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.

1989

1987

J. M. Fine, “Interpolants for Runge–Kutta–Nystrom methods,” Computing 39, 27–42 (1987).
[CrossRef]

1986

1982

Caldwell, J. B.

Fine, J. M.

J. M. Fine, “Interpolants for Runge–Kutta–Nystrom methods,” Computing 39, 27–42 (1987).
[CrossRef]

Forbes, G. W.

Ghatak, A. K.

Henrici, P.

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.

Hildebrand, F. B.

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

Kumar, D. V.

Moore, D. T.

Sharma, A.

Synge, J. L.

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.

Appl. Opt.

Computing

J. M. Fine, “Interpolants for Runge–Kutta–Nystrom methods,” Computing 39, 27–42 (1987).
[CrossRef]

J. Opt. Soc. Am. A

Other

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.

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.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (3)

Fig. 1
Fig. 1

Total transverse error at the image plane as a function of step size for Sharma’s method. The lens is a three-element photographic objective with an inhomogeneous front element. It has a focal length of 12.5 mm, is f/2.8, and has a half-field of view of 29.4°. The ray used in this example is initially parallel to the axis at a height of 2.3 mm. The dashed line shows the error that would result if the intersection error were zero. Notice that since this is a fourth-order scheme, the transfer error decreases by 4 orders of magnitude whenever the step size is reduced by 1 order of magnitude.

Fig. 2
Fig. 2

Schematic representation of the intermediate trajectories discussed in Sections 2 and 3.

Fig. 3
Fig. 3

These curves represent the total transverse error at the image plane as a function of step size for Sharma’s method and for the new, optimal method. The configuration of the ray and the optical system are identical to those described in Fig. 1. As remarked in the caption of Fig. 1, since the underlying method is of fourth order these curves have an average gradient equal to 4 on these log–log plots. Notice that the amplitude of the oscillations for Sharma’s method is constant on this plot, regardless of the step size, whereas the amplitude of the residual oscillations present in the optimal method decreases as the step size is decreased.

Equations (30)

Equations on this page are rendered with MathJax. Learn more.

d 2 r d t 2 = 1 2 n 2 ( r ) ,
R = ( x y z ) , T = ( d x / d t d y / d t d z / d t )
D = 1 2 ( n 2 / x n 2 / y n 2 / z )
d 2 R d t 2 = D [ R ( t ) ] .
R k + 1 = R k + Δ t [ T k + ( ) ( A + 2 B ) ] ,
T k + 1 = T k + ( ) ( A + 4 B + C ) ,
A = Δ t D ( R k ) , B = Δ t D [ R k + ½ Δ t T k + ( ) Δ t A ] , C = Δ t D [ R k + Δ t T k ( ½ ) Δ t B ] .
R k + 1 = R k + Δ t T k + 1 2 Δ t 2 D ( R k ) + 1 6 Δ t 3 d D ( R k ) d t + 1 24 Δ t 4 d 2 D ( R k ) d t 2 + O ( Δ t 5 ) = R ( Δ t ; k ) + O ( Δ t 5 ) ,
R ( t ; l ) = a 0 + a 1 t + a 2 t 2 + a 3 t 3 + 0 ( Δ t 4 ) ,
a 0 = R l , a 1 = T l , a 2 = [ 3 ( R l + 1 R l ) ( 2 T l + T l + 1 ) Δ t ] / ( Δ t ) 2 , a 3 = [ ( T l + 1 + T l ) Δ t 2 ( R l + 1 R l ) ] / ( Δ t ) 3 .
d T l d t = D ( R i ) ,
a 0 + a 1 t + a 2 t 2 + a 3 t 3 = R l + t T l + 1 2 t 2 D ( R l ) + 1 6 t 3 d D ( R l ) d t + 1 24 t 2 ( 2 Δ t t Δ t 2 ) d 2 D ( R l ) d t 2 + O ( Δ t 5 ) .
R ( t ; k ) = R k + t T k + 1 2 t 2 A Δ t + 1 6 t 3 ( 4 B C 3 A ) Δ t 2 + 1 24 t 4 [ 4 ( A + C 2 B ) ] Δ t 3 + O ( Δ t 5 ) .
A / Δ t = D ( R k ) , ( 4 B C 3 A ) Δ t 2 = d D ( R k ) d t + O ( Δ t 2 ) , [ 4 ( A C 2 B ) ] Δ t 3 = d 2 D ( R k ) d t 2 + O ( Δ t ) .
G = Δ t D ( R l + c 0 Δ t T l + ½ c 0 2 Δ t C ) ,
c 0 = 3 4 1 4 ( 9 + 16 τ τ 3 ) 1 / 2 .
T ( t int ; l ) = T l + b 1 A + b 2 B + b 3 C + b 4 G + O ( Δ t 5 ) ,
b 1 = 1 / 12 ( 2 c 0 3 ) ( τ 3 ) [ c 0 ( 4 τ 2 9 τ + 6 ) 3 τ ( τ 1 ) 2 ] , b 2 = 2 τ 3 ( c 0 1 ) [ 3 τ ( τ 2 + 13 ) 4 ( 5 τ 2 + 6 ) ] 3 [ τ ( 3 τ 4 ) + ( 3 2 τ ) ( 3 2 c 0 ) ] , b 3 = τ 2 [ c 0 ( 4 τ 3 ) + τ ( 2 3 τ ) ] 6 ( c 0 1 ) , b 4 = 1 / 12 τ ( τ 1 ) ( τ 3 ) 2 ( 2 c 0 3 ) .
F ( x , y , z ) = 0 ,
F 0 = F ( R l ) , F 1 = F ( R l + 1 ) , F ˙ 0 = d F [ R ( t ; l ) ] dt | t = 0 = F ( R l ) · T l .
F [ R ( t ; l ) ] = F 0 + t F ˙ 0 t 2 Q + O ( Δ t 3 ) ,
Q = F 0 F 1 + Δ t F ˙ 0 Δ t 2 .
t app = 2 F 0 F ˙ 0 sign ( F 0 ) ( F ˙ 0 2 + 4 F 0 Q ) 1 / 2 ,
sign ( x ) = { 1 x > 0 1 x < 0 .
t int = t app F [ R ( t app ; l ) ] { d F [ R ( t app ; l ) ] d t } 1 ,
d F [ R ( t app ; l ) ] d t = F [ R ( t app ; l ) ] · d R ( t app ; l ) d t
G=ΔtD[Rn+Δt(c0Tn+c1A+c2B+c3C)].
Tn+b1A+b2B+b3C+b4G
c0=34+14(9+16ττ3)1/2,
b2=2τ2[c0(64τ)+τ(3τ4)]3(2c01),b3=τ3(2c01)(6τ332τ2+37τ12)6[2c0(34τ)(6τ216τ+9)].

Metrics