Question

Write a MATLAB function that uses Simpson’s 1/3 method ( to integrate any function) to estimate an even
more accurate value of the integral using Richardson extrapolation. The function
should have the following inputs: the function name, the integration range, and the
number of intervals. It should output as an argument the numerical value of the
integral. The function should stop with an error message if an inappropriate number
of intervals are requested.

Answer 1

ONLY SIMPSON:

Matlab Function file:(my_simps13.m)

function intVal=my_simps13(foo,a,b,N)
% h is step size
h=(b-a)/N;
X=zeros(1,N+1);
Y=zeros(1,N+1);
for i=1:N+1
X(i)=a+(i-1)*h;
Y(i)=foo(X(i));
end
if rem(N,2)~=0
error('No. of segments must be in multiples of 2');
end
% I=h/3{(y0+yn)+4*(y1+y3+y5+y7+.......yn-1)+2*(y2+y4+y6+y8+....yn-1)}
% As the indices of vector in Mat Lab starts from 1
% yk is termed as yk+1 Example: y0 is termed as y1 and so on....
intVal=0;
for i=2:2:N
intVal=intVal+(4*h/3*(Y(i)));
end
for i=3:2:N
intVal=intVal+(2*h/3*(Y(i)));
end
intVal=intVal+((h/3)*(Y(N+1)+Y(1)));
end

Screenshot:

Save the above program as my_simps13.m in the working directory.

Matlab Script file:

clear;
clc;
% lower limit
a=0;
% upper limit
b=10;
% no. of intervals
N=10;
% function is ananymous so @ is not required
foo=@(x) 4*x+x.^3;
intVal=simpson_int(foo,a,b,N)

Screenshot:

Save the above program in the working directory and execute it.

Note:

Yoy can copy the code in command window also, saving script file is not needed

Result:

RICHARDSON INCLUDING:

Matlab Function file:(my_simps13.m)

function [intVal_simp,intVal_Rich]=my_simps13(foo,a,b,N)
% h is step size
h=(b-a)/N;
X=zeros(1,N+1);
Y=zeros(1,N+1);
for i=1:N+1
X(i)=a+(i-1)*h;
Y(i)=foo(X(i));
end
if rem(N,2)~=0
error('No. of segments must be in multiples of 2');
end
% I=h/3{(y0+yn)+4*(y1+y3+y5+y7+.......yn-1)+2*(y2+y4+y6+y8+....yn-1)}
% As the indices of vector in Mat Lab starts from 1
% yk is termed as yk+1 Example: y0 is termed as y1 and so on....
intVal_simp=0;
for i=2:2:N
intVal_simp=intVal_simp+(4*h/3*(Y(i)));
end
for i=3:2:N
intVal_simp=intVal_simp+(2*h/3*(Y(i)));
end
intVal_simp=intVal_simp+((h/3)*(Y(N+1)+Y(1)));
% I=h/4{(y0+yn)+2*(y1+y3+y5+y7+.......yn-1)+2*(y2+y4+y6+y8+....yn-1)}
% As the indices of vector in Mat Lab starts from 1
% yk is termed as yk+1 Example: y0 is termed as y1 and so on....
intVal_Rich=0;
for i=2:2:N
intVal_Rich=intVal_Rich+(2*h/3*(Y(i)));
end
for i=3:2:N
intVal_Rich=intVal_Rich+(2*h/3*(Y(i)));
end
intVal_Rich=intVal_Rich+((h/4)*(Y(N+1)+Y(1)));
end

Screenshot:

Save the above program as my_simps13.m in the working directory.

Matlab Script file:

clear;
clc;
% lower limit
a=0;
% upper limit
b=10;
% no. of intervals
N=10;
% function is ananymous so @ is not required
foo=@(x) 4*x+x.^3;
% exact value is 2700
EV=2700;
[intVal_simp,intVal_Rich]=my_simps13(foo,a,b,N)
fprintf('Error in Simpson1/3 is %.2f percntage\n',abs((EV-intVal_simp)/EV)*100);
fprintf('Error in Richardson is %.2f percentage\n',abs((EV-intVal_Rich)/EV)*100);

Screenshot:

Save the above program in the working directory and execute it.

Note:

Yoy can copy the code in command window also, saving script file is not needed

Result:

I hope this will help you, please give me thumbs-up.