17
$\begingroup$

I am looking for an accurate algorithm to calculate the error function

$$\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\ dt$$

I have tried using the following formula,

math97 second question example

(Handbook of Mathematical Functions, formula 7.1.26), but the results are not accurate enough for the application.

$\endgroup$
  • 1
    $\begingroup$ You may want to take a look at python's code.google.com/p/mpmath or other libraries that advertise a "multiple precision" feature. Also, this may be a better question for stack overflow instead, since it's more of a computer science thing. $\endgroup$ – Jon Bringhurst Jul 20 '10 at 20:26
  • $\begingroup$ @Jon: Nope, I'm not interested in a library, there is no such library for the language I'm writing in (yet). I need the mathematical algorithm. $\endgroup$ – badp Jul 20 '10 at 20:49
  • $\begingroup$ Have you tried numerical integration? Gaussian Quadrature is an accurate technique $\endgroup$ – Digital Gal Aug 28 '10 at 1:25
  • $\begingroup$ GQ is nice, but with (a number of) efficient methods for computing $\mathrm{erf}$ already known, I don't see the point. $\endgroup$ – J. M. is a poor mathematician Aug 29 '10 at 23:07
  • $\begingroup$ @Evgeny this was probably python but I have honestly no idea what I was up to seven years ago $\endgroup$ – badp Nov 9 '17 at 8:54
16
$\begingroup$

I am assuming that you need the error function only for real values. For complex arguments there are other approaches, more complicated than what I will be suggesting.

If you're going the Taylor series route, the best series to use is formula 7.1.6 in Abramowitz and Stegun. It is not as prone to subtractive cancellation as the series derived from integrating the power series for $\exp(-x^2)$. This is good only for "small" arguments. For large arguments, you can use either the asymptotic series or the continued fraction representations.

Otherwise, may I direct you to these papers by S. Winitzki that give nice approximations to the error function.


(added on 5/4/2011)

I wrote about the computation of the (complementary) error function (couched in different notation) in this answer to a CV question.

$\endgroup$
  • $\begingroup$ Assumption correct. :) $\endgroup$ – badp Jul 30 '10 at 20:02
  • $\begingroup$ +1 for the Winitzki reference; I've seen that 2nd paper before + it's a nice one. $\endgroup$ – Jason S Aug 6 '10 at 13:05
  • $\begingroup$ The only thing infuriating for me about the Winitzki paper is that apart from the Hermite-Pade ansatzes (or however you pluralize these things :P), he gives no indication on how he arrived at these simple approximants. $\endgroup$ – J. M. is a poor mathematician Aug 6 '10 at 13:10
  • $\begingroup$ Web Archive of Winitzki Paper: web.archive.org/web/20100601094000/https://…, it´s Ansätze $\endgroup$ – mondlos Jan 25 at 13:45
6
$\begingroup$

You can use a Taylor polynomial of sufficient degree to guarantee the accuracy that you need. (The Taylor series for erf(x) is given on the Wikipedia page to which you linked.) The Lagrange Remainder term can be used to bound the error in the Taylor series approximation.

$\endgroup$
4
$\begingroup$

Here's a link to the boost c++ math library documentation. They use their implementation of the incomplete gamma function, which in turn uses a mixed approach depending on the argument. It's all fairly well documented should you care to duplicate their method. And it looks like their error is within a few multiples of the machine epsilon.

Other than that, I would try the Taylor series. Numerical approximation might lead to a larger error term than the analytic one though, and it will only be valid in a neighborhood of 0. For larger values you could use the asymptotic series.

Another idea would be to restrict the domain to a closed interval. If you size it properly, then the function will appear constant with respect to your machine precision outside of this interval. Once you have a compact domain, you can know exactly how many Taylor terms you need, or you can use other types of spline interpolation. Chebyshev polynomials come to mind.

As for the problem that the language your writing in has no such library already: for me that is probably not as big of a deal as you think. Most languages seem to have a way to link in C functions, and if that is the case, then there is an open source implementation somewhere out there.

$\endgroup$
  • $\begingroup$ The naïve (alternating) Maclaurin series is not really that numerically sound; I had already mentioned in my answer the modified series that has much better properties for computing $\mathrm{erf}(x)$ near the origin. The (Laplace) continued fraction tends to be slightly easier to handle than the asymptotic series for medium-to-large arguments. $\endgroup$ – J. M. is a poor mathematician Sep 1 '11 at 10:34
  • $\begingroup$ If you're going for approximations of fixed degree near the origin, constructing a Padé approximant is slightly better than using a truncated Maclaurin series. $\endgroup$ – J. M. is a poor mathematician Sep 1 '11 at 10:35
  • $\begingroup$ I'll agree with that assessment. I thought about mentioning the numerical instability, but the post was already long. I think the best bet is to use a hybrid approach depending on the size of the argument. That way you can make an appropriate trade off of precision versus speed. Which method you use for which intervals is down to experimentation. $\endgroup$ – Tim Seguine Sep 1 '11 at 10:51
  • $\begingroup$ A lot of this comes down to the desired accuracy and how fast it must be... I think Chebyshev interpolation is worth looking into in any case $\endgroup$ – Tim Seguine Sep 1 '11 at 10:56
4
$\begingroup$

A simple way of computing error function is to use Kummer's equation of the form of,

$$ M(a,b,z)=\sum_{s=0}^{\infty} \frac{(a)_s}{(b)_s s!}z^s=1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)2!}z^2+...,\quad z \in \mathbb{C} $$

and

$$ M(a,b,z)=\sum_{s=0}^{\infty}\frac{(a)_s}{\Gamma{(b+s)}s!}z^s $$

and

$$ \text{erf}(z)=\frac{2z}{\sqrt{\pi}}M(.5,1.5,-z^2)=\frac{2z}{\sqrt{\pi}}e^{-z^2}M(1,1.5,z^2) $$

Here is the R code,

f<-function(a,b,z,maxt=5){
  s=1:maxt
  a=a+s
  b=b+s
  ss=1
  for(i in s){
    mt=prod(a[1:i]/b[1:i])
    ss=ss+mt *z^(2*i)/factorial(i)
  }
  ss=2*z/sqrt(pi)*exp(-z^2)*ss
  return(ss)
}

here is the plot (for real z) enter image description here

$\endgroup$
0
$\begingroup$

I'm surprised that nobody proposed the obvious option: Numerical quadrature. One of the simplest methods would be the composed trapezoidal rule, where you take equally spaced points $ 0=x_0, x_1, \cdots, x_n = x$, with spacing $h$, and obtain the estimate $$ \mathrm{erf} (x) \approx \frac{h}{\sqrt{\pi}}\left(1+e^{-x^2} + 2 \sum_{i=1}^{n-1} e^{-x_i^2}\right) $$

You can use the error formula to decide which $n$ to use for a particular target precision.

$\endgroup$
  • $\begingroup$ Could you also explore the computational complexity to reach the same accuracy as the other methods mentioned, starting with the error $10^{-7}$ method given in the question? $\endgroup$ – LutzL May 2 at 8:48
  • $\begingroup$ For this specific method, the error function would be given by $$|\varepsilon(x)| = \left| \frac{x^3}{12n^2} f''(\xi) \right| \leq \frac{x^3}{6n^2}$$ So, for a given $x$, you can choose $n$ so that any prescribed accuracy is reached. Performing Richardson extrapolation could accelerate the process (Romberg's method). I'm not sure this was your question though.@Lutz $\endgroup$ – PierreCarre May 2 at 10:28
  • $\begingroup$ So to get an error of $10^{-7}$ at $x=10$ you need $n\simeq 4\cdot 10^4$ for the trapzoidal rule or possibly $n\sim 10^2$ for the Simpson rule. Going down further will probably increase the constants too much due to insufficient sampling density. This is still a lot more effort than the dedicated approximation formulas, more so when $x$ is large. Any true improvement, variable step size, re-parametrization of the integral to get a more linear integrand etc. will lead to methods mentioned in the other answers. So your "obvious option" is not so obvious after all. $\endgroup$ – LutzL May 2 at 10:49
  • $\begingroup$ @LutzL I see your point and kind of agree, but I still think it is obvious, in the sense that is a very natural first approach that will work for a wide range of problems. Considering that for an error of $10^{-7}$ you don't need to compute for any $x>4$ (you can set the value to $1$) and using Romberg's method I expect this solution to be competitive, while having the advantage of requiring no extra work should you decide to make any changes in the integrand. $\endgroup$ – PierreCarre May 2 at 11:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.