Integral from a function that has a singularity

Hello I want to take an integral from the function that has a singularity(pole), by quadgk but this common don't give a right answer .for instance an integral of the function 1/(x-1) from(-2,2) anyone can guide me thanks

Answer 1

[Star Strider]

The best I can do is to rely on the Symbolic Math Toolbox, that came up with a piecewise result. Perhaps you can use this in your code:

syms x L U 
f(x) = 1/(x-1);
intf = int(f,x, L, U)

intf = 
piecewise(L <= 1 & 1 <= U, int(1/(x - 1), x, L, U), 1 < L | U < 1, log(U - 1) - log(L - 1))

The ‘L’ and ‘U’ variables are the lower and upper limits of integration, respectively.

[Kurt Hakan]

thanks for your answer. I would like to ask my question from other point of view. MATLAB Saied in help, the quadgk common can solve the integral of functions that have a singularity but here cant solve the problem.??!!

[Star Strider]

My pleasure.

In R2016b, the quadgk function will do the integration without error, but will throw Warnings. The quadgk function will also accept a 'Waypoints' argument, where you can ‘inform’ it about the singularities that exist within the limits of integration, and integrate around them if you wish.

For example:

f = @(x) 1./(x-1);
int1 = quadgk(f, -2, 2)
int2 = quadgk(f, -2, 2, 'Waypoints',[1, 1], 'MaxIntervalCount',1E+7)

int1 =
          12.9845608756444

int2 =
        -0.926665144495417

They both throw warnings. You can also divide the integral into two regions, from -2 to 1 and 1 to +2 and add them.

I am not certain there are any good ways to deal with such problems, although there are applied mathematicians here who I hope will add to this discussion.

[Kurt Hakan]

int2 = -0.926665144495417 is not correct answer and in we divide the(-2,2) to (-2,1)and (1,2) and take an integral the answer is

int3 =

   -1.1859
but the correct answer is (-log(3)=-1.09861).

???????

[Star Strider]

If you specify the limits of integration to go from -2 to +2 with 'Waypoints' around the singularity (doing complex contour integration), you get exactly that result!

int4 = quadgk(f, -2, 2, 'Waypoints',[1+1i, 1-1i], 'MaxIntervalCount',1E+7)

int4 =
         -1.09861228866811

[Kurt Hakan]

thanks

[Star Strider]

As always, my pleasure.

Wolfram Alpha evaluates it using the Cauchy principal value

Answer 2

[Paul]

According to quadgk, for an integrand with a singularity we should "split the interval and add the results of separate integrations with the singularities at the endpoints," and "write the integral as a sum of integrals over subintervals with the singular points as endpoints, compute them with quadgk, and add the results."

f = @(x) 1./(x-1);
I = quadgk(f,-2,1) + quadgk(f,1,2)

Warning: Minimum step size reached near x = 1; singularity possible.

Warning: Minimum step size reached near x = 1; singularity possible.

I = -1.1859

The warning messages are a concern. Also, quadgk states that: "it can integrate functions that behave at an endpoint c like log|x-c| or |x-c|^p for p >= -1/2." But in this case p = -1, so we should be concerned about the result.

We can also try integral the same way

I = integral(f,-2,1) + integral(f,1,2)

Warning: Minimum step size reached near x = 1. There may be a singularity, or the tolerances may be too tight for this problem.

Warning: Minimum step size reached near x = 1. There may be a singularity, or the tolerances may be too tight for this problem.

I = -1.1906

Same concerns with error mesages. Interestingly, integral doesn't make any statements about limitations on behavior of the integrand at the endpoints. I wonder if that's just an oversight in the doc insofar as quadgk states: "quadgk and integral use essentially the same integration method. You should generally use integral rather than quadgk."

Another numerical approach would be to integrate not quite to/from the singularity:

I = @(a) integral(f,-2,1-a) + integral(f,1+a,2);

Then we have for some small value of a

I(1e-5)

ans = -1.0986

which is the expected answer

-log(3)

ans = -1.0986

More generally, we'd have to play around with the value of a getting smaller and smaller until we thing we've indentified the limiting value of I as a->0.

However, for this problem, the split integration is insensitive to the value of a

[I(1e-7),I(1e-5),I(1e-3),I(1e-1),I(0.5)]

ans = 1×5

-1.0986 -1.0986 -1.0986 -1.0986 -1.0986

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We can see this result as follows

clearvars
syms x real
syms a positive
f(x) = 1/(x-1)

f(x) = 

An anti-derivative is

F(x) = int(f(x),x)

F(x) = 

Splitting the integral into two parts would look something like this as function of a

I(a) = int(f(x),x,-2,1-a,'Hold',true) + int(f(x),1+a,2,'Hold',true) %1

I(a) = 

Evaluate with the anti-derivative

I(a) = F(1-a) - F(-2) + F(2) - F(1+a)

I(a) = 

If we expand the result

I(a) = expand(I(a))

I(a) = 

we see that our split integral is actually independent of the value of a. Or we could just show directly

I(a) = int(f(x),x,-2,1-a) + int(f(x),1+a,2) % 2

I(a) = 

Hence if we take the limit of I(a) in (1) or (2) as a ->0 from the right (recall that a is positive) we get I(0) = -log(3), which is the Cauchy Principal Value of the original integral of interest

int(f(x),x,-2,2,'PrincipalValue',true)

ans =