3
$\begingroup$

It turns out that SciPy's Fresnel values are wrong for complex arguments and large enough absolute value. I'm trying to fix that.

The implementation is based on Zhang/Jin, Computation of special functions, which in turn is based on Abramowitz/Stegun, Handbook of Mathematical Functions. There we find for the Fresnel S integral (7.3.10) $$ S(z) = \frac{1}{2} - f(z) \cos\left(\frac{\pi}{2} z^2\right) - g(z) \sin\left(\frac{\pi}{2} z^2\right) $$ for all $z$ with the auxiliary functions (7.3.5), (7.3.6) $$ \begin{split} f(z) &= \left[\frac{1}{2} - S(z) \right] \cos\left(\frac{\pi}{2} z^2\right) - \left[\frac{1}{2} - C(z) \right] \sin\left(\frac{\pi}{2} z^2\right),\\ g(z) &= \left[\frac{1}{2} - C(z) \right] \cos\left(\frac{\pi}{2} z^2\right) + \left[\frac{1}{2} - S(z) \right] \sin\left(\frac{\pi}{2} z^2\right). \end{split} $$

Computation of $S$ for large values is done via the asymptotic expansion of $f$ (7.3.27) $$ \DeclareMathOperator\arg{arg} \pi z f(z)\sim 1 + \sum_{m=1}^\infty (-1)^m \frac{1\cdot 3\cdot\cdots \cdot (4m-1)}{(\pi z^2)^{2m}}, \quad z\to\infty, |\arg(z)|<\frac{\pi}{2}. $$ This is where I have problems understanding the approximation. Consider $f(iz)$; for the asymptotic $f_a$, we have $f_a(iz) = -if_a(z)$. However no such thing is true for $f$ itself. Numerical computation via the representations $$ \DeclareMathOperator\erf{erf} \begin{split} S(z) &= \frac{1 + i}{4} \left[\erf\left(\frac{1 + i}{2} \sqrt{\pi} z\right) - i \erf\left(\frac{1 - i}{2} \sqrt{\pi} z\right)\right]\\ C(z) &= \frac{1 - i}{4} \left[\erf\left(\frac{1 + i}{2} \sqrt{\pi} z\right) + i \erf\left(\frac{1 - i}{2} \sqrt{\pi} z\right)\right] \end{split} $$ shows that the approximation incorrect for everything off of the real axis, but the issue here could be numerical instability in the $\erf$ represenation too.

My current guess is that the above infinite sum is valid only in the area $|\arg{z}|<\pi/4$, which would already change how we compute the Fresnel integral values significantly.

To finish things off, here's a cplot of $f$ (for smaller $|z|$):

enter image description here

$\endgroup$
5
  • 1
    $\begingroup$ The signum function is not defined for complex x. I'd try the expansions in the Digital Library of Mathematical Functions 7.12.2-7. $\endgroup$
    – skbmoore
    Oct 27, 2021 at 14:54
  • $\begingroup$ I've checked those, and it seems to me that the asymptotic expansion of the auxiliary function is indeed wrong for values away from the real axis. I've created a gist here. $\endgroup$ Oct 27, 2021 at 16:25
  • $\begingroup$ I only glanced at the code, so what I say may not be true in actuality. It seems you are testing the algorithm with somewhat small values of |z|, like around 5 or so. The formulas you tried, and the complex analog I pointed you to, are asymptotic. For a moderate |z|, you will need a sophisticated way to truncate the infinite series and may need dozens of terms. You might be able to use acceleration methods to converge on an answer. For |z| small enough, the asymptotic expansion is meaningless. In that case it is better to use a error function representation that takes complex arguments. $\endgroup$
    – skbmoore
    Oct 27, 2021 at 17:04
  • $\begingroup$ Matlab, which is close to SciPy (I think) has a user-submitted Fresnel integral with complex argument z. I have not downloaded or tested it. mathworks.com/matlabcentral/fileexchange/18021-fresnel-integral $\endgroup$
    – skbmoore
    Oct 27, 2021 at 17:11
  • $\begingroup$ Thanks @skbmoore for the comments. I tested values around 5 since 4.5 is the point where algorithms (such as the one in scipy) switch to the asymptotic mode. Results don't get any better for larger $|z|$. I really do think the expansion is not valid for the argument range given on NIST. $\endgroup$ Oct 27, 2021 at 17:55

1 Answer 1

2
$\begingroup$

I will consider the function $\operatorname{f}(z)$, the treatment of $\operatorname{g}(z)$ is analogous. By http://dlmf.nist.gov/7.12.ii, we have $$ \operatorname{f}(z) = \frac{1}{{\pi z}}\sum\limits_{m = 0}^{N - 1} {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} + R_N^{(\operatorname{f})} (z) $$ where \begin{align*} R_N^{(\operatorname{f})} (z) & = \frac{{( - 1)^N }}{{\pi \sqrt 2 }}\int_0^{ + \infty } {\frac{{e^{ - \pi z^2 t/2} t^{2N - 1/2} }}{{1 + t^2 }}dt} \\ & = \frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}\frac{1}{{\Gamma \left( {2N + \frac{1}{2}} \right)}}\int_0^{ + \infty } {\frac{{e^{ - s} s^{2N - 1/2} }}{{1 + s^2 /(\pi z^2 /2)^2 }}ds} \\ & = \frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}\Pi _{2N + 1/2} (\pi z^2 /2), \end{align*} provided that $|\arg z|<\frac{\pi}{4}$ and $N\geq 0$. Here $\Pi_p(w)$ denotes one of Dingle's basic terminants: $$ \Pi _p (w) = \frac{1}{{\Gamma (p)}}\int_0^{ + \infty } {\frac{{e^{ - s} s^{p - 1} }}{{1 + (s/w)^2 }}ds} $$ for $|\arg w|<\frac{\pi}{2}$ and by analytic continuation elswhere. Using the expression for $R_N^{(\operatorname{f})} (z)$ in terms of this terminant, we can extend $R_N^{(\operatorname{f})} (z)$ to the universal covering of $\mathbb C \setminus \left\{ 0\right\}$. Now employing the estimates for the basic terminant established in https://doi.org/10.1007/s10440-017-0099-0, we obtain the bound \begin{align*} \left| {R_N^{(\operatorname{f})} (z)} \right| \le &\; \left| {\frac{1}{{\pi z}}( - 1)^N \left( {\frac{1}{2}} \right)_{2N} \frac{1}{{(\pi z^2 /2)^{2N} }}} \right| \\ & \times \begin{cases} 1 & \text{ if } \; \left|\arg z\right| \leq \frac{\pi}{8}, \\ \min\!\Big(\left|\csc ( 4\arg z)\right|,1 + \cfrac{1}{2}\chi(2N+1/2)\Big) & \text{ if } \; \frac{\pi}{8} < \left|\arg z\right| \leq \frac{\pi}{4}, \\ \cfrac{\sqrt {2\pi (2N + 1/2)} }{2\left| {\sin (2\arg z)} \right|^{2N+1/2} } + 1 + \cfrac{1}{2}\chi (2N +1/2) & \text{ if } \; \frac{\pi}{4} < \left|\arg z\right| < \frac{\pi}{2}. \end{cases} \end{align*} Here $\chi(p) =\sqrt{\pi}\Gamma(p/2+1)/\Gamma(p/2+1/2)$ for $p>0$. It is seen that the asymptotic expansion of $\operatorname{f}(z)$ is valid in every closed sub-sector of $|\arg z|<\frac{\pi}{2}$ in the sense of Poincaré. The Stokes lines are $\arg z =\pm \frac{\pi}{4}$ where terms are swiched on that remain exponentially small compared to any negative power of $z$ as long as we stay away from the rays $\arg z =\pm \frac{\pi}{2}$ (the anti-Stokes lines).

To obtain a better result and reveal the exponentially small terms, we can use the functional equation $$ \Pi _p (w) = \pm \pi i\frac{{e^{ \mp \frac{\pi }{2}ip} }}{{\Gamma (p)}}w^p e^{ \pm iw} + \Pi _p (we^{ \mp \pi i} ) $$ where $p>0$ and $w$ is any element of the universal covering of $\mathbb C \setminus \left\{ 0\right\}$ (the Riemann surface of the logarithm). With this functional equation, we find, after some algebra, $$ R_N^{(\operatorname{f})} (ze^{ \mp \frac{\pi }{2}i} ) = \pm iR_N^{(\operatorname{f})} (z) + \frac{{1 \mp i}}{2}e^{ \pm \frac{\pi }{2}iz^2 } . $$ This result is valid for all $N\geq 0$ and $z$ on the universal covering of $\mathbb C \setminus \left\{ 0\right\}$. You can see that if we omit the second term (which is exponentially small when $ 0 < \pm \arg z < \frac{\pi}{2}$) then, with $N=0$, we get the false result $$ \operatorname{f}(ze^{ \mp \frac{\pi }{2}i} ) = \pm i\operatorname{f}(z). $$ See also http://dlmf.nist.gov/7.4.

In summary, use $$ \operatorname{f}(z) \sim \frac{1}{{\pi z}}\sum\limits_{m = 0}^\infty {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} $$ when $\left| {\arg z} \right| \le \frac{\pi }{4}$, and use $$ \operatorname{f}(z) \sim \frac{{1 \pm i}}{2}e^{ \pm \frac{\pi }{2}iz^2 } + \frac{1}{{\pi z}}\sum\limits_{m = 0}^\infty {( - 1)^m \left( {\frac{1}{2}} \right)_{2m} \frac{1}{{(\pi z^2 /2)^{2m} }}} $$ when $\frac{\pi }{4} < \pm \arg z \le \frac{{3\pi }}{4}$. Of couse by relaying on the symmetry relation http://dlmf.nist.gov/7.4.E8, we can assume $|\arg z|\leq \frac{\pi}{2}$.

$\endgroup$
1
  • $\begingroup$ Thanks Gary for the extensive write-up! The error bounds are very instructive and indeed say "stay away from the imaginary axis". I've just PR'ed a fix in SciPy, using the symmetry relations in S and C themselves, which turned out to be somewhat easier. $\endgroup$ Oct 28, 2021 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.