1

Dealing with numerical calculations having exponential terms often becomes painful, thanks to overflow errors. For example, suppose you have a probability density, P(x)=C*exp(f(x)/k), where k is very small number, say of the order of 10^(-5).

To find the value of C one has to integrate P(x). Here comes the overflow error. I know it also depends on the form of f(x). But for this moment let us assume, f(x)=sin(x).

How to deal with such problems?

What are the tricks we may use to avoid them?

Is the severity of such problems language dependent? If yes, in which language should one write his code?

ddas
  • 199
  • 1
  • 15
  • 3
    Isn't the usual trick calculating with logarithms of the actual values and using `+` instead of `*`? – tobias_k Jul 27 '18 at 11:27
  • 1
    This sounds like a question much better suited to https://scicomp.stackexchange.com Especially asking which programming language to use is clearly off-topic here. It is quite broad as well, ask here about concrete and focused problems to solve, not for possible tricks. – Vladimir F Героям слава Jul 27 '18 at 11:45
  • Regarding the last question, raw python has no integer overflow (i.e. as long as you don't use libraries dependent on C, etc.). E.g. `2**500 returns 3273390607896141870013189696827599152216642046043064789483291368096133796404674554883270092325904157150886684127560071009217256545885393053328527589376L` – isaacbernat Jul 27 '18 at 11:50
  • 4
    Step 1: Go analytically as far as possible! In your case, the following integral is very handy `Integral[exp(i*b*sin(x)) cos^(2n)(x),{x,-Pi/2,Pi/2}] = Sqrt(Pi) (2/b)^n Gamma(n+1/2) J_n(b)` with `b` and `n` a complex numbers and `J_n(x)` the Bessel function of the first kind. (Ref Eq 3.915.3 in Gradshteyn and Ryzhik) – kvantour Jul 27 '18 at 12:21
  • 1
    For what kind of problems do you need such huge evaluations ??? This looks pathological. –  Jul 27 '18 at 19:15
  • 2
    Yet another approach... If we consider P(x) = D * exp( a * (sin(x) - 1) ) (where D is a normalization constant and a = 1/k), then the quadrature of exp( a * (sin(x) - 1)) may be easier. For general f(x), maybe consider P(x) = D * exp( a * (f(x) - c) ) with c = max f(x) etc... – roygvib Jul 27 '18 at 23:22
  • @YvesDaoust In non-equilibrium statistical mechanics, solution of the Fokker-Planck equation often needs such huge evaluations. – ddas Jul 28 '18 at 02:38
  • @tobias_k you're right. Logarithms may help. – ddas Jul 28 '18 at 02:41
  • 1
    @ddas A literature search on existing approaches might be helpful. This is outside my area of expertise, but a minimal search finds the following, for example: Brian T. Park and Vahe Petrosian. "Fokker-Planck equations of stochastic acceleration: A study of numerical methods." *The Astrophysical Journal Supplement Series* 103 (1996): 255. – njuffa Jul 28 '18 at 17:37
  • My [answer](https://stackoverflow.com/questions/48990778/newton-raphson-iteration-unable-to-iterate/49019017#49019017) to a similar question might be helpful. – Scientist Jul 30 '18 at 17:21

2 Answers2

1

As I mentioned in the comments, I strongly advise using analytical methods as far as you can. However, If you want to compute integrals of the form

I=Integral[Exp[f(x)],{x,a,b}]

Where f(x) could potentially overflow the exponent, then you might want to renormalize the system a bit in the following way:

Assume that f(c) is the maximum of f(x) in the domain [a,b], then you can write:

I=Exp[f(c)] Integral[Exp[f(x)-f(c)],{x,a,b}]

It's an ugly trick, but at least your exponents will be small in the integral.

note: just realized this is the comment of roygvib

kvantour
  • 25,269
  • 4
  • 47
  • 72
  • @roygvib Feel free to write this as an answer, and I'll delete this one! – kvantour Jul 30 '18 at 10:36
  • Hi, no problem. I tried a bit of numerical experiments, and it worked with "not too large" f(x). But when it comes to 10^5 etc, the integrand becomes like a delta function, which may be problematic for numerical quadrature (because of the too narrow width)... But anyway, some (pre)scaling of the problem may be helpful. (Otherwise, I guess a more neat approach (specific to the problem) may be necessary like the paper shown in the comments.) – roygvib Jul 30 '18 at 12:10
0

One of the options is to use GSL - GNU Scientific Library (python and fortran wrappers are available). There is a function gsl_sf_exp_e10_e that according to the documentation

computes the exponential \exp(x) using the gsl_sf_result_e10 type to return a result with extended range. This function may be useful if the value of \exp(x) would overflow the numeric range of double.

However, I would like to note, that it's slow due to additional checks during evaluation.

P.S. As it was said earlier, it's better to use analytical solutions where possible.

Bracula
  • 335
  • 1
  • 3
  • 14