MATLAB: Calculate derivative of erfcx(x) function

derivatediff()erfcx()MATLABSymbolic Math Toolboxsymbolic toolbox

I want to know how i can calculater the derivative of an erfcx(x) function with the symbolic math toolbox. The following is what i tried so far:

>> syms x real
>> y = erfcx(x)
Undefined function 'erfcx' for input arguments of type 'sym'.

I also tried doing this with an anonymous function:

>> testing = @(x) erfcx(x)
testing =
  function_handle with value:
    @(x)erfcx(x)

but then i don't know how to use the diff() function, as i don't know what variable i should give: diff(testing(x)) <- instead of this x.

Best Answer

Try

x = 0:0.1:5;
analytical_derivative = -2/sqrt(pi) + 2.*x.*erfcx(x);
numerical_derivative = (erfcx(x+1.0e-6)-erfcx(x))*1e6;
plot(x,analytical_derivative,x,numerical_derivative)

So my guess is that the derivative of erfcx is erfcx'(x) = -2/sqrt(pi) + 2*x*erfcx(x) (and you can easily deduce this result by differentiating exp(x^2)*erfcx(x) using the product rule).

Related Solutions

MATLAB: Derivative of a Function

If you have the Symbolic Toolbox:

>> g = @(x) exp(x^2-2)
g = 
    @(x)exp(x^2-2)
>> syms x
>> diff(g(x))
ans =
2*x*exp(x^2 - 2)

MATLAB: Inverse of erfcx (scaled complementary error function)

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

Best Answer

Related Solutions

MATLAB: Derivative of a Function

MATLAB: Inverse of erfcx (scaled complementary error function)

Related Question