How to accurately calculate the error function $\operatorname{erf}(x)$ with a computer?

Question

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,

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

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

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)

Answer 5

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.

Answer 6

You could get good theoretical approximations of the error function using $$\text{erf}\left(x\right)\sim \text{sgn}(x)\sqrt{1-\exp\Big(-\frac 4 {\pi}\,\frac{1+P_n(x^2)}{1+Q_n(x^2)}\,x^2 \Big)}$$ For example $$P_1(x^2)=\frac{10-\pi ^2}{5 (\pi -3) \pi }x^2\qquad Q_1(x^2)=\frac{120-60 \pi +7 \pi ^2}{15 (\pi -3) \pi }$$ while $$P_2(x^2)=\frac{\left(105840-110880 \pi +37800 \pi ^2-4260 \pi ^3+69 \pi ^4-17 \pi ^5\right) }{3 \pi \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^2+$$ $$\frac{-2116800+1270080 \pi ^2-504000 \pi ^3+48510 \pi ^4+503 \pi ^6}{315 \pi ^2 \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^4$$

$$Q_2(x^2)=\frac{60480-70560 \pi +27720 \pi ^2-3600 \pi ^3-143 \pi ^4+43 \pi ^5}{\pi \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^2+$$ $$\frac{-6350400+8467200 \pi -4515840 \pi ^2+1192800 \pi ^3-145320 \pi ^4+2380 \pi ^5+793 \pi ^6}{105 \pi ^2 \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^4$$ and so on.

The first one gives a maximum absolute error of $0.000150$ while the second gives a maximum absolute error of $0.000012$.

With regard to the infinite norms, they are respectively $3.04\times 10^{-8}$ and $1.20\times 10^{-10}$.