f(x) =
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
采纳的回答
Star Strider
2017-1-29
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.
6 个评论
Kurt Hakan
2017-1-29
Star Strider
2017-1-29
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
2017-1-29
编辑：Kurt Hakan
2017-1-29
Star Strider
2017-1-29
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
2017-1-29
Star Strider
2017-1-29
As always, my pleasure.
Wolfram Alpha evaluates it using the Cauchy principal value
更多回答（1 个）
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)
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.
I = integral(f,-2,1) + integral(f,1,2)
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)
which is the expected answer
-log(3)
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)]
We can see this result as follows
clearvars
syms x real
syms a positive
f(x) = 1/(x-1)
An anti-derivative is
F(x) = int(f(x),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
Evaluate with the anti-derivative
I(a) = F(1-a) - F(-2) + F(2) - F(1+a)
If we expand the result
I(a) = expand(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
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)
类别
在 帮助中心 和 File Exchange 中查找有关 Creating, Deleting, and Querying Graphics Objects 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
