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

Previous Post Next Post