# 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.

Simpson’s 1/3 rule

**TASK 3**

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

where*V *is the shear force, and x is the length in distance along the beam. We know that *V *= d*M*/*d*x, and*M *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:

- Plot the force f against height
*z*, assuming*H*=30. - Use quad or integral functions in MATLAB to determine the force on the mast, assuming
*H*=30. - 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.
- 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

where*B *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

where*U *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 |

- 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*and*Q*.

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