Problem with numerical integration of error function

maverick280857

Hi,

For a small problem I'm working on in quantum field theory, I have to numerically evaluate the integral

[tex]\int_{-\infty}^{t}du e^{i\omega u}erf\left(\frac{u}{\sqrt{2}\sigma}\right)[/tex]

where [itex]erf(z)[/itex] is the error function.

Now, I have to replace [itex]-\infty[/itex] by some large negative number, so I effectively end up evaluating

[tex]\int_{-t_H}^{t}du e^{i\omega u}erf\left(\frac{u}{\sqrt{2}\sigma}\right)[/tex]

If I use Matlab or Mathematica to evaluate this integral numerically, I get warnings indicating that the integrand is singular, or highly oscillatory. Secondly, the choice of [itex]t_H[/itex] seems critical, and since the integrand itself doesn't fall off asymptotically, its not clear to me how [itex]t_H[/itex] should be chosen in terms of [itex]\sigma[/itex].

Note that in the [itex]t_H, t \rightarrow \infty[/itex] limit, this is just the Fourier Transform of the error function, which is well defined for [itex]\omega \neq 0[/itex]. Incidentally, [itex]\omega \neq 0[/itex] is ensured in my physical problem.

Any suggestions for how this integral could be numerically evaluated?

Secondly, what is a good way to evaluate the error function of a complex number, i.e. erf(a + ib) where a and b are real, and i = sqrt(-1)?

Thanks in advance!

lurflurf

The Fourier transform of unity involves a delta function so your function does as well. It would also help to know the domain of your parameters.
Assuming σ>0 we can write
erf(u/(sqrt(2)σ)=-1+(1+erf(u/(sqrt(2)σ))
(1+erf(u/(sqrt(2)σ))->0 rapidly for t_H large

This separates the delta function part from the numerical integration. If t is large we need to do something similar for it. You could also write the integral in terms of special functions like erf.

JJacquelin

This integral (with -infinity lower bound) cannot be numerically evaluated because it is not convergent.
For u tending to -infinity, the erf function tends to -1. The function to be integrated is equivalent to exp(iwu) which primitive is exp(iwu)/iw = (-i cos(wu)+sin(wu))/w +C. So the integral is oscillating.
If the lower bound is a large negative number, but not -infinity, the integral can be evaluated. But in practice, numerical integration would involve some difficulties due to a too large number of periods. In that case, it is more convenient to integrate analytically the function (This is possible in terms of erf and exp functions) and then numerically evaluate the formula which has been obtained.

LaurieAG

That may be true for linear constructs but possibly not for cyclic quantum field constructs.

http://www.math.wisc.edu/~keisler/calc.html

In the second version HJ Keisler says,
The new version just says that there are finite limit improper integrals that converge and other integrals with infinite limits that diverge (and are indefinite).

This quote and figure was in relation to a cyclic function with a changing variable and a sub part with infinite limits. So divergent indefinite integrals can equal 0 if their infinite limits are equal and they are part of a higher level cycle but convergent integrals are only improper and equal to 0 when they have finite limits.

maverick280857

Thanks for the replies everyone!

The integral can be expressed in terms of error functions. Specifically it evaluates to

[tex]\left[-\pi\delta(\omega)-\frac{i e^{-\sigma^2\omega^2/2}}{\omega}\left(erf\left(\frac{t}{\sqrt{2}\sigma}\right)e^{(\sigma^2\o mega^2 + 2i\omega t)/2} + erfc\left(\frac{t-i\sigma^2\omega}{\sqrt{2}\sigma}\right) - 2\right)\right][/tex]

The problem is the numerical evaluation of the error function of a complex argument.

If I replace the lower limit ([itex]-\infty[/itex]) by some large negative number, [itex]-T_{0}[/itex], my results are sensitive to the value of [itex]T_{0}[/itex] relative to [itex]\sigma[/itex]. Also, both MATLAB as well as Mathematica return warnings when evaluating the integral numerically even for finite [itex]T_0[/itex].

By cyclic, do you mean periodic?

JJacquelin

In addition, I would like to propose a method in order to give an approximate of the integral.
The joint page was written in hurry without checking. It will require to carry out more checks.

Attached Thumbnails

 

maverick280857

Thanks JJacquelin, this seems interesting. I'll check it out.

JJacquelin

A little change in notations f(t,T) for clarification in attachment :

Attached Thumbnails

 

jackmell

Mathematica version 8 computes the integral but does report that it may not converge. However, the value reported for the Fourier transform agrees with the actual transform result. Here's the code I use and a plot of the real and imaginary parts of the integral evaluated for t from 0 to 20. Are these plots consistent with your expectations?

Other than using Mathematica or other CAS, evaluate it just like any other contour integral: since the integrand is entire, it's path-independent so make the path the easiest. For example, to evaluate:

[tex]Erf(1+i)=\frac{2}{\sqrt{\pi}} \int_0^{1+i} e^{-t^2}dt[/tex]

just parameterize a straight-line path from the origin to the point 1+i as t=t+it so that the integral becomes:

[tex]Erf(1+i)=\frac{2}{\sqrt{\pi}} (1+i)\int_0^1 e^{-(t+it)^2} dt[/tex]

Attached Thumbnails

 

JJacquelin

A more simple way in order to separate the non-convegent part of the integral :

Attached Thumbnails