Question

The Simpson numerical method tries to obtain a result closer to the analytic solution with and independent variable x, using the Simpson method (both 1/3 and 3/8 methods). Use a MATLAB script. The program should ask for the user to submit the function he wishes to integrate and the values of the a & b limits. Additionally it should ask the number of elements that would divide the function( sub-intervals).

Answer 1

clear all
close all
%function for which integration have to do
func=@(x) (sin(x.^2))./(1+x);
fprintf('function for which integration have to do f(x)=')
disp(func)
%all limits
a=0; b=1;
%number of subintervals
N=200;
val=simpson13(func,a,b,N);
fprintf('Integration using Simpson 1/3 method for N=20 with interval [0 1] is %f\n',val)

%number of subintervals
N=200;
val=simpson38(func,a,b,N);
fprintf('Integration using Simpson 3/8 method for N=20 with interval [0 1] is %f\n',val)

%%Matlab function for Simpson integration
function val=simpson13(func,a,b,N)
    % func is the function for integration
    % a is the lower limit of integration
    % b is the upper limit of integration
    % N number of rectangles to be used
    %splits interval a to b into N+1 subintervals
    xx=linspace(a,b,N+1);
    dx=xx(2)-xx(1); %x interval
    val=(dx/3)*(double(func(xx(1)))+double(func(xx(end))));
    %loop for Simpson integration
        for i=2:length(xx)-1
            xx1=xx(i);
            if mod(i,2)==0
                val=val+(dx/3)*4*double(func(xx1));
            else
                val=val+(dx/3)*2*double(func(xx1));     
            end
        end
end

%%Matlab function for Simpson 3/8 Method
function val=simpson38(f,a,b,n)
    % f is the function for integration
    % a is the lower limit of integration
    % b is the upper limit of integration
    % n is the number of trapizoidal interval in [a,b]
  
    %splits interval a to b into n+1 subintervals
    xx=linspace(a,b,n+1);
    dx=(xx(2)-xx(1)); %x interval
    val=f(a)+f(b);
    %loop for simpson integration
        for i=2:n
            if mod(i-1,3)==0
                val=val+2*double(f(xx(i)));
            else
                val=val+3*double(f(xx(i)));
            end
        end
    %result using simpson integration method
    val=(3/8)*dx*val;
end
  
%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%