How to solve an equation with term beyond realmax?
2 views (last 30 days)
Show older comments
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,
0 Comments
Accepted Answer
Roger Stafford
on 18 Aug 2014
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).
0 Comments
More Answers (1)
Roger Stafford
on 17 Aug 2014
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-...)
2 Comments
Roger Stafford
on 19 Aug 2014
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.
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!