Published May 5, 20214 minute read

This post documents the implementation of the error function erf(z)\mathrm{erf}(z) used in my complex function plotter. For small zz, I use the asympotic series by Abramowitz & Stegun. For large zz, I use a custom expansion around the 45° line.

Let ω:=eiπ/4\omega \defeq e^{i \pi/4}, and parameterize the complex plane via z=rω+siωz = r \omega + s \, i \omega. Consider the piecewise linear contour traveling from 00 to rωr \omega, and then from rωr \omega to rω+siωr \omega + s \, i \omega. Using this contour, we compute

π2erf(z)=0zex2dx=0re(tω)2ωdt+0se(rω+tiω)2iωdt=ω[C(r)iS(r)]+iωeir20se2rt+it2dt,\begin{aligned} \frac{\sqrt{\pi}}{2} \, \mathrm{erf}(z) &= \int_0^z e^{-x^2} \dd{x} \\ &= \int_0^r e^{-(t \omega)^2} \omega \dd{t} + \int_0^s e^{-(r \omega + t i \omega)^2} i \omega \dd{t} \\ &= \omega \, [C(r) - i \, S(r)] + i \omega e^{-ir^2} \int_0^s e^{2rt + it^2} \dd{t}, \end{aligned}

where S(x)S(x) and C(x)C(x) are the Fresnel integrals. From Wikipedia, we have the asymptotic expansion

C(x)iS(x)=π2ω[1+O(x4)]cosx2isinx22x[12x2i]C(x) - i \, S(x) = \frac{\sqrt{\pi}}{2} \overline{\omega} - [1 + \O(x^{-4})] \, \frac{\cos{x^2} - i \sin{x^2}}{2x} \, \left[\frac{1}{2x^2} - i\right]

for large positive xx. The second term is more annoying. For large a>0a > 0, consider the integral

Ia(x):=eaxxeateit2dt=axxea(ts)eit2dsdt.\begin{aligned} I_a(x) &\defeq e^{-ax} \int_{-\infty}^x e^{at} e^{it^2} \dd{t} \\ &= a \int_{-\infty}^x \int_x^\infty e^{a \, (t-s)} e^{it^2} \dd{s} \dd{t}. \end{aligned}Making the change of variables (u,v):=(a(st),xt)(u, v) \defeq (a \, (s-t), x-t), we obtainIa(x)=0eu0u/aei(xv)2dvdu.\begin{aligned} I_a(x) &= \int_0^\infty e^{-u} \int_0^{u/a} e^{i \, (x-v)^2} \dd{v} \dd{u}. \end{aligned}

The inner integral can be evaluated as a sum of four Fresnel integrals, and is therefore bounded by a constant. Then for large aa, it is fruitful to perform the Taylor expansion

0u/aei(xv)2dv=eix20u/ae2ixv+iv2dv=eix20u/a[12ixv+iv22x2v2+O(v3)]dv=eix2[uaixu2a2+(i2x2)u33a3+O(a4)].\begin{aligned} \int_0^{u/a} e^{i \, (x-v)^2} \dd{v} &= e^{ix^2} \int_0^{u/a} e^{-2ixv + iv^2} \dd{v} \\ &= e^{ix^2} \int_0^{u/a} [1 - 2ixv + iv^2 - 2x^2v^2 + \O(v^3)] \dd{v} \\ &= e^{ix^2} \, \left[\frac{u}{a} - \frac{ixu^2}{a^2} + \frac{(i - 2x^2) \, u^3}{3a^3} + \O(a^{-4})\right]. \end{aligned}

Evaluating the standard gamma integrals, we thus obtain

Ia(x)=eix2[a12ixa2+(2i4x2)a3+O(a4)].\begin{aligned} I_a(x) &= e^{-ix^2} \left[a^{-1} - 2ixa^{-2} + (2i - 4x^2) \, a^{-3} + \O(a^{-4})\right]. \end{aligned}

We were a bit handwavy with the justification for Fubini. Indeed, the above expansion works well for x>0x > 0, but fails catastrophically for x<0x < 0. From the definition of Ia(x)I_a(x), it is evident that small errors in the integral are exponentially amplified for negative xx. Empirically, however, we find that the approximation

0xeateit2dt=eaxIa(x)Ia(0)eax+ix2[a12ixa2+(2i4x2)a3]a12ia3\begin{aligned} \int_0^x e^{at} e^{it^2} \dd{t} &= e^{ax} \, I_a(x) - I_a(0) \\ &\approx e^{ax+ix^2} \left[a^{-1} - 2ixa^{-2} + (2i - 4x^2) \, a^{-3}\right] - a^{-1} - 2i a^{-3} \end{aligned}

works very well, and could probably be derived rigorously with more effort. Combining all the above, we have

erf(z)=2πω[C(r)iS(r)]+2πiωeir2[e2rsI2r(s)I2r(0)]12πcosr2isinr22r[12ri]+iω2πeir2[e2rs+is2(2r12isr2+(i+2s2)r3)2r1ir3],\begin{aligned} \mathrm{erf}(z) &= \frac{2}{\sqrt{\pi}} \omega \, [C(r) - i \, S(r)] + \frac{2}{\sqrt{\pi}} i \omega e^{-ir^2} [e^{2rs} I_{2r}(s) - I_{2r}(0)] \\ &\approx 1 - \frac{2}{\sqrt{\pi}} \frac{\cos{r^2} - i \sin{r^2}}{2r} \left[\frac{1}{2r} - i\right] \\ &+ \frac{i \omega}{2 \sqrt{\pi}} e^{-ir^2} \left[e^{2rs+is^2} \left(2r^{-1} - 2isr^{-2} + (i + 2s^2) \, r^{-3} \right) - 2r^{-1} - i r^{-3}\right], \end{aligned}

which is the implementation I use in my complex function plotter. (I go up to fourth order in the app.)