# Simple integration in Matlab is off

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.

• Can you state in a formula what you are trying to integrate, just to make the code easier to read? – Wolfgang Bangerth Apr 14 '15 at 16:15

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$.

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)$.

• So it looks right? I was thinking since the error is so predictable and consistent that it was an error in the integral. I was expecting to the have the error be a bit more random and even have times where my final integral is less than the analytic value. – drjrm3 Apr 14 '15 at 18:08
• I think it looks fine. It's not a great way to compute 3d integrals in general, but it should work for this problem. For the rectangle rule the sign of the error will just depend on the slope of the function you're integrating. $r^2$ is always increasing on $[0,1]$ but $sin(\theta)$ goes up then down on $[0,\pi]$. Apparently because of the way you've chosen your rectangles and grid points the net effect in this case is a consistent sign on the error in your result. – Doug Lipinski Apr 14 '15 at 18:28