Numerical integration improves when step size is smaller, but in which direction?

Question

I am implementing a double integration of a function; numerically in Matlab. What I basically do is:

number_of_elements = 10;
x = linspace(0,3,number_of_elements);
y = linspace(0,3,number_of_elements);
dx = (3-0)/number_of_elements;
dy = (3-0)/number_of_elements;

for i = 1:number_of_elements
    for j = 1:number_of_elements
        f_j(j) = (x(i)^2 + y(j)^2) * dx * dy;
    end
    f_i(i) =  trapz(f_j) % trapezoidal rule
end

f = trapz(f_i) % trapezoidal rule

To do its tests, I am comparing my code's result to an analytically calculated result of simple functions.

For e.g: I am integrating the function: function

Calculating analytically using pen and paper, its result is 54. If my implementation is right, than I would expect my calculated value to be near 54 and would get nearer and nearer to 54 as I make my step sizes of integration smaller , i.e. the number_of_elements larger.

Good news that this is exactly what happens: results. On x axis is the number of elements(which increases with smaller step size). On y axis in the function value in blue, and real value in red.

Its all good. But then I tried an another function: function2. Here dr and dgama is similar to dx and dy as in 'function1'

The results are: results (function2). The limits for the first integral are R_hub = 0.3 to R_tip =1.25. The calculated value is in yellow and real value in black.

Its right that the results for both cases improves as the steps get smaller.

But why does for larger step sizes, the numerically calculated values starts from a higher value than the true value in the first case and then from a lower value than the true value in the second case.

Thanks in advance.

Answer 1

If your function $f$ is convex, such as $f(x,y) =x^2+y^2$, the trapezoid rule will overestimate the true integral, since the height of points on the interior of each trapezoid will be greater than the value of $f$ at the corresponding point. Similarly if $-f$ is convex the numerical integral will underestimate the true value.

For more complicated functions that contain a mix of convex, concave, and saddle regions, such as your second function, it's hard to say anything meaningful generically about the sign of the error. For a concrete $f$ you can compute it in terms of $f$'s higher-order derivatives by Taylor-expanding $f$, as sketched in the error analysis section on Wikipedia.