# Numerical approximation of a complex integral with a nested exponential

I've been working on a maths problem as part of my Physics PhD; but have been stumped by the following integral. All I need to know is a numeric approximation to the integral (along with an estimation of the error); but if an analytic solution exists, that would be fantastic.

$$\int_{-\infty}^{\infty} G(x,\sigma) e^{(-i(\delta x - \nu \exp(-\kappa x)/\kappa))}dx$$

where $G(x,\sigma)$ is a Gaussian distribution in $x$ with standard deviation $\sigma$ and mean zero; and where $\delta$, $\nu$ and $\kappa$ are real constants.

I have tried several approaches, including method of steepest descent/stationary phase and assymptotic expansions in $x$; but every method apart from a brute force ODE solver or polygon approach [both of which are extremely slow, such that I have not actually allowed them to complete] have proven to be very unreliable. I have considered contour integration also, but since this integral does not have any poles, I am not sure how helpful this would be.

I have not found there to be any simplifying relationships between the real constants, but here are some typical values:

$$\begin{eqnarray} \delta&\approx&5\times10^6\\ \nu&\approx&5\times10^6\textrm{ And slightly different from \delta}\\ \kappa&\epsilon& \mathbb{R^+}\\ \sigma&\gg&\frac{1}{\kappa} \end{eqnarray}$$

Anyway, if anyone has any ideas on how such an integral might be solved, I would be ever so grateful.

-
I am assuming $\nu$ and $\kappa$ are constants. And what is that $\Delta x$? – Fixed Point Oct 23 '12 at 5:21
$\Delta$ is just another constant as well; sorry. Added to the description above (and replaced with $\delta$ for clarity). – Matthew Oct 23 '12 at 5:28
Are there any values of the constants for which you have been able to compute the integral numerically? And what are typical values for those four constants? – Fixed Point Oct 23 '12 at 5:49
I have added some typical values into the question :). – Matthew Oct 23 '12 at 6:02

The reason you can't use the stationary phase technique here is that there's no point of stationary phase – the derivative of the phase

\begin{align} \phi(x)&=\delta x-\nu\exp(-\kappa x)/\kappa\;, \\ \phi'(x)&=\delta+\nu\exp(-\kappa x)\;, \end{align}

is everywhere positive. You can turn that to your advantage though, since it allows you to change variables to $\phi$, yielding

$$I=\int_{-\infty}^\infty \frac{G(x(\phi),\sigma)}{\phi'}\exp(-\mathrm i\phi)\,\mathrm d\phi\;.$$

This is a Fourier component of a slowly varying function. For positive $x$ we have $\phi\sim\delta x$, whereas for negative $x$ we have $\phi\sim-\nu\exp(-\kappa x)/\kappa$ and thus $x/\sigma\sim\log(-\phi)/(\kappa\sigma)$, so the integral should decay exponentially with $\delta$ and $\kappa\sigma$. You can get upper bounds for its magnitude by integrating by parts; the boundary terms vanish because of the Gaussian, so you get

$$|I|=\left|\int_{-\infty}^\infty \frac{G(x(\phi),\sigma)}{\phi'}\exp(-\mathrm i\phi)\,\mathrm d\phi\right|\le \int_{-\infty}^\infty\left|\frac{\mathrm d^n}{\mathrm d\phi^n}\frac{G(x(\phi),\sigma)}{\phi'}\right|\,\mathrm d\phi$$

for all $n\in\mathbb N$. Since each derivative yields a large inverse factor, this should lead to very small bounds on the magnitude of the integral.

These bounds are easier to evaluate if we transform back to $x$. For instance, for $n=0$ we have

$$|I|\le\int_{-\infty}^\infty G(x,\sigma)\,\mathrm dx=1\;,$$

for $n=1$, using $|\phi'|\lt\delta$ and $|\phi''/\phi'|\lt\kappa$,

\begin{align} |I|&\le\int_{-\infty}^\infty\left|\frac{G'(x(\phi),\sigma)}{\phi'^2}-\frac{G(x(\phi),\sigma)\phi''}{\phi'^3}\right|\,\mathrm d\phi \\ &=\int_{-\infty}^\infty\left|\frac{G'(x,\sigma)}{\phi'}-\frac{G(x,\sigma)\phi''}{\phi'^2}\right|\,\mathrm dx \\ &\le\frac1\delta\int_{-\infty}^\infty\left(|G'(x,\sigma)|+\kappa|G(x,\sigma)|\right)\,\mathrm dx \\ &=\frac1\delta\left(2G(0,\sigma)+\kappa\right) \\ &=\frac1\delta\left(\frac2{\sqrt{2\pi}\sigma}+\kappa\right)\;, \end{align}

and for $n=2$, using $|\phi'''/\phi'|\le\kappa^2$,

\begin{align} |I|&\le\int_{-\infty}^\infty\left|\frac{G''}{\phi'^3}-3\frac{G'\phi''}{\phi'^4}+3\frac{G\phi''^2}{\phi'^5}-\frac{G\phi'''}{\phi'^4}\right|\,\mathrm d\phi \\ &=\int_{-\infty}^\infty\left|\frac{G''}{\phi'^2}-3\frac{G'\phi''}{\phi'^3}+3\frac{G\phi''^2}{\phi'^4}-\frac{G\phi'''}{\phi'^3}\right|\,\mathrm dx \\ &\le\frac1{\delta^2}\int_{-\infty}^\infty\left(|G''|+3\kappa |G'|+3\kappa^2|G|+\kappa^2|G|\right)\,\mathrm dx \\ &\le\frac1{\delta^2}\left(\frac2{\sigma^2}+\frac{6\kappa}{\sqrt{2\pi}\sigma}+4\kappa^2\right)\;. \end{align}

We can go on like this until the cows come home, that is, until the numerical factors from the derivatives begin to overwhelm the factors of $\kappa/\delta$, so if $\kappa\ll\delta$, the integral is for all intents and purposes zero. If $\kappa\sim\delta$, you may have to work a bit harder to get better bounds.

-
Thanks very much joriki :). I've gone through your maths, and agree with everything. I had not thought to transform into $\phi$, so thank you. I have yet to check if I can use it for my purposes, though. If it turns out that I cannot compute the integral directly, what would be lovely is a positive (non-zero) lower bound... but that might be impossible to find. I'll give it some more thought... but if you already see such a thing, I would love to know :). Thanks again! ^^. – Matthew Oct 23 '12 at 23:43
Oh... and thanks for pointing out my mistake with the stationary phase method. I mistakenly assumed I need only expand around $\phi(x) = 0$. Whoops :). – Matthew Oct 23 '12 at 23:51
+1: Nice job! – AD. Oct 24 '12 at 6:48
As I've been playing with this integral further, I discovered that the assumption used to generate the expressions of the integrals (e.g. $|\phi^\prime| < \delta$) are not quite right (it should be $|\phi^\prime| > \delta$); or so it seems to me. This might actually be quite nice :). – Matthew Oct 25 '12 at 23:32
Note that the results following are correct; however :). – Matthew Oct 25 '12 at 23:39

Here are some steps in order to arrive at an estimate.

Break up the exponential using the Euler identity $\exp(it)=\cos t + i\sin t$ $$\exp(-i(\delta x -\nu\exp(-\kappa x)/\kappa) =\\ \cos(\delta x -\nu\exp(-\kappa x)/\kappa)-i\sin(\delta x -\nu\exp(-\kappa x)/\kappa)$$ Next, we get two real integrals which we estimate separately, now the hard work begins which I must leave to you. However, note that for large $x$ $$\cos(\delta x -\nu\exp(-\kappa x)/\kappa)-i\sin(\delta x -\nu\exp(-\kappa x)/\kappa)\sim \\ \cos(\delta x)-i\sin(\delta x)$$

-
Thanks AD :). The problem is that I can guarantee that the exponential term dominates in at least part of the integration domain; and so I cannot make that kind of approximation. Furthermore, because the exponential blows up, the oscillation becomes very rapid, which almost certainly makes it nearly identical to zero in this region. The trick is to know where to truncate it to get a known error. Anyway... thank you again! :). – Matthew Oct 23 '12 at 22:58
@Matthew Yes of course! That was not what I meant, this is for $x$ close to $+\infty$, when $x$ is close to $-\infty$ you will as you say have a rapid oscillation. Also, this is just to get started. – AD. Oct 24 '12 at 6:45
As you clearly said. Thank you again for your kind assistance. :). – Matthew Oct 24 '12 at 23:01

Good answer by AD, it would have been a nice complement to my idea. I was going to suggest (now that you have two real integrals on an infinite domain) you can just use a transformation of variables, a map from $(-\infty,\infty)$ to $(-1,1)$ will change the interval of integration to a finite one (your integrand will change too of course) and then just use any standard technique. I was going to suggest then to just play around with the two real (transformed) integrands to figure out which technique is best for them and then just use one of the many canned routines available out there (unless you want to have a chapter of your thesis on how you came up with a new/novel way to integrate functions of this kind ;-) ). Because of the sines and cosines we suspect oscillations, I would have suggested looking at Gaussian quadrature with chebyshev nodes.

This was all before I started playing with the two real integrands myself. Holy shit these things are oscillatory especially if I use numbers you suggested. The graph is going across orders of magnitudes wildly fast. For example, I used $\sigma = 3700, \delta = 5000000, \nu=50001000, \kappa=10$. Mathematica couldn't even plot it. Well it spit out a picture but it was complete garbage. It cannot find an anti-derivative (there goes that) but numerically it instantly spit out an answer on the order of $10^{-14}$ so it thinks that the integral is zero which it really does look like it might be. Then I forced the numerical integrator to use a working precision of 10 digits (it carries at least 10 digits in all of its computations) and give me an answer that is at least accurate to five digits (these are...ehmmm...usually very modest goals for a world class software like mathematica) and here is what it spat out at me

"NIntegrate::mtdfb: Numerical integration with LevinRule failed. The integration continues with Method -> GaussKronrodRule. >>

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-0.4464083980}. NIntegrate obtained -5.478908926*10^-7+1.469367460*10^-6 I and 0.00001382317328. for the integral and error estimates. >>"

Just look at the integral estimate and its error in the last line. The error is like two orders of magnitude HIGHER than the error estimates. This output is from trying to integrate the original complex integral you presented so the answer is a complex number.

Now the only things I can think of is, first take comfort in the fact that this really is a hard problem. Don't feel bad for having trouble with this. If mathematica is having trouble, then it really is a tricky problem. I was going to suggest maybe some simplification like taylor expansion and then truncation but taking derivatives will be horrendous. Maybe restricting the interval of integration because of $G(x,\sigma)$ the integrands converge rather fast to zero as you go out to infinity in both directions. After chopping the interval using fourier expansion might also be interesting. It seems "natural" because of the oscillations.

Time to walk over to the math department and collaborate with the numerics people. Sorry couldn't help.

-
Thanks Prince Ali :). I went through a similar process as did you; trying Mathematica and Sage (with various algorithms); and this problem just didn't want to co-operate. Thanks again for your help :). – Matthew Oct 23 '12 at 22:55

Here's a different approach, using contour integration. You wrote that you're not sure how helpful this would be because the integrand doesn't have any poles; but contour integration isn't only useful for applying the residue theorem; it's also useful for shifting the contour of integration, and especially when the integrand oscillates, shifting the contour can transform that into exponential damping.

We can't close the contour with a semicircle at infinity because both the Gaussian and the exponential in the exponent would blow up in our face. But nevertheless we can turn some of the rapidly oscillating phase factors into rapidly decaying exponentials.

We have two distinct regions: For negative $x$ the exponent is dominated by the exponential, and for positive $x$ it's dominated by $\delta x$. Let's deal first with the exponential.

This blows up when $0\lt\kappa\def\im{\operatorname{Im}}\im x\lt\pi\bmod2\pi$ and is most damping at $\kappa\im x=-\pi/2$. Thus, a first step is to shift the contour to $x=t-\frac\pi{2\kappa}\mathrm i$. (More precisely, we can integrate around the strip $0\ge\im x\ge-\frac\pi{2\kappa}$, and since the short contributions at infinity vanish because of the Gaussian and the strip contains no poles, the original integral and the shifted integral must be equal to cancel.)

The result is

\begin{align} &\int_{-\infty}^\infty G\left(t-\frac\pi{2\kappa}\mathrm i,\sigma\right)\exp\left(-\mathrm i\left(\delta\left(t-\frac\pi{2\kappa}\mathrm i\right)-\nu\exp\left(-\kappa \left(t-\frac\pi{2\kappa}\mathrm i\right)\right)/\kappa\right)\right)\,\mathrm dt \\ =&\exp\left(-\frac12\left(\frac\pi{2\kappa\sigma}\right)^2-\frac{\pi\delta}{2\kappa}\right)\int_{-\infty}^\infty G\left(t,\sigma\right)\exp\left(\mathrm i\left(\frac\pi{2\kappa\sigma}-\delta\right) t-\nu\exp\left(-\kappa t\right)/\kappa\right)\,\mathrm dt\;. \end{align}

As might have been expected from the bounds in my other answer, we now have a large damping factor on the outside that falls off exponentially with $\delta/\kappa$, and the exponential in the exponent has also turned into a strong damping factor that cuts the integral off for negative $t$ at $|t|\sim1/\kappa\ll\sigma$.

But we still have an oscillatory phase factor damped only by the Gaussian for positive $t$, and though we can bound it by $1$ to get a good bound given by the factors before the integral, we'd have to integrate over many periods to get an exact value. But for positive $t$ we don't have to worry about the exponential blowing up, so we can not only shift the contour but actually rotate it by $\pi/4$ to get the most rapid decay from the phase factor without making the Gaussian blow up. This works because we can close the contour in an eighth-circle at infinity because we always have either the Gaussian or the phase factor damping the integrand exponentially. So use $t=\mathrm e^{-\mathrm i\pi/4}u$ to obtain

\begin{align} &\exp\left(-\frac12\left(\frac\pi{2\kappa\sigma}\right)^2-\frac{\pi\delta}{2\kappa}\right)\left(\int_{-\infty}^0 G\left(t,\sigma\right)\exp\left(\mathrm i\left(\frac\pi{2\kappa\sigma}-\delta\right) t-\nu\exp\left(-\kappa t\right)/\kappa\right)\,\mathrm dt+\right. \\&\left.\frac1{\sqrt{2\pi}\sigma}\int_0^\infty \exp\left(\frac12\mathrm i\left(\frac u\sigma\right)^2+e^{\mathrm i\pi/4}\left(\frac\pi{2\kappa\sigma}-\delta\right) u-\nu\exp\left(-\kappa\mathrm e^{-\mathrm i\pi/4}u\right)/\kappa\right)\,\mathrm du\right)\;. \end{align}

The Gaussian has become entirely oscillatory, and the exponential and the linear phase factor have become partly oscillatory, partly decaying. Since $\delta$ is very large, this part of the integral now also falls off very rapidly. You should be able to numerically integrate this part quite accurately in this form; in fact you can get an easy estimate by keeping only linear terms in $u$ over the range over which $\exp\left(-\delta u/\sqrt2\right)$ is appreciable:

$$\frac1{\sqrt{2\pi}\sigma}\int_0^\infty \exp\left(e^{\mathrm i\pi/4}\left(\frac\pi{2\kappa\sigma}-\delta\right) u-\nu/\kappa+\nu\mathrm e^{-\mathrm i\pi/4}u\right)\,\mathrm du=\frac{\exp\left(-\nu/\kappa-\mathrm i\pi/4\right)}{\sqrt{2\pi}\sigma(\pi/(2\kappa\sigma)-\delta-\mathrm i\nu)}\;.$$

-