Numerical evaluation of integral gives warning message

Dear all,

I am trying to compute the following integral in MATLAB:

with corresponding code is:

p1_minusSign = @(r) exp(-2*pi * integral(@(x) x/(1-x^4/r^4), r, Inf, 'ArrayValued', 1)) * exp(-pi*r^2)*2*pi*r;
p2_minusSign = integral(p1_minusSign, 0, Inf, 'ArrayValued', 1),

Unfortunately, MATLAB gives me a bunch of warning messages of the following form: "Warning: Minimum step size reached near x = 5833.72. There may be a singularity, or the tolerances may be too tight for this problem."

As far as I can tell, I don't see any potential singularity in my expression. In contrast, if we just change the sign in the demoninator of the inside integral, it computes perfectly:

with corresponding code:

p1_plusSign = @(r) exp(-2*pi * integral(@(x) x/(1+x^4/r^4), r, Inf, 'ArrayValued', 1)) * exp(-pi*r^2)*2*pi*r;
p2_plusSign = integral(p1_plusSign, 0, Inf, 'ArrayValued', 1),

Any ideas on what the error might be and how to correct it?

Many thanks in advance for your help.

[Ameer Hamza]

It seems that there is actually a singularity in the inner integrand for x=r. For example, if r=1 the integrand become

x/(1-x^4) and 'x' changes from r=1 to inf

and at x=1

1/(1-1) = 1/0

This singularity will create problems in the computation of integral. Not to mention that in outer integral r goes from 0 to inf. So when r=0 and x=0, you will have

0/(1-0/0)

[Torsten]

Obviously, the inner integral has a singularity at x=r.

[John D'Errico]

Um, you CANNOT factor out -1 from the term 1-x^4/r^4, to get 1+x^4/r^4. Mathematics does not work that way.

Is there a singularity? OF COURSE THERE IS!

When x=r, you have

1- r^4/r^r = 1-1 = 0

So you have a divide by zero. But when x=r=0, then you have something that looks like 0/0 at one point in the domain of your integrals.

[Gunnar]

Ah, yes of course, you are right. Any ideas on how I should proceed to compute this integral?

I tried adding a small perturbation, x=r+delta, but I still got NaN as output. Maybe I should try to doing piecewise integration?

[Gunnar]

@John D'Errico: Thanks for your answer. I was not trying to go against the laws of mathematics, you misunderstood me.

[John D'Errico]

Defying the laws of mathematics will see you get severely punished. I am left wondering what is the fine?

Seriously, you cannot solve it. The inner integral is not finite. Piecewise integration is not going to solve the problem.

[Torsten]

Shouldn't the limits of the integral in the exp(...) expression be x=0 to x=r instead of x=r to x=Inf ?

[Gunnar]

@Torsten: No, the limits are fine I think. Basically, with respect to a a point in 2D, we are integrating with respect to an outer region from r to Inf, and an inner region from 0 to r.

Answer 1

[John D'Errico]

Ignoring the fact that you cannot just factor out -1 from that term in the integral, and get what you think you got, we have the added problem that your inner integral is not even finite.

syms x r

K = x/(1 - x^4/r^4)
K =
-x/(x^4/r^4 - 1)

Now, pick some arbitrary value for r. Lets say 1.

int(subs(K,r,1),[1,inf])
ans =
-Inf

But if you change the denominator in the simple way that is NOT mathematically correct, as you wrote, then it is of course finite, as expected.

K2 = x/(1 + x^4/r^4)
K2 =
x/(x^4/r^4 + 1)

int(subs(K2,r,1),[1,inf])
ans =
pi/8

There is a fundamental difference between those two kernels. So, K2 is a mountain that we can indeed climb. (Sorry about that. Could not resist it.) It would probably be a HW assignment in first year calc, to show the former case is not finite?

Anyway, the answer is your problem lacks a finite solution.

[John D'Errico]

Looking still at the inner integral, since the entire thing is irrelevant, as long as the inner integral is not finite,

suppose we transform the problem

int(x/(1-x^4/r^4),[r,inf])

Suppose we make it simpler?

u = x/r

Then du=dx/r, and the transformed integral becomes

u*r/(1-u^4)*r*du

with limits of [1,inf].

So, we can factor out r^2, since r is constant wrt u. Your inner integral becomes

r^2*int(u/(1-u^4),[1,inf])

But that is still -inf.

syms u
int(u/(1-u^4),[1,inf])
ans =
-Inf

[Gunnar]

What about doing complex contour integration around the singular point for my inside integral similar to what was proposed here: Integral from a function that has a singularity

Would it work in my case too? Will have to refresh my Complex Analysis.

Another option is maybe to construct a piecewise function like I was saying earlier using the Symbolic Math Toolbox Integral from a function that has a singularity #2

Also, what I found a bit puzzling about my problem is that the singular point is an end point. From what I gather, MATLAB actually recommends putting singular points on the boundary of the domain: Singularity on interior of integration domain

[John D'Errico]

I might accept that your double integral could possibly exist, even though the inner integral does not. Part of me worries if that can be true. But I certainly cannot argue that it will never be possible to resolve this problem without some significant thought on the matter. All I can do at the moment is argue that you cannot solve the problem as you have tried to formulate it.

The approach of constructing a piecewise solution does not apply here. For example, some integration kernels can have a singularity in them, yet still have a finite integral. A beta function is a simple case. Given the correct set of arguments, the beta function kernel can be singular at the end point, yet have a finite integral. Thus we can write this

int(x^(-1/2),[0,1])
ans =
2

and get the expected finite result, even though a numerical integration tool might balk at the singularity of the function on that domain. The integral still exists there as a finite value.

In your problem however, as it is written, the inner integral has no finite result. The problem is not that the numerical integration tool cannot resolve the issue. It is that for any value of r, the inner integral apparently does not exist.

[Gunnar]

Ok, let's say the inner integral as it is written diverges. Let's therefore change the lower limit of the inside integral from x=r:Inf to x=r+1:Inf:

with associated code:

p1_minusSign = @(r) exp(-2*pi * integral(@(x) x/(1-x^4/r^4), r+1, Inf, 'ArrayValued', 1)) * exp(-pi*r^2)*2*pi*r;
p2_minusSign = integral(p1_minusSign, 0, Inf, 'ArrayValued', 1),

This corresponds to adding a small perturbation. This is just to play around with the expression a bit. Now the error message is different: Warning: Infinite or Not-a-Number value encountered.

Where is the error coming from now? As far as I can tell the singularity at x=r is no longer, yet we get NaN as output.

[Torsten]

Execute this code and you'll see the reason since you will have to take exp(-2*pi*F):

p1_minusSign = @(r)integral(@(x) x/(1-x^4/r^4), r+1, Inf, 'ArrayValued',1);
r=0:0.1:20;
for i=1:numel(r)
  F(i)=p1_minusSign(r(i));
end
plot(r,F);

[Gunnar]

Right, it grows unbounded/diverges. In the end it's due to the fact that the sign in the exponential function is positive.

[Torsten]

That's why I asked about the limits of the inner integral ...

[John D'Errico]

And my point is, you cannot compute the inner integral numerically. It is inf. Once it becomes inf, nothing you do with the result will be usable.

[Walter Roberson]

With the +1 the inner integral goes to

(1/4)*r^2*ln((2*r+1)/(2*r^2+2*r+1))

As r goes to infinity the ratio goes to 0, leading to ln(0) which is negative infinity. But you have exp(-2*pi*that) so you are getting exp(infinity)

Answer 2

[Walter Roberson]

Maple says that the inner integral is infinite for non-negative r.

When you substitute that in to the outer expression, you get

 int(exp(-Pi*r^2)*r*infinity, r, 0, infinity)

The exp(-Pi*r^2) term is nonnegative, and r is nonnegative, so we can see by inspection that the result will be infinity.