Simple integration in Matlab is off

Question

This is an incredibly mundane integration, but for some reason I can't spot the error. I am just performing a simple 3d integration in spherical coordinates which should return the value of 1.0:

Nr = N;
Nt = round(pi*N);
Np = round(2*pi*N);

rs = linspace(R0, Rf, Nr);
ts = linspace(0, pi, Nt);
ps = linspace(0, 2*pi, Np);

dr = rs(2)-rs(1);
dt = ts(2)-ts(1);
dp = ps(2)-ps(1);

%fprintf('Total number of nodes: %i\n', Nr*Nt*Np);

C = 1/((4/3)*pi);

fint = 0.0;
for ir = 2:Nr
  r = rs(ir);
  r2dr = r*r*dr;
  for it = 1:Nt-1
    t = ts(it);
    sintdt = sin(t)*dt;
    for ip = 1:Np-1
      p = ps(ip);
      fint = fint + C*r2dr*sintdt*dp;
    end 
  end 
end

This is the same thing as $$ \int_0^{2\pi} \int_0^\pi \int_0^1 C r^2 \sin(\theta) dr\, d\theta\, d\phi$$

where I take $$ C = \frac{1}{(4/3)\pi}$$ since that is the reciprocal of the volume of a sphere.

However, I find that my returned integral value has a consistent error which is approximately (1/N)/(1.0 - numeric_integral) ~= 0.66. So I'm making a mistake somewhere but I can't spot where or why.

Answer 1

You're using a first order accurate integration technique (the rectangle rule) and your error is proportional to $1/N\propto\Delta r,\Delta \phi,\Delta \theta$. This is exactly the type of convergence behavior you should expect.

If you use a more accurate integration scheme you will see faster convergence. You can easily implement the trapezoid rule for this integral and I would expect to see second order convergence: $\text{Error}\propto(1/N)^2$.

Edit to add more details:

The rectangle rule you're using is identical to the following (check your answers with this, they should match up to rounding errors):

%create a 3d grid of values
[PHI,THETA,R] = ndgrid(ps,ts,rs);

%function to be integrated
F = C*R.^2.*sin(THETA);
fint = sum(sum(sum(F(1:end-1,1:end-1,2:end))))*dr*dt*dp;

while the trapezoid rule can be implemented as

%create a 3d grid of values
[PHI,THETA,R] = ndgrid(ps,ts,rs);

%function to be integrated
F = C*R.^2.*sin(THETA);

% INTEGRATE:
% r direction
I = sum( F(:,:,1:end-1)+F(:,:,2:end), 3 );
% theta direction
I = sum( I(:,1:end-1)+I(:,2:end), 2 );
% phi direction
I = sum( I(1:end-1)+I(2:end) );
fint = I*dr*dp*dt/8;

Computing the errors for different $N$ and plotting on a log-log plot reveals their order of accuracy. As expected, the trapezoid rule is $\mathcal{O}(1/N^2)$ while the rectangle rule is $\mathcal{O}(1/N)$.