exp(infty1) * erfc(infty2) = NaN

Hi

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

Answer 1

[John D'Errico]

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.

[Walter Roberson]

log(erfc(x)) for large x approaches -infinity. I can't seem to get further than roughly 4.5E9 before the loss of precision overtakes any sensible result.

[John D'Errico]

Of course it will. I think you missed my point. First of all, we were told that x would be on the order of 100-200.

Regardless, the point is you can use a good approximation. And since you are taking the product of terms, one of which goes to +infinity, they will cancel each other out, to some extent. Then you can quickly learn whether the term is even worth the effort to compute accurately.

For example, take only the first term of the asymptotic expansion for erfc. That reduces to:

exp(-x^2)/(x*sqrt(pi))

So, for x==200, we will see:

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

vpa(exp(-x.^2)/x/sqrt(pi))
ans =
4.6894179115094684782688697934864e-17375

As you can see, that first term is really quite accurate.

Now, how about the product exp(y)*erfc(x), (with y=2e4)? We can compute the approximate log of that product very simply, as

y + -x^2 - log(x) - log(sqrt(pi))

And if you wanted to add in a few extra terms in the expansion, it would still be easy as pie (or is that easy as pi?) to gain 16 correct digits for large x.

y = 2e4;
x = 200;
y + -x^2 - log(x) - log(sqrt(pi))
ans =
       -20006

Of course, the exponential of that is incredibly small. But it would allow you to decide if there was even any reason to include that in a sum. Again, I have no idea how many digits are required in this, and I agree the odds are small that this would ever be significant. But you can now trivially compare the magnitude of that term with the other terms.

Suppose x were just a wee bit smaller?

y = 2e4;
x = 141;
y + -x^2 - log(x) - log(sqrt(pi))
ans =
       113.48

Now that term is quite large.

exp(ans)
ans =
   1.9198e+49

Answer 2

[Bastian Andersen]

Thanks for the extensive and fast answers! The last term is significant for certain ranges of the time t, so it will be required (e.g. for exp(20000)*erfc(140). The first thing I tried was also just to delete the term, but it provided bad results :) For other time steps, where x in erfc(x) becomes very large, it would however be possible to neglect, which may simplify the expression. That's the reason I think you could limit x to be in the range of x = 140..200, as t is the only changing variable (atm at least).

It is used for chromatography, where peaks of values different from zero arises at specific time intervals. I think they are this problematic to calculate as they approximate being shocks. It should be fine for terms to drop out when they are clearly insignificant. I have to iterate over a_i to fit the equation to something similar to this, where all other values are close to zero.

Do you recommend one of your mentioned methods above the others? Using the symbolic toolbox will probably increase the computational time, but if I specify a limit on whether the terms should be calculated or not I think it is fine. Regarding the approximation with the logarithm, I don't see how you handle the other terms when singling out the third and last? They can't just be separated. Using the asymptotic expansion seems reasonable, as it will always be in the asymptotic range. Something similar should be possible for the exponential term.

[John D'Errico]

Why did you feel it necessary to use an answer instead of a comment?

I don't see the problem here. You can easily compute the magnitude of each term, using an approximation and/or logs as necessary if there is an underflow or overflow.

That is, the first term is just an erfc. Again, you can compute the argument. If it is large, then again, we have a simple expansion if erfc will crap out.

The second term is the product of some stuff with an exponential. Again, it is trivial to either compute or see the relative magnitude.

The third term we have discussed already, but as I've shown, it is easy to compute the log. And it is easy to compute the logs of all three terms.

A simple trick would be to subtract off the maximum of the logs of all three terms, essentially scaling the problem by an additive shift of the logs. After you compute the sum of the exponentials of the three now shifted terms, then you have a known scale factor out front.

If the result will be so small that the sum will be an underflow or an overflow, then you can decide what to do. I don't know what you would do if the sum was going to underflow double precision arithmetic anyway, or overflow it for that matter. But that is for you to decide.

Finally, you ask which method is best. Again, it is trivial to compute the logs of each term, and do so in double precision. You can do virtually everything using doubles, until the very end when you can decide if the result will underflow or overflow a double. It does not sound as if you really need high precision arithmetic for this. But either sym or hpf would be adequate, although relatively slow. What you do with the result, and the precision required is your choice.

[Bastian Andersen]

Because it's 1 am here, so I just click any bottoms nearby (: Thanks for the help, seems like I got what I needed. Will look it through tomorrow

[John D'Errico]

:)