# MATLAB: Numerical Integration Formulae (Trapezoidal and Simpson’s rule) for Equidistant x co-ordinates.

Problem – 1:
Write programs to evaluate the numerical value of I by Trapezoidal, Simpson’s 1/3 rule and Simpson’s 3/8 rule and compare the result

Solution:
(a)  Input: for h = pi/8
clc
clear all
close all

f = @(x) sin(x);
a = 0;
b = pi/2;

h = pi/8;
n = (b-a)/h;

t = (f(a)+f(b));

for i = 1:n-1
t = t+2*f(a+i*h);
end
A_T = 0.5*h*t

s = (f(a)+f(b));
for j = 1:2:n-1
s = s + 4*f(a+j*h);
end
for j = 2:2:n-2
s = s + 2*f(a+j*h);
end
A_S = (h/3)*s

s1 = (f(a)+f(b));

for i = 1:3:n-2
s1=s1+3*f(a+i*h);
end
for i = 2:3:n-1
s1 = s1+3*f(a+i*h);
end
for i = 3:3:n-3
s1=s1+2*f(a+i*h);
end
A_S1 = (3*h/8)*s1
Output:
A_T     =          0.9871
A_S     =          1.0001
A_S1   =          0.6287
Input: for h = pi/4
clc
clear all
close all

f = @(x) sin(x);
a = 0;
b = pi/2;

h = pi/4;
n = (b-a)/h;

t = (f(a)+f(b));

for i = 1:n-1
t = t+2*f(a+i*h);
end
A_T = 0.5*h*t

s = (f(a)+f(b));
for j = 1:2:n-1
s = s + 4*f(a+j*h);
end
for j = 2:2:n-2
s = s + 2*f(a+j*h);
end
A_S = (h/3)*s

s1 = (f(a)+f(b));

for i = 1:3:n-2
s1=s1+3*f(a+i*h);
end
for i = 2:3:n-1
s1 = s1+3*f(a+i*h);
end
for i = 3:3:n-3
s1=s1+2*f(a+i*h);
end
A_S1 = (3*h/8)*s1

Output:
A_T = 0.9481
A_S = 1.0023
A_S1 = 0.2945

(b)  Input: for h = 0.5
clc
clear all
close all

f = @(x) 1/(1+x);
a = 0;
b = 1;

h = 0.5;
n = (b-a)/h;

t = (f(a)+f(b));

for i = 1:n-1
t = t+2*f(a+i*h);
end
A_T = 0.5*h*t

s = (f(a)+f(b));
for j = 1:2:n-1
s = s + 4*f(a+j*h);
end
for j = 2:2:n-2
s = s + 2*f(a+j*h);
end
A_S = (h/3)*s

s1 = (f(a)+f(b));

for i = 1:3:n-2
s1=s1+3*f(a+i*h);
end
for i = 2:3:n-1
s1 = s1+3*f(a+i*h);
end
for i = 3:3:n-3
s1=s1+2*f(a+i*h);
end
A_S1 = (3*h/8)*s1

Output:
A_T = 0.7083
A_S = 0.6944
A_S1 = 0.2813

Input: for h = 0.25
clc
clear all
close all

f = @(x) 1/(1+x);
a = 0;
b = 1;

h = 0.25;
n = (b-a)/h;

t = (f(a)+f(b));

for i = 1:n-1
t = t+2*f(a+i*h);
end
A_T = 0.5*h*t

s = (f(a)+f(b));
for j = 1:2:n-1
s = s + 4*f(a+j*h);
end
for j = 2:2:n-2
s = s + 2*f(a+j*h);
end
A_S = (h/3)*s

s1 = (f(a)+f(b));

for i = 1:3:n-2
s1=s1+3*f(a+i*h);
end
for i = 2:3:n-1
s1 = s1+3*f(a+i*h);
end
for i = 3:3:n-3
s1=s1+2*f(a+i*h);
end
A_S1 = (3*h/8)*s1

Output:
A_T = 0.6970
A_S = 0.6933
A_S1 = 0.5531
Input: for h = 0.125
clc
clear all
close all

f = @(x) 1/(1+x);
a = 0;
b = 1;

h = 0.125;
n = (b-a)/h;

t = (f(a)+f(b));

for i = 1:n-1
t = t+2*f(a+i*h);
end
A_T = 0.5*h*t

s = (f(a)+f(b));
for j = 1:2:n-1
s = s + 4*f(a+j*h);
end
for j = 2:2:n-2
s = s + 2*f(a+j*h);
end
A_S = (h/3)*s

s1 = (f(a)+f(b));

for i = 1:3:n-2
s1=s1+3*f(a+i*h);
end
for i = 2:3:n-1
s1 = s1+3*f(a+i*h);
end
for i = 3:3:n-3
s1=s1+2*f(a+i*h);
end
A_S1 = (3*h/8)*s1

Output:
A_T = 0.6941
A_S = 0.6932
A_S1 = 0.5563

Problem – 2:
Add command in the program of problem 1 as necessary to obtain an integration of a function y=f(x), from its set of data points given below:
 x x0=1 x1=2 x2=3 x3=4 x4=5 x5=6 x6=7 y y0=2 y1=5 y2=10 y3=17 y4=26 y5=37 y6=50
Solution:
clc
clear all
close all

x = [1 2 3 4 5 6 7];
y = [2 5 10 17 26 37 50];

a = x(1,1);
b = x(1,7);
n = b - a;
h = (b - a)/n;

t = 0.5*(y(1,1)+y(1,7));

for i = 2:n
t = t+y(1,i);
end
A_T = h*t

s1 = (y(1,1)+y(1,7));

for i = 2:2:n
s1 = s1+4*y(1,i);
end
for i = 3:2:n-1
s1 = s1 + 2 * y(1,i);
end
A_S1 = (h/3)*s1

s2 = (y(1,1)+y(1,7));
for i = 2:3:n-1
s2=s2+3*y(1,i);
end
for i = 3:3:n
s2 = s2+3*y(1,i);
end
for i = 4:3:n-2
s2=s2+2*y(1,i);
end
A_S2 = (3*h/8)*s2
Output:
A_T = 121
A_S1 = 120
A_S2 = 120