Published May 5, 2021 ｜4 minute read
This post documents the implementation of the error function e r f ( z ) \mathrm{erf}(z) erf ( z ) complex function plotter . For small z z z Abramowitz & Stegun . For large z z z
Let ω : = e i π / 4 \omega := e^{i\pi/4} ω := e iπ /4 z = r ω + s i ω z = r\omega + s \, i \omega z = r ω + s iω 0 0 0 r ω r \omega r ω r ω r \omega r ω r ω + s i ω r \omega + s \, i \omega r ω + s iω
π 2 e r f ( z ) = ∫ 0 z e − x 2 d x = ∫ 0 r e − ( t ω ) 2 ω d t + ∫ 0 s e − ( r ω + t i ω ) 2 i ω d t = ω [ C ( r ) − i S ( r ) ] + i ω e − i r 2 ∫ 0 s e 2 r t + i t 2 d t , \begin{aligned}
\frac{\sqrt{\pi}}{2} \, \mathrm{erf}(z) &=
\int_0^z e^{-x^2} \, \mathrm{d}x \\
&= \int_0^r e^{-(t\omega)^2} \omega \, \mathrm{d}t
+ \int_0^s e^{-(r \omega + t i \omega)^2} i \omega \, \mathrm{d}t \\
&= \omega \, [C(r) - i \, S(r)] + i \omega e^{-ir^2} \int_0^s e^{2rt + it^2} \, \mathrm{d}t,
\end{aligned} 2 π erf ( z ) = ∫ 0 z e − x 2 d x = ∫ 0 r e − ( t ω ) 2 ω d t + ∫ 0 s e − ( r ω + t iω ) 2 iω d t = ω [ C ( r ) − i S ( r )] + iω e − i r 2 ∫ 0 s e 2 r t + i t 2 d t , where S ( x ) S(x) S ( x ) C ( x ) C(x) C ( x )
C ( x ) − i S ( x ) = π 2 ω − 1 − [ 1 + O ( x − 4 ) ] cos x 2 − i sin x 2 2 x [ 1 2 x 2 − i ] C(x) - i \, S(x) = \frac{\sqrt{\pi}}{2} \omega^{-1} - [1 + \mathcal{O}(x^{-4})] \frac{\cos{x^2} - i \sin{x^2}}{2x} \left[\frac{1}{2x^2} - i\right] C ( x ) − i S ( x ) = 2 π ω − 1 − [ 1 + O ( x − 4 )] 2 x cos x 2 − i sin x 2 [ 2 x 2 1 − i ] for large positive x x x a > 0 a > 0 a > 0
I a ( x ) : = e − a x ∫ − ∞ x e a t e i t 2 d t = a ∫ − ∞ x ∫ x ∞ e a ( t − s ) e i t 2 d s d t . \begin{aligned}
I_a(x) &:= e^{-ax} \int_{-\infty}^x e^{at} e^{it^2} \, \mathrm{d}t \\
&= a \int_{-\infty}^x \int_x^{\infty} e^{a \, (t-s)} e^{it^2} \, \mathrm{d}s \, \mathrm{d}t.
\end{aligned} I a ( x ) := e − a x ∫ − ∞ x e a t e i t 2 d t = a ∫ − ∞ x ∫ x ∞ e a ( t − s ) e i t 2 d s d t . Making the change of variables ( u , v ) : = ( a ( s − t ) , x − t ) (u, v) := (a \, (s-t), \, x-t) ( u , v ) := ( a ( s − t ) , x − t )
I a ( x ) = ∫ 0 ∞ e − u ∫ 0 u / a e i ( x − v ) 2 d v d u . I_a(x) = \int_0^\infty e^{-u} \int_0^{u/a} e^{i \, (x-v)^2} \, \mathrm{d}v \, \mathrm{d}u. I a ( x ) = ∫ 0 ∞ e − u ∫ 0 u / a e i ( x − v ) 2 d v d u . The inner integral can be evaluated as a sum of four Fresnel integrals, and is therefore bounded by a constant. Then for large a a a
∫ 0 u / a e i ( x − v ) 2 d v = e i x 2 ∫ 0 u / a e − 2 i x v + i v 2 d v = e i x 2 ∫ 0 u / a [ 1 − 2 i x v + i v 2 − 2 x 2 v 2 + O ( v 3 ) ] d v = e i x 2 [ u a − i x u 2 a 2 + ( i − 2 x 2 ) u 3 3 a 3 + O ( a − 4 ) ] . \begin{aligned}
\int_0^{u/a} e^{i \, (x-v)^2} \, \mathrm{d}v &= e^{ix^2} \int_0^{u/a} e^{-2ixv + iv^2} \, \mathrm{d}v \\
&= e^{ix^2} \int_0^{u/a} [1 - 2ixv + iv^2 - 2x^2 v^2 + \mathcal{O}(v^3)] \, \mathrm{d}v \\
&= e^{ix^2} \left[\frac{u}{a} - \frac{ixu^2}{a^2} + \frac{(i - 2x^2) \, u^3}{3a^3} + \mathcal{O}(a^{-4})\right].
\end{aligned} ∫ 0 u / a e i ( x − v ) 2 d v = e i x 2 ∫ 0 u / a e − 2 i xv + i v 2 d v = e i x 2 ∫ 0 u / a [ 1 − 2 i xv + i v 2 − 2 x 2 v 2 + O ( v 3 )] d v = e i x 2 [ a u − a 2 i x u 2 + 3 a 3 ( i − 2 x 2 ) u 3 + O ( a − 4 ) ] . Evaluating the standard gamma integrals, we thus obtain
I a ( x ) = e − i x 2 [ a − 1 − 2 i x a − 2 + ( 2 i − 4 x 2 ) a − 3 + O ( a − 4 ) ] . I_a(x) = e^{-ix^2} \, [a^{-1} - 2ixa^{-2} + (2i-4x^2) \, a^{-3} + \mathcal{O}(a^{-4})]. I a ( x ) = e − i x 2 [ a − 1 − 2 i x a − 2 + ( 2 i − 4 x 2 ) a − 3 + O ( a − 4 )] . We were a bit handwavy with the justification for Fubini. Indeed, the above expansion works well for x > 0 x>0 x > 0 x < 0 x < 0 x < 0 I a ( x ) I_a(x) I a ( x ) x x x
∫ 0 x e a t e i t 2 d t = e a x I a ( x ) − I a ( 0 ) ≈ e a x + i x 2 [ a − 1 − 2 i x a − 2 + ( 2 i − 4 x 2 ) a − 3 ] − a − 1 − 2 i a − 3 \begin{aligned}
\int_0^x e^{at} e^{it^2} \, \mathrm{d}t &= e^{ax} \, I_a(x) - I_a(0) \\
&\approx e^{ax + ix^2} \, [a^{-1} - 2ixa^{-2} + (2i - 4x^2) \, a^{-3}] - a^{-1} - 2ia^{-3}
\end{aligned} ∫ 0 x e a t e i t 2 d t = e a x I a ( x ) − I a ( 0 ) ≈ e a x + i x 2 [ a − 1 − 2 i x a − 2 + ( 2 i − 4 x 2 ) a − 3 ] − a − 1 − 2 i a − 3 works very well, and could probably be derived rigorously with more effort. Combining all the above, we have
e r f ( z ) = 2 π ω [ C ( r ) − i S ( r ) ] + 2 π i ω e − i r 2 [ e 2 r s I 2 r ( s ) − I 2 r ( 0 ) ] ≈ 1 − 2 π cos r 2 − i sin r 2 2 r [ 1 2 r − i ] + i ω 2 π e − i r 2 [ e 2 r s + i s 2 ( 2 r − 1 − 2 i s r − 2 + ( i + 2 s 2 ) r − 3 ) − 2 r − 1 − i r − 3 ] , \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} - ir^{-3}\right],
\end{aligned} erf ( z ) = π 2 ω [ C ( r ) − i S ( r )] + π 2 iω e − i r 2 [ e 2 r s I 2 r ( s ) − I 2 r ( 0 )] ≈ 1 − π 2 2 r cos r 2 − i sin r 2 [ 2 r 1 − i ] + 2 π iω e − i r 2 [ e 2 r s + i s 2 ( 2 r − 1 − 2 i s r − 2 + ( i + 2 s 2 ) r − 3 ) − 2 r − 1 − i r − 3 ] , which is the implementation I use in my complex function plotter . (I go up to fourth order in the app.)
