Trapezium Rule value not correct

Hello, I posted on here earlier about a problem with my function but I have managed to fix that problem. However the value for my trapezium rule is wrong and I was wondering what I had done wrong? Also my graph won't create the stairs and steps that I have created, the area underneath the graph is just filled in red instead of creating rectangles unless there is a better way to plot the rectangle rule? I have attached my code to this question, help would be fantastic as the code is due in for hand in on friday and I am very stuck!!!

[Geoff Hayes]

Beth - Were you comparing your output with that from the MATLAB trapz function? Note that if you just input the y vector from your code and run the command as trapz(y) then the result will not match with yours since trapz assumes that unit spacing was used to determine the values of y. Instead, try trapz(x,y) and the answer will be similar to yours.

[Geoff Hayes]

What code is there (in ni2.m) to plot the steps? I just see code to plot the function.

[beth]

I have to create an equation for the trapezium rule rather than use the trapz function. Ah I seem to have deleted the steps and stairs by accident, I had originally put stairs(x,y,'r'); steps(x,y,'r');

I typed this in after the line plot(x,y).

Answer 1

[the cyclist]

For your default program, the exact value is 2500 and the trapezium value is about 2313. That does not seem wrong to me, for 50 bins. If I use 5000 bins, it gets about 2498. Again, seems fine to me.

Regarding plotting ... you seem to only be plotting the exact value curve. So, it is impossible to provide any advice on why your code may be wrong.

If you want to plot the individual rectangles, I would suggest using the bar() function.

[beth]

ahh okay i see thank you. however I am still having a problem with my graph as it won't give me the stairs and steps function, it only fills the area underneath the graph. Any ideas?

[the cyclist]

I was running into a problem with mismatched vector lengths, but I did a quick fix and see what I think you want.

% Implentation of Rectangle, Trapezium and Simpson's Rule
% a: lower bound, b: upper bound, N: number of strips

clear,
%Default program
for run = menu('Run default program?','Yes','No');
      if run == 1;
          a = 0;
          b = 10;
          N = 50;
          h = (b-a)/N;   
          x = linspace(a,b, N+1);
          usedfx = 'x.^3';
          y=eval(usedfx);

        elseif run == 2;
            fprintf(' \n Numerical Integration');
            fprintf(' \n Input values for integral boundaries a and b with a<b')
            fprintf(' \n a = '); a= input('');
            fprintf(' \n b = '); b= input('');

            if a > b;
            fprintf(' \n a has to be less than the value of b')
            fprintf(' \n Input new values for a and b')
            fprintf(' \n a = '); a= input('');
            fprintf(' \n b = '); b= input('');
            end

            fprintf(' \n Input the number of strips')
            fprintf(' \n N = '); N= input('');
            if mod(N,2)> 0;
               N = N+1;
            else
            end

            x = linspace(a,b, N+1);
            fx = input('type in f(x) as a string ');
            usedfx = strrep(fx,'^','.^');
            h = (b-a)/N;
        end
   end

   % Rectangle Rule
   yr = eval(usedfx);
   rectsum = 0;
   for i = 1:N-1
       rectsum = rectsum + yr(i);
   end
   Rectangle_Rule = h*rectsum;

   % Trapezium Rule
   trapsum = 0;
  for i = 2:N-1;
       trapsum = trapsum + 2*y(i);
   end
   Trapezium_Rule = h/2*(y(1)+ trapsum + y(end));

   % Simpson's Rule
   simp4 = 0;
   simp2 = 0;

   for i = 1:2:N-1
       simp4 = simp4 + 4*y(i);
   end

   for i = 2:2:N-2
       simp2 = simp2 + 2*y(i);
   end
   Simpsons_Rule = h/3*(y(1)+ simp4 + simp2 + y(end));

   % Exact Value
   exact4 = 0;
   exact2 = 0;
   Ns = 50;
   hs = (b-a)/Ns;
   x = linspace(a,b, Ns+1);
   ys = eval(usedfx);
   for i = 1:2:Ns-1
       exact4 = exact4 + 4*ys(i);
   end

   for i = 2:2:Ns-2
       exact2 = exact2 + 2*ys(i);
   end
   Exact_Value = hs/3*(ys(1)+ exact4 + exact2 + ys(end));

   % Plotting the graph
   figure
   hold on
   plot(x,ys);
   stairs(x,yr,'r');
   steps(x,yr,'r');
   xlabel('x', 'FontWeight', 'bold');
   ylabel('fx', 'FontWeight', 'bold');
   title('Numerical integration of function fx' , 'FontWeight', 'bold');
   xstring = num2str(Exact_Value);
   text(6,200,['Exact Value = ' xstring]);
   xstring2 = num2str(Simpsons_Rule);
   text(6,150,['Simpsons Rule = ' xstring2]);
   xstring3 = num2str(Trapezium_Rule);
   text(6,100,['Trapezium Rule = ' xstring3]);
   xstring4 = num2str(Rectangle_Rule);
   text(6, 50,['Rectangle Rule = ' xstring4]);
   axis 'auto xy'

[beth]

thank you so much, the graph is looking a lot better. however when running the default program the value for the trapezium rule is higher than the value of the accurate result and simpsons rule result. any thoughts on why?

[beth]

also wouldn't y(i) not be recognised as we haven't defined y?