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

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])

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

Sign in to comment.

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 6 Comments

## bym (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_139420

## Dominic (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_139705

## bym (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_139768

## Dominic (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_139806

## Dominic (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_139809

## Charles (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#comment_145704

Sign in to comment.