All the information I found about inverse error function $\text{erf}^{-1}(z)$ was about $z\in\mathbb{R}$. Also I found some Taylor expansions for it, but as the function is unbounded near $z=\pm1$, these expansions only converge in the disk $|z|<1$.
I want to look at real and imaginary parts of this function for $z\in\mathbb{C}$. I tried using Mathematica's InverseErf[z], but it appeared to only support real arguments. I then tried using FindRoot to determine values and plot, but I got somewhat strange results, which don't disappear when I increase WorkingPrecision (tried up to 50 decimal places).
Here's what I got for real (left) and imaginary (right) parts on DensityPlot from FindRoot:
Very light and very dark regions correspond to values outside of $-3<f<3$ range.
As one can see, using Taylor series won't help me with these strange regions, because all they are outside of disk of convergence.
So, my questions are:
- (answered) How can one compute $\text{erf}^{-1}(z)$ for $z\in\mathbb{C}$ including those $|z|\ge1$?
- Are there any packages which are able to compute it without me having to implement the algorithm?
- Are properties of this function for complex arguments described anywhere?
Let's denote $\text{erf}^{-1}(z)=\text{ierf}(z)$ to remove potentially misleading $^{-1}$ "power".
As $\text{ierf}(z)$ is an inverse of an integral of an elementary function, we can easily find its derivative:
$$\frac{d}{dz}\text{ierf}(z)=\left.\frac1{\frac{d}{dq}\text{erf}(q)}\right\rvert_{q=\text{ierf}(z)}=\frac{\sqrt{\pi}}2e^{\text{ierf}\,^2(z)}$$
We can see that the derivative has only the value of the function at the point of differentiation as a parameter. This means that all higher-order derivatives will also depend only on value of this function in a single point. Thus, we could expand it to series in some far enough from singularities point to get the function outside of $|z|<1$. But we'll need to find the value of this function at that point. Of course, we could find it by some root finding algorithm, but now let's go another way.
Since we know the derivative of our function, and it only depends on the value of the function at the point of differentiation, we can make a differential equation, which the function must satisfy:
$$\frac{df(z)}{dz}=\frac{\sqrt{\pi}}2e^{f^2(z)}$$
Now, we know that since $\text{erf}(0)=0$, then its inverse satisfies $\text{ierf}(0)=0$. This allows us to write a Cauchy problem for the first order ODE:
$$\cases{\frac{df(z)}{dz}=\frac{\sqrt{\pi}}2e^{f^2(z)} \\ f(0)=0}$$
This equation will give us all the values of $f$ for $-1<x<1$. Of course, this is not what we need, so let's instead define a function $g(t)=f(v(t))$, where $v(t)$ is a function, which parametrizes some contour in the complex plane. Our differential equation will then be:
$$\frac{dg(t)}{dt}=v'(t)\frac{\sqrt{\pi}}2e^{g^2(t)}$$
Let's see what happens if this contour is the line $v(t)=(1+i)t$. On the following plot blue and magenta are real and imaginary parts of $g(t)$, and black and green are real and imaginary parts of corresponding partial sum (up to $n=150$) of series expansion at $z_0=0$, for comparison:
We can see that Cauchy problem approach is more fruitful in the sense of giving the solution without limiting our input values range.
Let's now select a circle, which starts at $v(0)=0$ passes through $v(\pi)=2$, and ends back at $v(2\pi)=0$, thus avoiding (and going round) the singularity at $z=1$, i.e. $v(t)=1-e^{-it}$. On the next image colors denote the same.
At the left side of the plot, where series solution has not yet started to diverge, both solutions agree. But at the endpoint of the circle the solutions disagree. This is because we've travelled around a branch point. If go in another direction (i.e. selecting $v(t)=1-e^{it}$), we'll get yet another value at endpoint (this time with negative real and imaginary parts).
Here's how this looks on parametric plot along with partial sums of expansions one at $z=0$ and another at $z=6+0.1i$ (left image is real part, right is imaginary part, red and green curves are solutions of Cauchy problem, rendered according to values of $v$):
Now that we can evaluate $\text{ierf}(z)$ at some points, let's generalize this procedure to access any point in complex plane. For this, let's define $v(t)=a \tanh(t)e^{it}$. Denote the point we want to find $\text{ierf}$ of $p$. Our endpoint will be $T=\arg(p)$ ($\arg$ is defined here to have range $(-\pi,\pi]$), but we have to adjust it to select correct path, choosing a branch, for which travel from $0$ to our $p$ is nearest; if $p\ge0$, we select path such that $\Im [v]\ge0$, otherwise we make $\Im [v]<0$. So, final $T$ will be:
$$T=\begin{cases}\arg(p)-\pi&0\le \arg(p)<\frac\pi2\\ \arg(p)+\pi&-\frac\pi2\le\arg(p)<0\\ \arg(p)&\text{otherwise}\end{cases}$$
Next, we need to define $a$. It'll just be such that $v(T)=p$. Here's how our path will look example $p=1.5+0.9i$:
Now, if we make some executable code from this and try to render $\text{ierf}(z)$, here's how the function will look (left real, right imaginary parts):
These plots have some stains, I think it's because of some numerical difficulties Mathematica had while solving Cauchy problem for some points. Nevertheless, the final goal to be able to compute and render $\text{ierf}(z),\;z\in\mathbb{C}$ is met.