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]

Any suggestions for how this integral could be numerically evaluated?

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]

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.

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,

If the limit exists the improper integral is said to converge. Otherwise the integral is said to diverge.

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).

Figure 4.4.7
We do not know how to find the indefinite integrals in this example. Nevertheless the answer is 0 because on changing variables both limits of integration become the same.

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!

You could also write the integral in terms of special functions like erf.

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.

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\omega^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].

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

By cyclic, do you mean periodic?

 

[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.

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.

 

[maverick280857]

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.

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

 

[JJacquelin]

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

 

[jackmell]

Any suggestions for how this integral could be numerically evaluated?[/B]


Thanks in advance!

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?

Clear[t]
\[Omega] = 1; 
\[Sigma] = 1; 
myf[(t_)?NumericQ] := NIntegrate[
   Exp[I*\[Omega]*u]*Erf[u/(Sqrt[2]*\[Sigma])], 
   {u, -Infinity, t}]
mytable = Table[{t, myf[t]}, {t, 0, 20, 0.25}]; 
ListPlot[{({#1[[1]], Re[#1[[2]]]} & ) /@ mytable, 
   ({#1[[1]], Im[#1[[2]]]} & ) /@ mytable}, 
  Joined -> True]

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)?

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]

 

[JJacquelin]

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