MATLAB: Exp(infty1) * erfc(infty2) = NaN

erfcinfinity

My problem is that I have the following analytical solution to a flow problem

The meaning of the parameters are as such not relevant to the problem, but the problem arises from the last term of A_i. The solution works perfectly for specific sets of parameters, but for the ones I need it for there is a problem. The exponential goes towards infinity, ~exp(2*10^4), where the complimentary error function goes towards zero, erfc(x) where x is in the interval 140 to infinity (however can MAYBE be limited to the interval of 140-200). How do I handle this? Matlab's limit of doubles seems to be 10^(+-308), which is not nearly enough to solve the problem, where the exp term is approx equal to 10^8685

Best Answer

You have several choices.

1. Use the symbolic Toolbox, which can handle high precision computations.

x = sym(200);
vpa(erfc(x))
ans =
4.6893592959836014980356997769923e-17375

2. Use my HPF toolbox, which can also handle high precision computations. I thought I implemented erfc, but its been a while. It should be doable though. I do know there is what looks like a nicely convergent expansion for large x for erfc. So if it turned out I did not implement erfc for large arguments, it would be pretty easy to write (I think.)

3. Use approximations, numerical analysis, etc. For example, I'd not be surprised if there is not a method to compute the log(erfc(x)), for large x. No matter what, you will still probably need a tool that can provide wide ranging computations on lots of digits, because you are adding terms there. Will you need thousands of significant digits? Or can some terms simply drop out when they are clearly insignificant? Again this is a numerical analysis question as much as anything.

Let me expand on my last comment. If the last term is insignifcant compared to the others, AND you are willing to work in double precision, you would merely compute the log of that last term. A poor approximation for the log of that expression would suffice there to know IF that term is indeed insignificant.

In the end, you will need to be careful and use good numerical methods no matter which approach you take here. Remember that the issue is not only the dynamic range of a double, but the number of significant digits it carries in the mantissa, and how many correct digits you really need in your result.

Related Solutions

MATLAB: How to solve an equation with term beyond realmax

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).

MATLAB: How to solve this summation with exponential in Matlab

If you are not careful, you will need to use a tool with literally thusands of digits of precision, as those terms will vary by thousands of digits in magnitude. Worse, if you even try to compute any of those terms, you will find that double precision either underflows or overflows.

vpa(factorial(sym(1247)))
ans =
8.458643820551953682887859707083e+3320

As such, in order to compute that series, you will need to use symbolic computations, or use a tool like my HPF. I'll do it using HPF, although I could have as easily used the symbolc toolbox too.

DefaultNumberOfDigits 100
S = hpf(0);
E = exp(hpf(-1211.5));
T = hpf(1);
for n = 0:1247
  S = S + E*T;
  T = T*1211.5/(n+1);
end
S
S =
0.8494300028975462166085239594512340095256328874564113849900813900015574782646329472378321829350836040

That is the lazy way to solve the problem. It involves no real effort, no thought about how to compute the result.

Can you solve the problem using double precision? Well yes. But again, you MUST be careful. And you must understand how to work with double precision in a way that does not underflow or overflow before you have gotten the result you need. The trick is often to use logs on these problems.

S = 0;
T = -1211.5;
for n = 0:1247
    S = S + exp(T);
    T = T + log(1211.5) - log(n+1);
end
S
S =
         0.849430002892729

In this sum, you should realize that almost every single term is effectively zero. So by using natural logs, I computed the log of every term, then exponentiate only after I have computed the log of the complete term.

As you should see, the two computations agree reasonably well. The double precision term is wrong in the last few digits, but you also need ot consider that every term in that series varies by as much as roughly 3 decimal digits. So that is probably as good as I can easily do using doubles.

Best Answer

Related Solutions

MATLAB: How to solve an equation with term beyond realmax

MATLAB: How to solve this summation with exponential in Matlab

Related Question