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

Question

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

Answer 1

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.

Answer 2

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] 2 a1 + 
          3 Sqrt[2] a3 (g^2)  )/((4 + 2 Sqrt[2] a1 g + 2 a2 g^2 + 
          Sqrt[2] a3 g^3 + a4 g^4)^5)  /. g -> (mean - x)/sd,
        ( mean - x)/sd >= 0},
     {(  
        4 g a2 + 4 a4 g^3  + Sqrt[2] 2 a1 + 
         3 Sqrt[2] a3 (g^2)  )/((4 + 2 Sqrt[2] a1 g + 2 a2 g^2 + 
         Sqrt[2] 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.