How to solve an equation with term beyond realmax?

I run some equation which include exponent term multiplied by erfc term:

A=exp(B)*erfc(C)

For certain values the exponent term returns Inf and I can't solve the equation. Since the erfc term tend to zero I can't do log for the equation because the log(erfc) term will result in -Inf.

How can I overcome this problem and solve the equation?

Thanks,

Answer 1

[Roger Stafford]

The following is a somewhat improved version of my previous answer. It guarantees that you will not get a NaN for the product of exp(B)*erfc(C), (though of course you can get an infinity or a zero if the actual product is overflowed or underflowed.) For C > 25.54 a seven term asymptotic expansion is used in place of erfc(C) except that the exp(-C^2) part is omitted. Otherwise erfc(C) is computed. Then the needed corrective multiplicative factor, either exp(B) or exp(B-C^2), is broken up into n factors if necessary.

Let B and C be defined.

   if C <= 26.54  % Use erfc
     P = erfc(C);
     B2 = B;
   else           % Use asymptotic expansion
     P = 1;
     for k = 11:-2:1  % This generates seven terms
       P = 1-k/2/C^2*P;
     end
     P = P/sqrt(pi)/C;
     B2 = B-C^2;
   end
   n = max(ceil(B2/log(.99*realmax)),1);
   for k = 1:n  % If necessary make correction in n separate factors
     P = P*exp(B2/n);
   end

The value of P is the required product of exp(B)*erfc(C).

Answer 2

[Roger Stafford]

Your trouble would occur when both B and C are too large. In your function which evaluates the product exp(B)*erfc(C) you could make a special case for the circumstance of having both B and C too large and replace it with a continued fraction approximation for erfc(C):

   exp(B)*erfc(C) = C/sqrt(pi)*exp(B-C^2)*1/(C^2+1/2/(1+1/(C^2+3/2)))

Combining the B and C in exp(B-C^2) would presumably avoid your overflow/underflow difficulty. Here I have carried out the continued fraction to three quotients, but you would have to decide how many are really needed for accuracy in this circumstance. You can see this continued fraction defined to as many quotient levels as you need in the "Continued Fraction Expansion" section at

   http://en.wikipedia.org/wiki/Error_function

There is also another approximation for erfc(C) in the Asmptotic Expansion Section of this article in which you could also combine B and C in exp(B-C^2) in a similar manner to the above.

   exp(B)*erfc(C) = 1/sqrt(pi)*exp(B-C^2)*(1/C-1/2/C^3)+3/4/C^5-...)

[Roi]

Thanks, It's a very elegant solution. Just a small question to make sure:

For which range the continued fraction & asmptotic Expansions are valid for?

Is it for large C i.e. erfc(C)==>0 only? right?

(I couldn't find conclusive answer)

Thanks,

Roi

[Roger Stafford]

The use of the asymptotic expansion (it's not a continued fraction) is for the range when C is greater than 26.54 where the answer to erfc(C) would be zero or perilously close to zero. By omitting the exp(-C^2) factor in the expansion and later making an appropriate correction using the value B-C^2 it makes it possible for you to have a value of B which would ordinarily cause exp(B) to produce an infinity, and yet obtain a perfectly finite non-zero answer for what is equivalent to the product exp(B)*erfc(C). In other words it avoids the infinity-times-zero case which would produce a NaN for you.