# Are integrals with coalescing saddles available in the core Mathematica or in standard packages?

I was having a look at the DLMF chapter on integrals with coalescing saddles, and there are a few things that I'd like to try out with the functions it describes, including e.g. the Pearcey integral, $$\Psi_{2}\left(\mathbf{x}\right)=P(x_{2},x_{1})=\int_{-\infty}^{\infty} \exp\left(\mathrm{i}(t^{4}+x_{2}t^{2}+x_{1}t)\right)\mathrm{d}t,$$ or the elliptic and hyperbolic umbilic catastrophe integrals, $$\Psi^{(\mathrm{E})}\left(x,y,z\right)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \exp\left(i\left(s^3-3⁢s⁢t^2+z⁢(s^2+t^2)+y⁢t+x⁢s \right)\right)\mathrm{d}s\,\mathrm{d}t$$ and $$\Psi^{(\mathrm{H})}\left(x,y,z\right)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \exp\left(i\left(s^3+t^3+zst+yt+xs,\right)% \right)\mathrm{d}s\,\mathrm{d}t,$$ and I was wondering whether they are implemented in the core Wolfram language and in Mathematica (which doesn't appear to be the case) or in reasonably standard packages. What good resources are available for accessing these functions in Mathematica? I'm interested in pushing these things fairly hard, so I would like to e.g. be able to calculate them at arbitrary complex-valued parameters, as well as fairly deep into the asymptotic, highly-oscillatory region.

• Only for some special cases. E.g. the first integral equals $\frac{1}{4} (-1)^{3/8} e^{-\frac{i}{8}} \pi H_{-\frac{1}{4}}^{(1)}\left(\frac{1}{8}\right)$ if $x_2=1,x_1=0$. The one equals $\frac{\sqrt[8]{-1} \pi e^{-\frac{1}{8} \left(i \text{x2}^2\right)} \left(\left| \text{x2}\right| J_{-\frac{1}{4}}\left(\frac{\text{x2}^2}{8}\right)+(-1)^{3/4} \text{x2} J_{\frac{1}{4}}\left(\frac{\text{x2}^2}{8}\right)\right)}{2 \sqrt{2} \sqrt{\left| \text{x2}\right| }}$ if $x2\in \mathbb R ,x_1=0$. Jan 12, 2018 at 18:06
• @user64494 Sure - though the reduction of order is clear from the integral itself, which can be written as $\int_0^\infty e^{i(u^2+x_2u)}/\sqrt{u}\, \mathrm du$ when $x_1=0$, lowering the order of the catastrophe to $K=0$. The primary interest is for the parameter choices where the functions don't reduce to a lower order of catastrophe, though. Jan 12, 2018 at 18:19

I've done some preliminary work on numerically evaluating these saddle-point integrals a few years back. The basic idea, as you may have surmised, is to convert these into contour integrals that are rapidly convergent. Unfortunately, my attempts for the umbilic catastrophe integrals and the swallowtail catastrophe integral have been a bit hit-and-miss, so I'm going to defer posting about those for later.

In the case of Pearcey's integral, however, Connor and Farrelly give an expression that happens to be convenient to implement in Mathematica:

SetAttributes[pearcey, Listable];
pearcey[x_?NumericQ, y_?NumericQ] /; Precision[{x, y}] < ∞ :=
With[{prec = Precision[{x, y}]}, 2 Exp[I π/8]
NIntegrate[Exp[y t^2 Exp[3 I π/4] - t^4] Cosh[x t Exp[5 I π/8]],
{t, 0, ∞}, Method -> "Trapezoidal", WorkingPrecision -> prec]]


Here's how to reproduce the plots in the DLMF:

(* DLMF color scheme from http://dlmf.nist.gov/help/vrml/aboutcolor *)
DLMFContinuousColorPhase[u_, s_: 1, b_: 1] := Module[{rgb},
rgb = Clip[{1, -1, -1} Abs[{8, 4, 8} Mod[u/(2 π), 1] - {9, 3, 11}/2] +
{-3, 3, 5}/2, {0, 1}];
Apply[RGBColor, b (1 + s (rgb - 1))]]

{DensityPlot[Abs[pearcey[x, y]], {x, -10, 10}, {y, -46/5, 2},
ColorFunction -> (Hue[2 (1 - #)/3] &), PlotPoints -> 75],
DensityPlot[Arg[pearcey[x, y]], {x, -10, 10}, {y, -46/5, 2},
ColorFunction -> DLMFContinuousColorPhase, ColorFunctionScaling -> False,
Exclusions -> None, PlotPoints -> 75]} // GraphicsRow


These evaluations are a bit slow, however; I have yet to come up with a more efficient implementation.

• I'm a bit wary of the numerical integrations, particularly as you go to the $y\ll -1$ region and the integrals get very oscillatory - I've found the simpler NIntegrate implementations to suddenly go belly-up. Can you comment on that regime? Mar 4, 2018 at 15:19
• "I'm a bit wary of the numerical integrations" - me too; that's why I felt reluctant to share my implementations for the umbilic and swallowtail integrals (altho I have managed to reproduce all the plots in the DLMF). I haven't had a problem with the Connor-Farrelly formula, tho. If you manage to break it, I'd be interested in looking into it. Mar 4, 2018 at 15:25
• Well, I'm running those same plots down to $y=-18$ and it's already complaining about NIntegrate::ncvi. I suspect that any efficient implementation in that regime will have to use explicit asymptotics or at least alter its integration contour so that it looks more similar to a steepest-descent contour. Mar 4, 2018 at 15:32
• "alter its integration contour so that it looks more similar to a steepest-descent contour." - that actually is what I was trying to do for the other integrals, with limited success. I'll look into the literature again for Pearcey, since this is the most well-studied of the lot. Mar 4, 2018 at 15:36
• Update: yeah, this didn't work. Mar 4, 2018 at 16:01