Triangles appearing in numerical integration

Question

I have created a simple program in Matlab, which calculates the area of a polar function using rectangles (numerical integration, image 1). The formula used it's the one on the image 2.

$result = \frac{1}{2} * \sum_1^N \big(r(i)^2 * resolution\big)$

Then the program plots the area in function of the resolution, and here is where the magic happens.

The red line it's the exact value.

Why do these triangles appear?

The Matlab's code it's the following:

func = @(x) x*0 + 3;
area = (1/2) * integral(@(x) func(x).^2, 0, 2*pi);
int = 500000;
result = zeros(1, int);

for dif=1:1:int
    tmp = pi - pi*(dif-1)/int;
    result(dif) = 0;
    for i=0:tmp:2*pi
        result(dif) = result(dif) + (tmp*func(i)^2);
    end
    result(dif) = result(dif)/2;
end

plot(result);
hold on
x = 1:1:int;
plot(x, area + x*0);
hold off

EDIT:

Following what Ymir had said, if:

$result(tmp) = \frac{( \frac{2\pi}{tmp} + z)*tmp*3^2}{2}$

we can isolate z:

$z = \frac{2*r - 18*\pi}{9*tmp}$

and then map it: But I still don't know the why, can you explain better what z is?

EDIT 2:

I think I figured out what z is.

We can define the function's result as this:

$r = \frac{cycles*tmp*3^2}{2} $

Where cycles should be $\frac{2\pi}{tmp}$, but due to non-exact values it doesn't. Instead,

$cycles = 1 + fix\big(\frac{2\pi}{tmp}\big)$

(where fix(a) it's the truncated value of a)

Mapping that it results on the same function.

for dif=1:1:int
    tmp = pi - pi*(dif-1)/int;
    z(dif) = 1+fix(2*pi/tmp);
    result(dif) = (9*tmp*z(dif))/2;
end

And we can see that when cycles changes the function loses it's continuity (as Ymir said).

FINAL EDIT:

My final conclusion was that the problem was on the code itself. So, I have modified it.

area = 9*pi; % real integration value
int = 2000;
result = zeros(1, int);

% divide a circle in n parts
% for n = a:b:c
% equivalent to C code: for(n = a; n <= c; n+=b)
for n=1:1:int
    resolution = 2*pi/n;
    result(n) = (9*resolution*(n+1))/2;
end

% graph
plot(result);
hold on
x = 1:1:int;
plot(x, area + x*0);
hold off

And that's it, the expected graph :)

Answer 1

We''ll examine your code:

func = @(x) x*0 + 3;

The function you have does not depend x, so it's basically the constants 3.

result(dif) = 0;
for i=0:tmp:2*pi
    result(dif) = result(dif) + (tmp*func(i)^2);
end
result(dif) = result(dif)/2;

That code block is the same as:

   result(dif) = 0;
for i=0:tmp:2*pi
    result(dif) = result(dif) + (tmp*3^2);
end
result(dif) = result(dif)/2;

Pay attention that because of the code line:

    tmp = pi - pi*(dif-1)/int;

Every iteration $result(dif)$ value decreases. The decrease takes a linear manner. But pay attention that the code line: for i=0:tmp:2*pi And because tmp deceases, the number of iteration increases. Therefor, you add more terms to your equation, so with the help of some math we could find the equation you have actually calculated.

Before we do that, pay attention that result is a vector with only one row. Thus automatically, Every column is assigned with an x (horizontal axis) number from 1 to 5000, by increments of 1. Column's value is the y (vertical axis) values that you see.

$(\frac{2\pi}{tmp}+z)\cdot tmp*3^2$

But z depends on the iteration. For example:

1:9:2, will result in 5 iteration. And 1:10:2, also will result in 5 iteration. Only 1:11:2 will result in 6. The numbers you computed - tmp, are not so round, and will result in some drastic changes, as you add one more/one less term to the expression.You should find z that I talk about, plot it, and you'll see that's what causes the problem.

I will show you another exmple:

$1:9:0.2343$ will result in 35 iterations. but $9/0.2343=38$, so you'll habe to cut 3 iterations.

Now we can properly define z: z is the correction to the value $2\pi /tmp$, so we get how many iterations the loop keeps adding values to result(dif).

I guess that you doubled $fun(i)$ by tmp, because you guessed that there are $2\pi /tmp$ iterations, but that's my point: there are less or more iterations. it depends.

Furthermore, I believe that your method to calculate the integral is incorrect. Maybe I'm not familiar with this technique. But if you want to learn more about how to calculate integrals numerically I advise you to read D. Kincaid and W. Cheney, “Numerical Analysis, Mathematics of Scientific Computing. Particularly I guess you will find much interest in Newton-Cotes methods.

So to conclude: I showed you why the triangles appear, and I believe you should use different method to calculate the integral. Your question was about the triangles, so my answer is related to it. You are very welcome to open an discussion on integration methods and discus what were you supposed to do.

Hope it solves the problem. If you find something that is unclear, please let me know, and I will fix it.