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

Question

Dominic 2013-3-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas

评论： 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 个评论
显示 4更早的评论隐藏 4更早的评论

Dominic 2013-3-29

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

Charles 2013-4-26

let's have a look at rest of the code

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Richard Treuren 2014-4-17

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#answer_133643

编辑：Richard Treuren 2014-4-17

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Richard Treuren 2014-4-17

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Answer 2

bym 2013-3-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#answer_80528

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Dominic 2013-3-31

在 MATLAB Online 中打开

Thanks

I'm not sure about some of the number. In this program, the user has to define the function, the lower limit of the integral, the upper limit of the integral, and the number of slices.

Below is my attempt to use your code with the variables I have set in my program.

% Variable list
% a    =  lower limit
% b    =  upper limit
% n    =  number of strips
% dx   =  strip width
% y_xs =  function (string form)
% y_x  =  function (see next line)
% y_xs = input('Enter your function: ', 's'); y_x = str2func(['@(x)', y_xs])
% x_smooth = array for plotting a smooth graph [a:dx/100:b]
% x_s  =  points of x for the curve
% y_s  =  points of y for the curve
% x_arcs = sample points 
% y_arcs = sample function
% 4. Create the strips (with parabolic arcs)
x_s = linspace(a,b);          % create x-data
y_s = y_xs(x_s);              % create y-data
x_arcs = linspace(a,b,n);     % sample points
y_arcs = y_xs(x_arcs);        % sample function
plot(x_arcs, y_arcs, 'b')     % plot sample points
xp = reshape(x_s, [], n)';    % prepare for plotting
xp( , )=x(end);               % duplicate end point
xp(   )=xp(     )             % duplicate end points
c(   )=0;                     % preallocate
for k = ?:?
    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

As you can see, some parts are empty as I'm not sure how can I fit the vraiables in there. I'm trying to figure it out now, but I would appreciate it a lot if you can help again as it will definitely save me a lot of time.

Thank you.

Dominic 2013-4-2

在 MATLAB Online 中打开

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 2013-4-2

please refer to the variable list above

请先登录，再进行评论。

Answer 3

Dominic 2013-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/68845-simpson-s-rule-illustration-how-to-create-those-parabolas#answer_80969

Bump. Please see comments above.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Dominic 2013-4-6

bump.

请先登录，再进行评论。

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（3 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（3 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论