MATLAB: Inverse of erfcx (scaled complementary error function)

erfcx()

could anyone help me to find inverse of erfcx (scaled complementary error function)?

Best Answer

As I recall, erfcx solves the problem that erfc goes to zero roughly like exp(x^2) for large x. For large x, erfcx is roughly
1/(x*sqrt(pi))
As we can see, this is true.
ezplot(@(x) 1./(x*sqrt(pi)),[1,10])
hold on
ezplot(@erfcx,[1,10])
And for negative x, erfcx blows up, and does so extremely fast. So you can only seriously be interested in inverting values of erfcx less than or equal to 1, where x will be non-negative.
So to invert y=erfcx(x) for a given value of y, a trivial solution is to use fzero as long as y lies in the interval [0.25,1].
For y less than say 0.25, I'd just construct a simple Newton iteration scheme, which will be even faster than fzero here.
y0 = erfcx(3)
y0 =
          0.17900115118139
>> x = 1/(y0*sqrt(pi))
x =
          3.15187684450162
>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99324120463763
>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99998665644266
>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
          2.99999999994798
>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
                         3
>> x = x - (erfcx(x) - y0)/(2*x*erfcx(x) - 2/sqrt(pi))
x =
     3
Quadratic convergence can be so nice. In fact, a quick test shows this even works quite well for larger values of y, so you need not even bother with fzero over part of the range.
In fact, if you are really lazy, don't even bother with a convergence test. Just run the Newton loop always for about a half dozen iterations for any value of y no larger then 1 and you should get full double precision accuracy.
Even better, all of this is easily vectorized.
function x = erfcxinv(y)
  % erfcx inverse, for y no larger than 1.
  x = 1./(y*sqrt(pi));
  for n = 1:6
    x = x - (erfcx(x) - y)./(2*x.*erfcx(x) - 2/sqrt(pi));
  end
end
That should work pretty cleanly. No checks to ensure that y is not greater an 1, in which case the 6 Newton iterations will probably not be sufficient, and one would want to have a tolerance test for convergence.
erfcxinv([1 .5 .25 .125])
ans =
     -1.59466173914173e-17         0.769079771061314          2.05154342412484           4.4052266527517
>> erfcx(ans)
ans =
     1  0.5  0.25  0.125
Related Question