Numerical blowup problem in a fractional function in R

Ciao,

I am working with this function in R:

betaFun = function(x){
  if(x == 0){
    return(0.5)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

The function is smooth and well defined for every x (at least from a theoretical point of view) and in 0 the limit approach to 0.5 (you can convince yourself about this by using Hopital theorem).

I have the following problem:

i.e. the fact that, due to the limit, R wrongly compute the values and I get a blowup in 0.

Here I report the numerical issue:

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13)  
sapply(x, betaFun)

[1] 5.000083e-01 5.000442e-01 2.220446e+00 0.000000e+00 0.000000e+00 1.111111e+10

As you can see the evaluation is pretty weird, in particular last one. I thought that I could solve this problem by defining the missing value in 0 (as you can see from the code) but it is not true.

Do you know how can I solve this numerical blow up problem?

I need high precision for this function since I have to invert it around 0. I will do it using nleqslv function from nleqslv library. Of course the inversion will return wrong solutions if the function has numerical problems.

Your problem is that you take the quotient of two numbers with very small absolute values. Such numbers are only represented to floating point precision.

You don't specify why you need these function values for x values close to zero. One easy option would be coercion to high precision numbers:

library(Rmpfr)  
betaFun = function(x){
  x <- mpfr(as.character(x), precBits = 256) 
  #if x is calculated, you should switch to high precision numbers for its calculation
  #this step could be removed then

  #do calculation with high precision, 
  #then coerce to normal precision (assuming that is necessary)
  ifelse(x == 0, 0.5, as((1 + exp(x) * (x - 1)) / (x * (exp(x) - 1)), "numeric"))
}  

x = c(1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-13, 0) 
betaFun(x)
#[1] 0.5000083 0.5000001 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000

Numerical simulation of blowup in nonlocal reaction–diffusion , Part of the blowup theory for the time fractional ordinary differential Paper [18] analyzes this moving mesh method for the problems with linear C.J. Budd, R. Carretero-Gonz�lez, R.D. RussellPrecise computations of� $$\begin{align*} \left|\psi(\alpha,\beta) - \psi(1,1)\right| &= \left|(\alpha -1) \left[ -\frac1{\theta_3^2} (r^{1/\theta_3}\ln r) \sin\left( \frac{3\pi}{4\theta_3

I think that you are losing accuracy in the evaluation of exp(x)-1 for x close to 0. In C if I evaluate your function as

double  f2( double x)
{   return (x==0)   ? 0.5
            : (x*exp(x) - expm1(x))/( x*expm1(x));
}

The problem goes away. Here expm1 is a math library function that computes exp(x) - 1, without losing accuracy for small x. I'm afraid I don't know if R has this, but you'd hope it would.

I think, though, that you would be better to test for |x| was sufficiently small, rather than 0.0. The point is that for small enough x both x*exp(x) and expm1(x) will be, as doubles, x, so their difference will be 0. To keep maximum accuracy may need to add a linear term to the 0.5 you return. I've not worked out precisely what 'sufficiently small should be, but it's somewhere around 1e-16 I think.

Blow-up of error estimates in time-fractional initial-boundary value , The numerical analysis of differential equations containing fractional Let Ω be a bounded domain in Rd, where d∈{1,2,3}. A technical comment: in the problems and papers that we consider the function g(x,⋅) is� Blow-up and global solutions for a class of time fractional nonlinear reaction–diffusion equation with weakly spatial source Author links open overlay panel Jianxiong Cao a b Guojie Song c Jie Wang a Qihong Shi a Sujing Sun d

As you notice, you are encountering the problem near zero. The roots of both the numerator and denominator are zero. And as the OP mentioned, using L'Hôpitcal, you notice that in that f(x) = 1/2.

From a numerical point of view, things go slightly different. Floating points will always have an error as not every Real number can be represented as a floating point number. For example:

exp(1E-3)  -1 = 0.0010005001667083845973138522822409868               # numeric
exp(1/1000)-1 = 0.001000500166708341668055753993058311563076200580... # true
                                  ^

The problem in evaluating numerically exp(1E-3)-1 already starts at the beginning, i.e. 1E-3

1E-3 = x   = 0.0010000000000000000208166817117216851
exp(x)     = 1.0010005001667083845973138522822409868
exp(x) - 1 = 0.0010005001667083845973138522822409868
  1. 1E-3 cannot be represented as a floating point, and is accurate upto 17 digits.
  2. IEEE will give the closest floating point value possible to the true value of x, which already has an error due to (1). Still exp(x) is only accurate upto 17 digits.
  3. By subtracting 1, we get a bunch of zero's in the beginning, and now our result is only accurate upto 14 digits.

So now that we know that we cannot represent everything exactly as a floating point, you should realize that near zero, it becomes a bit awkward and both numerator and denominator become less and less accurate, especially near 1E-13.

numerator_numeric(1E-13) = 1.1102230246251565E-16
numerator_true(1E-13)    = 5.00000000000033333333333...E-27

Generally, what you do near such a point is use a Taylor expansion around zero, and the normal function everywhere else:

betaFun = function(x){
  if(-1E-1 < x && x < 1E-1){
    return(0.5 + x/12. - x^3/720. + x^5/30240.)
  }
  return( ( 1+exp(x)*(x-1) )/( x*(exp(x)-1) ) )
}

The above expansion is accurate upto 13 digits for x in the small region

A numerical algorithm for blow-up problems revisited, Since the numerical blow-up time is defined by an infinite sum, which implies that one needs to compute infinite times to achieve blow-up, this method cannot be carried out in real computation. analysis of semilinear parabolic problems with blow-up solutions. Rev. R. Acad. Sync. Arch. Rational Mech. 3 Numerical blowup problem in a fractional function in R Feb 8 '19 3 Finding distinct pairs {x, y} that satisfies the equation 1/x + 1/y = 1/n with x, y, and n being whole numbers Jul 22 '19 3 Triangular matrix get() in one dimensional array Sep 7 '17

Numerical approach to fractional blow-up equations , Numerical approach to fractional blow-up equations with Atangana-Baleanu derivative We also suggest the Fourier spectral method as an alternate approach to the fractional derivatives, Caputo and Fabrizio tackled this problem effectively by drug-responsive and drug-resistance cells are denoted by r N ; r T ; r R with� We investigate pointwise upper bounds for nonnegative solutions u(x,t) of the nonlinear initial value problem 0≤(∂ t −Δ) α u≤u λ in R ⁿ ×R,n≥1, u=0in R ⁿ ×(−∞,0) where λ

(PDF) Numerical integration of blow-up problems on the basis of , The third method is based on adding to the original equation of a differential constraint, Some examples of nonclassical blow-up problems are considered. culation of the right-hand side of (5) for fractional values of γbecomes impossible Transformation or Function Max. interval Stepsize Grid points. In literatures, the fractional derivatives are widely applied to some nonlinear problems, such as the fractional Schrödinger equation [10,11], the fractional Ginzburg-Landau equation [12][13][14

[PDF] Blow-up dynamics for the aggregation equation with , length scale, it is a challenging numerical problem to capture the mass function M(r, t) := ∫ r. 0 u(r)rd−1dr, which satisfies a local PDE when K is Newtonian; our numerical method to the fractional porous medium equation� In this paper, the fractional Broer&#x2013;Kaup (BK) system is investigated by studying its novel computational wave solutions. These solutions are constructed by applying two recent analytical schemes (modified Khater method and sech&#x2013;tanh function expansion method). The BK system simulates the bidirectional propagation of long waves in shallow water. Moreover, it is used to study the

Comments
  • I don't see that with curve((1+exp(x)*(x-1) )/( x*(exp(x)-1)), -10, 10, n = 1e5). Please provide a full reproducible example.
  • @Roland be aware that the points you select are only upto 1E-4. The OP computes 1E-13
  • Could you please make the plot compatible with the formula used? For the domain [-10,10] as depicted, the given function does not come that close to the asymptotic lines, the asymptotic lines are at y=0 and y=1, and the singularity is much tighter around x=0, where the value is 0.5, not 0.
  • Thx, is works fine except for x = 0. It is enough to add a if-statement to manually assign value at x = 0.
  • Cool, it keeps the vectorial structure!
  • Yes, there is a expm1 function in R.