Assigning an analytical approximation to the error function erf(x)

Working with some iterative integral equations, I have Gaussian density functions involved therein. Integrating the Gaussian function I obtain the error function. When I take the second integration, Mathematica is not able to provide an analytical expression. Because apparently there is no closed-form expression for the integral of the error function.

Question: Is it possible the tell Mathematica to consider the error function as some given exact expression which is analytically integrable?

For example

$$\mathrm{erf}(x)=\frac{1}{(1+a_1x+a_2x^2+a_3x^3+a_4x^4)^4}$$

for some appropriate constants $a_1$, $a_2$, $a_3$, $a_4$?

• Yes, it's possible. Take a look at a closely related question Numerical underflow for a scaled error function. I guess therein you'll find answer to your question. – Artes Oct 22 '13 at 23:34
• hi @Artes thank you very much for your comment. I am looking at it right now. – Seyhmus Güngören Oct 22 '13 at 23:36
• I mean you can find only closed-form approximations, not an analytic closed-form formula. – Artes Oct 22 '13 at 23:38
• yes but my question is that, whenever a result comes as $erf$ can I tell mathematica to use a predefined function instead of it automatically. Can I interrupt the general understanding of mathematica from the error function? – Seyhmus Güngören Oct 22 '13 at 23:42
• And you can't adapt Integrate[Erf[t],{t,0,x}] -> x*Erf[x] + (Exp[-x^2] - 1)/Sqrt[Pi]?? – Ray Koopman Oct 24 '13 at 6:03

You could unprotect Erf and redefine it as you like it.

Unprotect[Erf];
(* Here you define the constants and your erf approximation *)
a1=...; etc.
Erf[x_] := 1/(1 + a1 x + a2 x^2 + a3 x^3 + a4 x^4)^4;
(* Here you perform your calculations with the approximate value of erf *)
myComplexCalculation = Erf[x]/2
(* When you are done, clear your definition and reprotect Erf *)
Clear[Erf];
Protect[Erf];

This can be dangerous but, as long as you know what you are doing, it will work.

A more local workaround is to enclose your 'approximate' calculation within a Block that blocks the definition of Erf:

Block[ {Erf},
(* Here you define the constants and your erf approximation *)
a1=...; etc.
Erf[x_] := 1/(1 + a1 x + a2 x^2 + a3 x^3 + a4 x^4)^4;

(* Here you perform your calculations with the approximate value of erf *)
myComplexCalculation = Erf[x]/2
]

Is this what you were looking for or was it the actual expression for the approximation?

As a more general approach to temporarily modifying built-in functions within the walls of a scoping construct, you might also want to have a look at the undocumented function InheritedBlock, described in this post on Stack Overflow. It appears that it will copy all definitions associated with the symbol you want to block, so that you can change, within the block only the part you need to. Basically, while

Block[ {builtinFun},
code redefining builtinFun from scratch
and using the new builtinFun
]

would give you a completely clean symbol for builtinFun to code from scratch, using

Internal`InheritedBlock[ {builtinFun},
code redefining parts of the inherited builtinFun
and using the modified builtinFun
]

would give you a completely modifiable copy of builtinFun, so that you can modify only the parts you need to modify.

• It seems like this was what I was looking for. I will check and let you know. – Seyhmus Güngören Oct 23 '13 at 10:35

Altrnate to redefining Erf you might define your own approximate probability distribution function and work with that , eg..

cvals = {a1 -> .278393, a2 -> .230389, a3 -> .000972, a4 -> .078108}

(* erf approximation from wikipedia , take for what its worth *)
approxerf[x_] :=
Sign[x] (
1 - 1/(1 + a1 Abs[x] + a2 x^2 + a3 Abs[x]^3 + a4 x^4)^4) /. cvals

approxpdf[mean_, sd_] = (512/sd) Piecewise[{
{(
4 g a2 + 4 a4 g^3  + Sqrt 2 a1 +
3 Sqrt a3 (g^2)  )/((4 + 2 Sqrt a1 g + 2 a2 g^2 +
Sqrt a3 g^3 + a4 g^4)^5)  /. g -> (mean - x)/sd,
( mean - x)/sd >= 0},
{(
4 g a2 + 4 a4 g^3  + Sqrt 2 a1 +
3 Sqrt a3 (g^2)  )/((4 + 2 Sqrt a1 g + 2 a2 g^2 +
Sqrt a3 g^3 + a4 g^4)^5)  /. g -> -((mean - x)/sd),
True}}] /. cvals;

mynormal[mean_, sd_] =
ProbabilityDistribution[
approxpdf[mean, sd], {x, -Infinity, Infinity}];

Plot[ { CDF[mynormal[2, 5], x] ,
CDF[NormalDistribution[2, 5], x]}  , {x, -20, 20}]

This should be equivalent - more a matter of semantics. Do you want an approximate solution to an exact problem or an exact solution to an approximate problem.

Note this performs pretty slowly..but I assume the other approach has the same issue as that analytic expression isn't particularly nice to integrate.