integration of (1/x-x)^alpha

Hi,

here is my code

function y = f(r, alpha)
F = @(x) (1./x-x).^alpha;
y = integral(F,0,r);
end

Then if i type

f(1/2,0.2)

i get

ans =

     0.707008459146274

but the problems start when alpha gets closer to 1, for example:

f(1/2, 0.8)
Warning: Infinite or Not-a-Number value
encountered. 
> In funfun\private\integralCalc>iterateScalarValued at 349
  In funfun\private\integralCalc>vadapt at 132
  In funfun\private\integralCalc at 75
  In integral at 88
  In f at 3 

ans =

     Inf

How can I compute this integral? Thank you for reading

Answer 1

[Roger Stafford]

Yes, you are right, Camille. If alpha < 1, your function is integrable, but the infinity in its integrand at x = 0 is apparently proving to be too difficult for matlab's 'integral' function to evaluate it properly.

However, if you make the following simple change of variable, you should encounter no such difficulty. To simplify notation use 'a' instead of 'alpha'. Define variable z as:

   z = x^(1-a)

Then we have

   dz = (1-a)/x^a*dx

and

   (1/x-x)^a*dx = (1-x^2)^a/x^a*dx = (1-z^(2/(1-a)))^a/(1-a)*dz

Hence you should evaluate the integral

   F = @(z) (1-z^(2/(1-a)))^a/(1-a);
   y = integral(F,0,r^(1-a));

Matlab should have no trouble with this integral provided a < 1, since it has no singularities in its integrand.

Answer 2

[Sean de Wolski]

1./1-1 yields a 0 and then 0 raised to any negative power is inf. What do you expect at 1?

[Camille]

Thank you for your response. The integral is well-defined for: 0 < r < 1 and alpha < 1 yet it won't compute when r=1/2 and alpha=0.8 (and in general when alpha gets too close to 1). So I'm looking for other means of computation.