Fresnel Integral Computation Techniques

The Fresnel integrals and their simultaneous parametric plot, the clothoid, have numerous applications including; but not limited to, optics and electromagnetic theory [8, 15, 16, 17, 21], robotics [2, 7, 12, 13, 14], civil engineering [4, 9, 20] and biology [18]. There have been numerous works over the past 70 years computing and numerically approximating Fresnel integrals. Boersma established approximations using the Lanczos tau-method [3] and Cody computed rational Chebyshev approximations using the Remes algorithm [5]. Another approach includes a spreadsheet computation by Mielenz [10]; which is based on successive improvements of known relational approximations. Mielenz also gives an improvement of his work [11], where the accuracy is less then $1.e$-$9$. More recently, Alazah, Chandler-Wilde and LaPorte propose a method to compute these integrals via a modified trapezoid rule [1].

Alazah et al. remark after some experimentation that a truncation of the Taylor series is more efficient and accurate than their new method for a small argument [1]. We build on this philosophy and add to it by endorsing asymptotic expansion for large arguments (see [10] for an earlier example). We propose a technique combining truncated Taylor series, the modified trapezium rule and asymptotic expansion and test its validity for practical computation of the Fresnel integrals all over the real line. Cut-off points are found analytically yielding a technique valid for arbitrary precision.

The paper is organized as follows: Section 2 gives a brief introduction to Fresnel integrals. In Section 3 the approximations with appropriate bounds are introduced. Section 4 contains a systematic method for choosing the parameters to construct an approximation using results from Section 3. Section 5 contains some benchmarking and comparative results. Lastly, conclusions are provided in Section 6.

We define the Fresnel cosine and sine integrals, $\mathrm{C}(x)$ and $\mathrm{S}(x)$ respectively, as follows:

This is consistent with [1]. With this pair of integrals, $\mathrm{C}(x)$ and $\mathrm{S}(x)$ can be written as a complex exponential integral:

The Argand diagram of $\mathrm{G}(x)$ for real $x$ plots the clothoid. Also, both $\mathrm{C}(x)$ and $\mathrm{S}(x)$ are odd functions of $x$ and can be extended to analytic functions over the complex plane. $\mathrm{G}(x)$ can also be expressed as

with $\text{erf}(x)=\frac{2}{\pi}\int_{0}^{x}\exp(-t^{2})\,\mathrm{d}t$ being the error function.

In this section we list the approximations to be used for $\mathrm{G}(x)$ along with their corresponding error bounds. We start with the truncated Taylor series $\mathrm{T}_{N}(x)$;

$T_{N}(x)$ is an alternating series in both real and imaginary parts bounding these by the next term, we have the estimate

Modified trapezoid sum from [1], which is

where $A_{N}$ is defined in (9). A global bound for eq. (6), see eq. (67) from [1], is

Where $c_{N}$ is given as eq. (47) in [1]; which is

and $\beta$, $A_{N}$ are the constants

The asymptotic expansion can be derived directly as follows, first write

Applying repeated applications of integration by parts gives

with

For the error bound, estimating the integral in (10) gives,

The expansion in eq. (11) is connected to the generalized hypergeometric series ${}_{2}F_{0}(1/2,1,-2\mathrm{i}/(\pi x^{2}))$. The bound in eq. (12) is a much sharper then directly estimating the remainder term of the hypergeometric series, see 13.7.5 from NIST [6].

We start this section by making some observations. The Taylor expansion $\mathrm{T}_{N}$ is accurate near zero. As $|x|$ becomes increasingly large, more terms are needed to maintain accuracy. On the other hand, the asymptotic expansion $\mathrm{Q}_{N}$ is accurate for $|x|$ large and as $|x|$ become small it needs more terms to maintain accuracy. The modified trapezoid rule does have a global bound; see [1], but in general is more prone to round-off errors for $|x|$ too small or too large. Given the parity of the Fresnel integrals, we consider positive $x$ and split the non-negative reals into 3 sub-intervals with cut-off points $x_{1}$ and $x_{2}$. Orders of approximation $N_{1},N_{2},N_{3}$ are also selected for the Taylor polynomial, $\mathrm{T}_{N}$, the modified trapezoid rule, $\mathrm{G}_{N}$, and the asymptotic expansion, $\mathrm{Q}_{N}$, respectively. We give a methodical approach to choosing the parameters $N_{1},N_{2},N_{3}$, $x_{1}$, $x_{2}$ and we select these parameters for double precision. In doing so we balance the float point operations between the 3 approximations and maintain the accuracy of the approximation. This provides and accurate approximation, near the desired machine precision with no increase in computational time for any $x>0$. We remark that if one would like another precision, the parameters can be refined in the same manner as done below.

We start with equation (8), which is a global bound for the modified trapezoid rule. Choosing the parameter $N_{2}=12$ gives $E_{\mathrm{G}_{12}}\approx 1.0733e$-$17$ which is below double precision. Choosing the value $N_{2}=11$ gives $E_{\mathrm{G}_{11}}\approx 2.301e$-$16$ which is slightly above double precision, for this example we will choose $N_{2}=12$. As per the above discussion, for a given $N$, $\mathrm{T}_{N}$ proves to be a better approximation for small $x$, while $\mathrm{Q}_{N}$ is a better approximation for large $x$. We thus choose $N_{1}$ and $N_{3}$ based on the computational time. Doing so gives $N_{1}=14$, $N_{3}=12$. In deciding the orders $N_{1},N_{2},N_{3}$, we can now choose $x_{1}$ and $x_{2}$ based off the precision we want. Considering eq. (5) at $x_{1}=0.688$ gives $E_{\mathrm{T}_{14}}(0.688)\approx 2.078e$-$16$. Also, considering eq. (12) at $x_{2}=6.725$ gives $E_{\mathrm{Q}_{12}}(6.725)\approx 2.212e$-$16$.

In summary, we have the approximation

Using the the parity of $\mathrm{G}(x)$ the estimate $|\mathrm{G}(x)-\tilde{\mathrm{G}}(x)|<2^{-52}$ holds for all $x\in\mathbb{R}$. So we achieve a double precision approximation that has similar computational time for all $x\in\mathbb{R}$ while avoiding machine rounding errors.

Alazah et al. (see [1]) test the validity of their method for arguments between 0 and 1000 and endorse the truncated Taylor series for small arguments, a well known approach. Many other methods present in literature use polynomial or rational approximations, splitting the domain into several regions. As far as we know, the most accurate one is Cody [5] with relative errors $\leq 10^{-15.58}$. However the form for the intermediate region and cut-off points are found experimentally. Mielenz [10] uses an asymptotic expansion for large arguments equivalent to a particular case of the form we use fixing $N=6$. We confirm the precision of our technique here and test its efficiency.

Numerical experiments were run with MatLab ver. R2020a using a 5th generation Intel i7-5500U CPU @ 2.4GHz, with 16GB of PC3L-12800 DDR3L ram. We compare our approximation against two sources. First the native Frensel cosine and sine functions for MatLab, we denote this by $\mathrm{G}_{m}(x)$ = fresnelc$(x)+\mathrm{i}$fresnels$(x)$ with the fresnelc and fresnels cosine and sine functions. We compare with an approximation based off Mielenz [10], which uses a spline interpolant. This will be denoted as $\mathrm{G}_{s}(x)$ and can be found in the MatLab file exchange¹¹1John D’Errico (2020). FresnelS and FresnelC
(https://www.mathworks.com/matlabcentral/fileexchange/28765-fresnels-and-fresnelc),
MATLAB Central File Exchange..

Figure 1 shows the errors between the approximation $\tilde{\mathrm{G}}(x)$, $\mathrm{G}_{m}(x)$ and $\mathrm{G}_{s}(x)$. For each subinterval, partitioned by $x_{1}=0.688$ and $x_{2}=6.725$, each approximation is evaluated at 1000 random points. The average of the errors is in Table 2. Although $\mathrm{G}_{m}(x)$ and $\mathrm{G}_{s}(x)$ are approximations themselves, there is consistent behavior between these approximations. The error reflects the accuracy of the approximation for a the given subinterval.

Figure 2 shows a comparison for large $x$, we use 2000 random test points between the interval $[10,1.e9]$. We show a comparison between the errors generated by the native functions in MatLab, $\mathrm{G}_{m}(x)$ and the asymptotic expansion $\mathrm{Q}_{12}(x)$, the errors generated by $\mathrm{G}_{m}(x)$ and $\mathrm{G}_{s}(x)$. This illustrates the limit of precision of the approximation $\mathrm{G}_{m}(x)$. We note that $\mathrm{G}_{s}(x)$ and $\tilde{\mathrm{G}}(x)$ use the same asymptotic expansion for large $x$, but their implementations are different. It should also be noted that both approximations $\mathrm{G}_{s}(x)$ and $\tilde{\mathrm{G}}(x)$ are fast to compute, the time can be seen in Table (1). If one were to need a quicker implementation, one could precompute the coefficients for given parameters $N_{1},N_{2},N_{3}$.

In this work we have presented a systematic way to compute Fresnel integrals up to an arbitrary machine precision. Our numerical experiments verify this for double precision. If more accuracy is required, one can select the parameters in a similar fashion to what is done in Section 4. In Section 5 we are also able to justify the results from previous works where justification was heuristic. With a proper implementation, the method presented can improve the performance for computing Fresnel integrals in standard packages such as SciPy [19]. In computer aided geometric design, this technique can be used to plot clothoid segments for transition curves and clothoidal splines.

The underlying philosophy combining Taylor series, some clever form of quadrature and asymptotic expansion can be investigated for the computation of other special functions and more exotic curves arising in CAGD.