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

조회 수: 7(최근 30일)
Dominic 2013년 3월 27일
댓글: Richard Treuren 2014년 4월 17일
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표시숨기기 이전 댓글 수: 5
Charles 2013년 4월 26일
let's have a look at rest of the code

댓글을 달려면 로그인하십시오.

### 답변(3개)

Richard Treuren 2014년 4월 17일
편집: Richard Treuren 2014년 4월 17일
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표시숨기기 없음
Richard Treuren 2014년 4월 17일
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

댓글을 달려면 로그인하십시오.

bym 2013년 3월 31일
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표시숨기기 이전 댓글 수: 3
Dominic 2013년 4월 2일
please refer to the variable list above

댓글을 달려면 로그인하십시오.

Dominic 2013년 4월 4일
##### 댓글 수: 1표시숨기기 없음
Dominic 2013년 4월 6일
bump.

댓글을 달려면 로그인하십시오.

### 범주

Find more on Data Import from MATLAB in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by