import numpy as np
from scipy.special import erf, gamma
import matplotlib.pyplot as plt

The error function

What is the error function?

The error function is defined as

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

$\mathrm{erf}(x): \mathbb R \to \mathbb R$, $t \in \mathbb R$.

Let's first try to understand this by plotting the the integrand, $e^{-t^2}$. The error function is then the area under the curve from the vertical axis up to $x$.

x = 0.7

plt.axvline(x,  ls='--')
t = np.linspace(0,4,1000)
plt.plot(t, np.exp(-t**2))

# shade under the curve
# The area shaded is the error(x)
integrated_region = np.linspace(0,x,1000)
plt.fill_between(integrated_region,np.exp(-integrated_region**2),alpha=0.5)

<matplotlib.collections.FillBetweenPolyCollection at 0x11152f770>

Now let's try to approximate the value of error(x) with a Riemmann sum (simple numerical integral)

x = 0.7

t = np.linspace(0,x,100000)
dt = np.diff(t)[0]

erf_approx = 2/np.sqrt(np.pi) * np.sum(np.exp(-t**2)*dt)
print(erf_approx, erf(x))

0.6778075626862514 0.6778011938374184

Pretty close!

How does the error function itself looks like?

t = np.linspace(-5,5,1000)
plt.plot(t,erf(t))

[<matplotlib.lines.Line2D at 0x111549d00>]

Mental test: when $x$ becomes larger and larger, the area under the curve doesn't increase much because the integrand $\to 0$.

Approximation with a series

$$ \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty \frac{(-1)^n x^{2n + 1}}{n! (2n+1)} $$

The error function as the limit of a series

One can develop the error function as a power series. Alternatively, when computing a problem with series, one can perhaps recognize that the terms of the series can be rewritten in terms of the error function.

Here, we test numerically how good the error function series expansion is.

big = 10

n = np.arange(0,big,1,dtype=int)
x = 0.7
nfactorial = gamma(n+1)

erf_series = 2/np.sqrt(np.pi)*np.sum( ((-1)**n * x**(2*n+1) )/(nfactorial*(2*n+1))  )

erf_series, erf(x)

(np.float64(0.6778011938294721), np.float64(0.6778011938374184))

What if the argument is complex?

One can extend the definition for $z \in \mathbb C$ such that $$ \mathrm{erf}(z) = \frac{2}{\sqrt{\pi}}\int_0^z e^{-t^2} dt $$

Let's visualize this.

x = np.linspace(-10,10,1000)
y = np.linspace(-10,10,1000)
x, y = np.meshgrid(x,y,indexing="ij")
z = x + 1j*y
erf_eval = erf(z)

This is the real part of the complex error function. Note how in the horizontal axis, we get the same behavior we plotted previously for erf(x).

plt.pcolormesh(np.real(z),np.imag(z),np.real(erf_eval),vmin=-10,vmax=10,cmap="seismic_r")
plt.colorbar()
plt.grid()

Now this is the imaginary part. Note that on the horizontal axis, we get zeros, confirming that the real erf(x) gives out only reals.

plt.pcolormesh(np.real(z),np.imag(z),np.imag(erf_eval),vmin=-10,vmax=10,cmap="seismic_r")
plt.colorbar()
plt.grid()

From the definition of erf(x), one can see that for a complex argument $z$

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

and a change of variables inside of the integral such that the integration limit becomes real but then the integrand must become complex, there is a relationship between the error function and integrals involving $\cos u^2$ and $\sin u^2$ for real $u$ (after the change of variable, the exponential becomes complex). Those are the Fresnel integrals