Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Problem with numerical integration of error function

  1. Oct 6, 2011 #1
    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!
     
  2. jcsd
  3. Oct 6, 2011 #2

    lurflurf

    User Avatar
    Homework Helper

    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 tH 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.
     
    Last edited: Oct 7, 2011
  4. Oct 7, 2011 #3
    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.
     
  5. Oct 7, 2011 #4
    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.
     
  6. Oct 7, 2011 #5
    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\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].

    By cyclic, do you mean periodic?
     
  7. Oct 8, 2011 #6
    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 Files:

  8. Oct 8, 2011 #7
    Thanks JJacquelin, this seems interesting. I'll check it out.
     
  9. Oct 9, 2011 #8
    A little change in notations f(t,T) for clarification in attachment :
     

    Attached Files:

  10. Oct 9, 2011 #9
    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?

    Code (Text):

    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]
     
    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 Files:

    Last edited: Oct 9, 2011
  11. Oct 11, 2011 #10
    A more simple way in order to separate the non-convegent part of the integral :
     

    Attached Files:

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Problem with numerical integration of error function
Loading...