Numerical evaluation of multiple integral not working

Hi guys,

I'm trying to numerically calculate the integral after Equation (1) from here: Distances In Bounded Regions

The final purpose is to calculate the average distance between all points in a region.

The result of the integral should be 0.5214.

I'm calculating it by calling integral2 function like this:

clear all;
f = @(theta,r)(1-r.*sin(theta)).*(1-r.*cos(theta)).*r.^2;
ymax = @(theta)1./cos(theta);
xmax = pi/8;
8*integral2(f, 0, xmax, 0, ymax)

However I'm getting the wrong answer. What's the problem with my code?

Answer 1

[John D'Errico]

Really? Lets do it symbolically.

syms r theta
kernel = (1-r.*sin(theta)).*(1-r.*cos(theta)).*r.^2;

rint = int(kernel,r,[0,1/cos(theta)])
rint =
1/(12*cos(theta)^3) - sin(theta)/(20*cos(theta)^4)

sbar = 8*int(rint,theta,[0,pi/8])
sbar =
log(tan((5*pi)/16))/3 - 16/(15*(2^(1/2) + 2)^(3/2)) + (2*2^(1/2))/(3*(2^(1/2) + 2)^(3/2)) + 2/15

vpa(sbar)
ans =
0.24810023714331117778996075103314

So the symbolic toolbox questions the result in that paper.

pretty(sbar)
   /    / 5 pi \ \
log| tan| ---- | |
   \    \  16  / /            16                2 sqrt(2)        2
------------------ - ------------------- + ------------------ + --
         3                           3/2                  3/2   15
                     15 (sqrt(2) + 2)      3 (sqrt(2) + 2)

What did integral2 produce as a result?

f = @(theta,r)(1-r.*sin(theta)).*(1-r.*cos(theta)).*r.^2;
ymax = @(theta)1./cos(theta);
xmax = pi/8;
8*integral2(f, 0, xmax, 0, ymax)
ans =
       0.248100237144312

Which happens to be quite the same as that produced by the symbolic solution.

So I don't know where that article you found on the internet came from. But it seems to have at least one flawed claim in it. I won't spend the time to check their previous results, or the derivation that went into it. Anyone can post anything online, and with no validation of their results, you get what you pay for.

[Daniel]

John,

What a fantastic answer. Thank you so much for your help.

For the sake of doubts I'll try to distill what's wrong with the original equations.

Again, thank you. I thought I was losing my mind with this one.

[Daniel]

To provide a follow up:

The paper is indeed wrong. With xmax=pi/4 the result in the paper is achieved.

Thanks again John!

[John D'Errico]

I read over that formula about 5 times to make sure that it was written as it was in the paper. Its easy to not trust your own code. After all, how can something written in a paper be wrong? As I said, the internet needs proofreaders. Sadly, the pay is poor.