Discover MakerZone

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

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Simpson's Rule Illustration - How to create those parabolas?

Asked by Dominic on 27 Mar 2013
Latest activity Commented on by Richard Treuren on 17 Apr 2014 at 10:59

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

6 Comments

Dominic on 29 Mar 2013

Looking forward to read that. Cheers

Dominic on 29 Mar 2013

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

Charles on 26 Apr 2013

let's have a look at rest of the code

Dominic

Tags

Products

No products are associated with this question.

3 Answers

Answer by proecsm 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])

4 Comments

Dominic on 1 Apr 2013

Bump.

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
Dominic on 2 Apr 2013

please refer to the variable list above

proecsm
Answer by Dominic on 4 Apr 2013

Bump. Please see comments above.

1 Comment

Dominic on 6 Apr 2013

bump.

Dominic
Answer by Richard Treuren on 17 Apr 2014 at 10:35
Edited by Richard Treuren on 17 Apr 2014 at 10:40

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

1 Comment

Richard Treuren on 17 Apr 2014 at 10:59

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
Richard Treuren

Contact us