MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply TodayTo resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

Asked by Dominic on 27 Mar 2013

I managed to create the rectangle and trapezium strips, but stuck on the parabola strips for Simpson's Rule like the one shown below.

In this code, the user has to input the strip width, function, limits,

**Here's my code for the RECTANGLE STRIPS**

% 7. Display figure clf, hold on; plot([a b], [0 0], 'k' ), plot([0 0], [min( y_x(x) ) max( y_x(x))], 'k') % This shows a black line for the x-axis (y=0) and the y-axis (x=0). xlabel('x', 'FontWeight', 'bold'), ylabel('y(x)', 'FontWeight', 'bold') title(['Numeric integration through Rectangle Rule of y(x)=' , y_xs , ' with ', num2str(n), ' slices ||| Result is ', num2str(S_r) '.'] , 'FontWeight', 'bold')

% 8. To create the rectangular strips for x=a:dx:(b-dx); y_x(x); left = x; right = x+dx; bottom = 0; top = y_x(x); X = [left left right right]; Y = [bottom top top bottom]; %to create the shape fill(X,Y, 'b', 'FaceAlpha', 0.3) end

% 9. Display a smooth line in the graph x = a:dx/100:b; plot(x, y_x(x), 'k')

**Here's my code for the TRAPEZIUM STRIPS:**

% 7. Display figure clf, hold on; plot(x, y_x(x), 'k--') plot([a b], [0 0], 'k'), plot([0 0], [min( y_x(x) ) max( y_x(x))], 'k') xlabel('x', 'FontWeight', 'bold'), ylabel('y(x)', 'FontWeight', 'bold') title(['Numeric integration through Trapezium Rule of y(x)=' , y_xs , ' with ', num2str(n), ' slices ||| Result is ', num2str(S_t) '.'] , 'FontWeight', 'bold')

% 8. Display Trapezium strips for x=a:dx:(b-dx); y_x(x); left = x; right = x+dx; bottom = 0; top1 = y_x(x); top2 = y_x(x+dx); X = [left left right right]; Y = [bottom top1 top2 bottom]; %to create the shape fill(X,Y, 'b', 'FaceAlpha', 0.3) end

% 9. Display a smooth line in the graph x = a:dx/100:b; plot(x, y_x(x), 'b')

Comments on how to optimise and improve brevity this code would also be appreciated! Cheers

*No products are associated with this question.*

Answer by bym on 31 Mar 2013

here is what I have done

clc;clear; close all x = linspace(0,4*pi); % create data f = sin(x); xs = linspace(0,4*pi,11); % sample points fs = sin(xs); % sample function plot(x,f,'g') % plot function hold on plot(xs,fs,'bo') % plot sample points xp = reshape(x,[],5)'; % prepare for plotting xp(5,21)=x(end); % duplicate end point xp(1:4,21)=xp(2:5,1); % duplicate end points

c(5,3)=0; %preallocate for k = 1:5 c(k,:) = polyfit(xs(2*k-1:2*k+1),fs(2*k-1:2*k+1),2); %fit coefficients fill([xp(k,1) xp(k,:) xp(k,end)],[0 polyval(c(k,:),xp(k,:)) 0] ... ,'c', 'FaceAlpha',.1)

end ylim([-1.25,1.25])

Show 1 older comment

Dominic on 2 Apr 2013

This is the figure that is showing up with mine...

Here is the code, I am not sure what to replace the numerical values I am asking above, so I just copied what you gave. Surprise, surprise, it didn't show as the way I want it to.

Here is my try

% 4. Create the strips (with parabolic arcs) x_s = linspace(a,b); % create x-data y_s = y_x(x_s); % create y-data x_arcs = linspace(a,b,n+1); % sample points y_arcs = y_x(x_arcs); % sample function plot(x_arcs, y_arcs, 'bo-') % plot sample points xp = reshape(x_s, [], 5)'; % !!! prepare for plotting xp(5,21)=x(end); % !!! duplicate end point xp(1:4,21)=xp(2:5,1); % !!! duplicate end points

c(5, 3)=0; % !!! preallocate for k = 1:5 c(k,:)=polyfit(x_arcs(2*k-1:2*k+1), y_arcs(2*k-1:2*k+1), 2); %fit coefficients fill([xp(k,1) xp(k,:) xp(k,end)], [0 polyval(c(k,:), xp(k,:)) 0], 'c', 'FaceAlpha', 0.1) end

Answer by Dominic on 4 Apr 2013

Bump. Please see comments above.

Answer by Richard Treuren on 17 Apr 2014

Edited by Richard Treuren on 17 Apr 2014

I changed the script of proecsm a bit so that it does work for different step sizes and some other changes in the area plot

clc;clear; close all steps = 5; % number of steps x = linspace(0,2*pi,steps*12); % create data xs = linspace(0,2*pi,2*steps+1); % sample points f = sin(x); fs = sin(xs); % sample function

c(steps,3)=0; for k = 1:steps c(k,:) = polyfit(xs(2*k-1:2*k+1),fs(2*k-1:2*k+1),2); %fit coefficients hold on z = linspace(xs(2*k-1),xs(2*k+1),12); y = c(k,1).*z.^2+c(k,2).*z+c(k,3); area(z,y,'FaceColor',[0.6,1,1]) end

hold on plot(x,f,'r','LineWidth',2) % plot function hold on plot(xs,fs,'bo','LineWidth',2) % plot sample points

Richard Treuren on 17 Apr 2014

if you want to calculate the area (using simpson's rule) you can add the next lines to the script:

sim=0; for i=1:steps sim = sim + (xs(2*i+1)-xs(2*i-1))/6*(fs(2*i-1)+4*fs(2*i)+fs(2*i+1)); simp(i)= sim; %variable for plotting end sim

## 6 Comments

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_139420

shouldn't your result be close to 0?

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_139705

Yes, the value shown there is the area between the graph and the axis. I initially thought that was what we were being asked for then, but I've removed the ABS from the integration equation so the new values are now close to zero (something like 1e17 for sin(x)).

Do you happen to know how to create those curves on Simpson's Rule?

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_139768

you basically have to divide your function into sampling points that are divisible by 3 and fit each set of 3 to parabola using polyfit or equivalent. I will try to post something later today

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_139806

Looking forward to read that. Cheers

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_139809

Let me know if you want to have a look at the other parts of this code.

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/68845#comment_145704

let's have a look at rest of the code