Hi,

Analytical way for example is to represent exponent as Tailor series with an error that you can define according to your problem. Following integration step you will have an analytical solution to the problem with finite error.

Regards,

Jacob Nagler

With values like 1e-800 you must use arbitrary precision numbers. The standard double precision number cannot represent numbers smaller than about 1e-308.

Have you tried setting $maxExtraPrecision to a high number like 1000 or more to make sure your small values around 1e-800 are reperesentable by Mathematica?

Dear Anton,

Thanks for your very valuable clue. I forgot about $maxExtraPrecision. I have to setup it according to your proposition. I also found the similar reference:

https://reference.wolfram.com/language/ref/$MaxExtraPrecision.html

Dear Jacob,

I agree with your proposition because I thought about Tailor series. I have to calculate my example and compare with results obtained by using other methods .

The problem for generic g seems to be tuff. However, if you thought about Taylor series, then probably you may have more information on g (e.g. g is analytic or it can be approximated by a polynomial).

I will describe an ideal algorithm to solve your problem in the generic case, in the sense that I do not know whether all steps are feasible in practice. Nevertheless, it still seems interesting from a theoretical perspective.

I suggest to compute the integral of h(x)= g[x] Exp[n f[x]] using Gaussian quadrature $\sum\limits_{i = 1}{m} w_i h(x_i)$.

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

This method is exact whenever the integrand is a polynomial of degree less or equal to m. Therefore, it is suitable for your purpose if your function h(x) can be accurately aproximated by polynomials.The unique problem with this approach is that you will need to evaluate the  formula for high values of m to get high accuracy and I am unaware of any practical algorithm for computing the quadrature weights for such high values of m (some research would be required).

Nevertheless, once this problem is solved, then it can be rigorously proved (under few assumptions) that the computed  values produced by this method have good accuracy.

In the standard floating point model, it is easy to prove that, if

expr_1, expr_2, .... expr_k are expression which can be computed with high relative precision and all these expression have the same sign, then the sum

expr_1+ expr_2+ .... +expr_k can also be computed with high relative precision.

This is the case for the truncated Taylor series for the exponential function

1 + x + x^2/2 + ....x^k/k!, provided that x > 0.

Thus, if your f(x) < 0, you can compute exp[nf(x)] as 1/exp(-nf(x)) with truncated Taylor series and the computed value will be accurate (compared with the original Taylor truncated series). Then, by assuming that the Taylor series give  accurate values for

exp[nf(x)]  and that g(x) can be computed with high relative precision, then the values of the integrand h required by the quadrature formula can also be computed with high relative precision. Moreover, such values are ALL positive.The important characteristic of the Gaussian quadrature is that the quadrature weights are ALL positive. Thus, the porducts w_ih(x_i) are also ALL positive and the sum of this values can be also computed with high relative precision.

This is only an informal statement, but a rigorous proof can be formulated, provided that the we formulate the right assumptions.

Best regards,

André.

Dear Andre,

Thank you for your very interesting and valuable suggestions.  I will try to apply them to my integrals.  But it will be a problem when f[x] has singularity. According to your link (Wikipedia)

“The method is not, for example, suitable for functions with singularities”  

Radoslaw, this is true. But when the integrand have singularities, then you need to compute an improper integral and this seems to be quite hard in general. In fact, if your function have singularities, then even expanding it in Taylor series would make no sense.

Dear Andre,

I agree with you. But we can consider the following simple example. Let’s our range of x is [0, 1] and f[x] changes monotonically from 0 to –Infinity for f[1]. As a consequence our Exp[n f[x]] changing from 1 to 0. What can we do in a such situation. As an obstacle we should consider the fact that we know only approximation of function f[x] because it is calculated numerically in the previous stage of calculations (also difficult integrals). For calculating f[x] I try to use InterpolatingFunction f[x].

Dear Radoslaw,

I do not have familiarity with "Mathematica", but if you

are using a function labeled "InterpolatingFunction ", then

this interpolant is certainly (with overhelming probability)

continuos. Therefore, in this case there will be no problem in the

integration step.

However, if your original function f has singularities, then perhaps

would not be a good idea to approximate it with this function

"InterpolatingFunction". This seems to be a major concern.

Still, if you are using interpolation to approximate f, then the values of

f (or an approximation of f) must be known at some grid. In this case,

it seems better to try to interpolate the integrand "Exp[n f[x]]" directly,

because the exponential can "remove" the effects of the singularities of f.

A concrete example is just the one you already gave:

if f is monotonic and has only one singularity at 1, then the integrand Exp[n f[x]]

can be extended to a continuos function.

Dear Andre,

I thought about “exponential effect” like you. Maybe it is good way to solve integrals which are similar to my previous example. What do you think about Maple. Akira said that  this software could be more  friendly then Mathematica? I have a little experience with this application.

Dear Akira,

Thank you for your suggestion. I wrote to Andre about a  little experience with this application. Do you have any interesting links to Maple that could help me to solve my problem?

I have no experience with maple too.

Dear Radoslaw,

Maybe I can help you if you tell me some concrete information about the function f(x) and the interval (a,b), in order to construct some kind of Gaussian rules (with the maximal degree of precision).

Best regrads,   Gradimir

Dear Gradimir,

Thank you for willingness to help me.

My question is connected with modelling physical phenomena. It consists of a few steps.

In the previous steps I found numerical values for two functions g1[x] and g2[x] for x in [0, 1]. These functions have asymptotes for x=1. I do not know the analytical formula (I will approximate them in the future) but we know that they changes monotonically from 0 to Infinity (for x in [0, 1]).

I also found numerical values for f1[x] and f2[x] for x in [0, 1]. These functions also have asymptotes for x=1. They change monotonically  from 0 to –Infinity (for x in [0, 1]).  I also do not know the analytical formula.

I would like to find in this step of calculation the value of parameter alpha which satisfies the equation

Integration[x^3 g1[x*alpha] Exp[n f1[x]],{x,0,1/alpha}= alpha^(3/2) Integration[x^3 g2[x/Sqrt[alpha]] Exp[n f2[x]],{x,0,1}

As you see the function g1 and g2 are scaled by unknown parameter alpha. Also the upper limits of integration is function of alpha. The number n are equal 10, 100,…,10000. I can consider ther alpha>=1 and for my caluculation is in the range [1,50]

I found numerically (using InterpolationFunction[]) value alpha but from the other hand I know solution for n=Infinity (another theory). If I compare my results with boundary values then I see that for n=10 alpha is less then boundary (it is OK) but for n=100 or more (it depends on other parameters) is greater then boundary. The differences I can only explain by numerical errors in finding described integrals.

“Pozdravleniya ot nashego obshchego druga”  Jacek G.

Best wishes

Radoslaw

Dear Radoslaw,

In the attached file I have a question.

Best regards to you and Jacek,

Gradimir

Dear Gradimir,

I made one error in the notation of my equation. On the right side of equation should be alpha^(-3/2) (not alpha^(3/2). Rest of the equation is very good written. As for the calculating I am afraid about scaling functions f1[x], f2[x] in the equation and about asymptotes. In the future I think about Pade approximation parts of these functions.  

As a very little, warm up to approximation of physical problems (similar to discussed) can be treated my last paper “Approximation of the inverse Langevin function revisited”

(in my profile or http://link.springer.com/article/10.1007%2Fs00397-014-0802-2#).

Here there is also asymptote for x=1. It is a particular case where there are not external fields (=0).

Best regards

Radoslaw

Dear Qutaibeh,

Thank you for your tips. I think that can be useful in a lot of different cases.

Best regards

Radoslaw