Trapezoidal and Simpson’s Rules

Trapezoidal and Simpson’s Rules

COMPUTING FOR ENGINEERS

FOCUS
This lab focuses on the following items:
• Applying the trapezoidal and Simpson’s rules to determine the integral value both by hand andthrough MATLAB implementation
• Understanding the requirements and limitations of each integration method
Primary workshops involved:
• Workshop 17: Numerical integration – Trapezoidal rule
• Workshop 18: Numerical integration – Simpson’s rule

TASK 1

Consider the following function:

Evaluate f(x) for the following cases:
A. analytically
B. single application of the trapezoidal rule
C. composite trapezoidal rule using 4 segments
D. single application of the Simpson’s 1/3 rule
E. composite application of the Simpson’s 1/3 rule using 6 segments
F. Simpson’s 3/8 rule

TASK 2

Consider the following function:

Using the following parameters in your functions:
func: the function/equation that you are required to integrate
a, b: the integration limits
n: the number of points to be used for the integration
I: Integral estimate
A. Write a function capable of performing numerical integration of h(x) using the composite
trapezoidal rule
. Use your function to solve the equation.
B. Write a function capable of performing numerical integration of h(x) using the composite
Simpson’s 1/3 rule
. Use your function to solve the equation. 

TASK 3

An 11m beam is subjected to a load, and the shear force follows the equation

whereV is the shear force, and x is the length in distance along the beam. We know that V = dM/dx, andM is the bending moment. Integration yields the relationship

If M0 is zero and x=11, calculate M using
A. Analytical integration (pen and paper, then plug into MATLAB)
B. Composite trapezoidal rule using 13 points
C. Composite Simpson’s 1/3 rule using 13 points
Use fprintf to print the answers to each part. 

TASK 4
The force on a sailboat mast can be represented by the following function:

where z is the elevation above the deck and H is the height of the mast. The total force F exerted on themast can be determined by integrating this function over the height of the mast:

  1. Plot the force f against height z, assuming H=30.
  2. Use quad or integral functions in MATLAB to determine the force on the mast, assuming H=30.
  3. Determine the minimum number of segments required to achieve a percentage error of 0.01%or less when compared to the answer from part B, when using the composite trapezoidal rule.
  4. Repeat part C but for the composite Simpson’s 1/3 rule.

TASK 5

The cross-sectional area of a channel can be computed as

whereB is the total channel width (m), H is the depth (m), and y is the distance from the bank (m). In asimilar fashion, the average flow Q (m3/s) can be computed as

whereU is water velocity (m/s). The data for y, H and U are shown below.

0.5 1.3 1.25 1.8 1 0.25
0.03 0.06 0.05 0.13 0.11 0.02
  1. Plot H against Y, and U*H against Y on the same figure. Fit both sets of data with second-orderpolynomials and plot the fitted curves between y=0 and y=9 with increments of 0.1. Use thesepolynomials as function handles and use the composite trapezoidal rule to determine Ac and Q.Use 6 equally spaced points between the integral limits. Use fprintf to print the values of Ac andQ.
    B. Instead of using a second-order polynomial as a fit, you are to use the raw data. Integrate usingtrapezoidal segments. Use fprintf to print the values of Ac and Q.
    Hint: you can define polyval to a function. E.g. H_fn = @(x) polyval(p,x) 

Solution

TASK 1

comsimp.m 

function I=comsimp(func,a,b,n) %function to evaluate composite simpson’s integration

h=(b-a)/(n-1); %step size

x=a:h:b; %array of equally sized points

tp=0; %temp variable to store the summation

for j=2:2:size(x,2)-2

tp=tp+(4*func(x(j)))+(2*func(x(j+1)));

end

I=(h/3)*(tp+func(x(1))+func(x(2)));

end

comtrap.m 

function I=comtrap(func,a,b,n) %function to evaluate composite trapezoidal integration

h=(b-a)/(n-1); %step size

x=a:h:b; %array of equally sized points

tp=0; %temp variable to store the trapezoidal summation

for j=1:1:size(x,2)-1

tp=tp+func(x(j))+func(x(j+1));

end

I=(h/2)*tp;

end

comtrap2.m 

function I=comtrap2(x,y) %This function uses data sets instead of the function

tp=0; %temporary variable for summation

for j=1:1:size(y,2)-1

tp=tp+(0.5*(x(j+1)-x(j))*(y(j)+y(j+1)));

end

I=tp;

end 

ft2.m 

function y=ft2(x) %function file for task 2

y=1+exp(-x);

end

ft3.m 

function y=ft3(x) %function file for task 3

y=5+(0.25.*x.*x);

end 

ft4.m 

function y=ft4(x) %function file for task 4

y=200*(x/(5+x))*exp((-2*x)/30);

end 

lab9t2.m 

%Script to perform task 2

t2a=comtrap(@ft2,0,4,7); %implementing composite trapezoidal rule with 7 points

t2b=comsimp(@ft2,0,4,7); %implementing composite simpson’s rule with 7 points

fileID=fopen(‘lab9t2.txt’,’w’);

fprintf(fileID,’The value of the integral obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,t2a);

fprintf(fileID,’The value of the integral obtained by composite simpsons method using 7 points is %4.6f\r\n’,t2b);

fclose(fileID); 

lab9t3.m 

%Script for task 3

t3a=integral(@ft3,0,11); %Analytical Integral

t3b=comtrap(@ft3,0,11,13); %Composite Trapezoidal Integral with n=13

t3c=comsimp(@ft3,0,11,13); %Composite Simpson’s Integral with n=13

fileID=fopen(‘lab9t3.txt’,’w’);

fprintf(fileID,’The value of the integral obtained analytically is %4.6f\r\n’,t3a);

fprintf(fileID,’The value of the integral obtained by composite trapezoidal method using 13 points is %4.6f\r\n’,t3b);

fprintf(fileID,’The value of the integral obtained by composite simpsons method using 13 points is %4.6f\r\n’,t3c);

fclose(fileID); 

lab9t4.m 

%Script for task 4

z=0:30/100:30; %equally spaced points between z=0 and z=H=30

fz=ft4(z); %values of the function f(z)

figure(1)

plot(z,fz,’rx-‘),xlabel(‘elevation above the deck’),ylabel(‘force on the sailboat’),title(‘Force vs. Elevation’),grid on

saveas(gcf,’lab9t4.png’,’png’); %saves plot to file

t4b=integral(@ft4,0,30); %analytical integral

n1=3;

n2=3;

t4c1=comtrap(@ft4,0,30,n1); %initial trapezoidal integral with n=3

t4c2=comsimp(@ft4,0,30,n2); %initial simpson integral with n=3

while abs((t4c1-t4b)/t4b)>0.01 %maximum error criterion

n1=n1+1; %decreasing step size with each check

t4c1=comtrap(@ft4,0,30,n1);

end

while abs((t4c2-t4b)/t4b)>0.01 %maximum error criterion

n2=n2+1; %decreasing step size with each check

t4c1=comtrap(@ft4,0,30,n2);

end

fileID=fopen(‘lab9t4.txt’,’w’);

fprintf(fileID,’The minimum number of points required to achieve 0.01% error using trapezoidal integral is %d\r\n’,t4c1);

fprintf(fileID,’The minimum number of points required to achieve 0.01% error using simpson integral is %d\r\n’,t4c2);

fclose(fileID);

lab9t5.m 

%Script for task 5

y=[0 2 4 5 6 9];

H=[0.5 1.3 1.25 1.8 1 0.25];

U=[0.03 0.06 0.05 0.13 0.11 0.02];

Y=H.*U; %element-wise multiplication

figure(1)

plot(y,H,’rx-‘,y,Y,’ko-‘),xlabel(‘distance from the bank’),title(‘Data Plots’),legend(‘channel depth(H)’,’flow flux’),grid on

saveas(gcf,’lab9t5_1.png’,’png’);

t5a=polyfit(y,H,2); %second-order polynomial fit for H

t5b=polyfit(y,Y,2); %second-order polynomial fit for UH

t51=polyval(t5a,y); %evaluating polynomial fit

t52=polyval(t5b,y);

figure(2)

plot(y,H,’bo’,y,t51,’r-‘,y,Y,’bo’,y,t52,’k-‘),xlabel(‘distance from the bank’), title(‘Plots of polynomial fits for H(y) and U(y)H(y)’),grid on

saveas(gcf,’lab9t5_2.png’,’png’);

[email protected](x) t5a(1)*x*x+t5a(2)*x+t5a(3); %constructing polynomial splines

[email protected](x) t5b(1)*x*x+t5b(2)*x+t5b(3);

r5a=comtrap(ft51,0,9,7); %Composite Trapezoidal Integral on Polynomial Fit H with n=7

r5b=comtrap(ft52,0,9,7); %Composite Trapezoidal Integral on Polynomial Fit UH with n=7

r5c=comtrap2(y,H); %Integration by Trapezoidal Segments for H(y)

r5d=comtrap2(y,Y); %Integration by Trapezoidal Segments for U(y)H(y)

fileID=fopen(‘lab9t5.txt’,’w’);

fprintf(fileID,’The value of Ac obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,r5a);

fprintf(fileID,’The value of Q obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,r5b);

fprintf(fileID,’The value of Ac obtained by integrating trapezoidal segments is %4.6f\r\n’,r5c);

fprintf(fileID,’The value of Q obtained by integrating trapezoidal segments is %4.6f\r\n’,r5d);

fclose(fileID);