Trapezoidal Estimator for the integral

I tried here to find a trapezoidal estimate of an area between different intervals(Delta_x); however, errors found are not correct numbers according to others solutions(I checked). I couldnt find my mistake.

Could you look at it?

clc; clear;
f1=@(x) 1+sin(pi*x/100)-sin(pi*x/25);
Area1_actual=integral(f1,0,100);
%%Trapezoidal estimator
c=0;
for Delta_x=[10 5 2.5 1 1/2.5 1/5 1/10];
    x = 0:Delta_x:100;
        for j=1:(length(x)-1)
            Trap_area = Delta_x/2*(f1(x(j))+f1(x(j+1)));
            mat_traparea_Tr(j,1) = Trap_area;
        end
    c=c+1;
    Total_area_Tr(c,1) = sum(mat_traparea_Tr);
    Error_Tr(c,1) = Total_area_Tr(c,1)-Area1_actual;
end
Delta_x=[10 5 2.5 1 1/2.5 1/5 1/10];
plot(Delta_x,abs(Error_Tr),'b-'); title('Trapezoidal estimator')
xlabel('DeltaX'), ylabel('Errors'); set(gca, 'XDir','reverse')

Answer 1

[David Goodmanson]

Hi Furkan,

case 1: you are calculating Area2 with the trapezoidal approximation, but comparing it to Area1_actual when computing the error. Comparing it to Area2_actual gives much more satisfactory results.

case 2: Area2_actual is computed from 0 to 10, while your trapezoidal calculation is from 1 to 10.

It's true, the devil is in the details.

[FURKAN CEVAHIR]

Oh, just edited, you're right but this was just the mistake while copying. I think the problem is different thing.

[John D'Errico]

After editing was done to correct the error found, what do you think to be the problem?

Total_area_Tr
Total_area_Tr =
           163.13751514675
          163.531023680874
          163.629248948393
          163.656741162872
          163.661139476512
           163.66176779711
          163.661924876872
>> Area1_actual
Area1_actual =
          163.661977236758
>> Error_Tr
Error_Tr =
        -0.524462090007688
        -0.130953555884616
       -0.0327282883654334
      -0.00523607388655023
     -0.000837760245900654
     -0.000209439648017451
     -5.23598861263963e-05

So you have an integral that is actually not that poorly estimated, going from increments of 10 down to a fairly small number. Each estimate is a small under-estimate, not uncommon for a trapezoidal method. But are your estimates correct? That is trivial to test.

trapz(0:10:100,f1(0:10:100))
ans =
           163.13751514675

So the estimate for Delta_x(1) is perfectly correct. So, I'll repeat my question. Why do you think it is wrong? I see no such evidence.

My recommendation is you accept the answer given by David. As it is now, your code now seems to work correctly.

[FURKAN CEVAHIR]

For example, this is another script that aims to give the same estimates, but giving different results I dont understand why.

clc; clear;
f1=@(x) 1+sin(pi*x/100)-sin(pi*x/25);
Area1_actual=integral(f1,0,100);
Delta_x=[10 5 2.5 1 1/2.5 1/5 1/10];
Error_Tr_f1=zeros(length(Delta_x),1);
for i=1:length(Delta_x)
    dx = Delta_x(i);
    x=1:Delta_x(i):100;
    %%Trapezoidal estimator
        Trap_area_f1=zeros(length(Delta_x)-1,1);
        for j=1:length(x)-1
            Trap_area_f1(j) = dx/2*(f1(x(j))+f1(x(j+1)));
        end
    Error_Tr_f1(i) = sum(Trap_area_f1)-Area1_actual;
end

[David Goodmanson]

Hi Furkan,

see addendum to my original answer.

[FURKAN CEVAHIR]

Do you suggest any solution for that?

[John D'Errico]

It aims to do so. But is incorrect in what it does. I see at least two clear errors.

First, a minor one:

Trap_area_f1=zeros(length(Delta_x)-1,1);

If you would initialize a vector, it makes sense to initialize it to the correct length. Otherwise, preallocation is just a waste of time. Better would be:

Trap_area_f1=zeros(length(x)-1,1);

Of course, the error there is a minor one, since it does not actually create a bug. Just poorly written, inefficient code.

More important?

x=1:Delta_x(i):100;

I might wonder what are the limits of integration? Do you think that might be important? ;-) (Hint.)

You can write lots of scripts that give different answers. Most randomly written scripts would give the wrong result, of course. As I showed rather carefully however, the results I showed in my comment (as generated by the code you wrote after you fixed the error found by David) were actually correct.

[FURKAN CEVAHIR]

So, to make sure, you say that my script with this form in the question is completely true, right?