Question

Write matlab program to compute ∫f(x)dx lower bound a upper bound b using simpson method and n points. Then, by recomputing with n/2 points and using richardson method, display the exact error and approximate error. (Test with a=0 b=pi f(x)=sin(x))

Answer 1

%%Matlab code for finding integration using 4 different method
clear all
close all
%functions for which integration have to do
fun1=@(x) sin(x);

%Number of steps in each methods
N=40;
%Limit of integration for function
a1=0; b1=pi;

fprintf('\n\nFor the function f(x)= ')
disp(fun1)
%Exact integration for function 1
I_ext = integral(fun1,a1,b1);
fprintf('\nExact integration value for a=%f to b=%f is %f.\n',a1,b1,I_ext)

%Simpson integration for function 1
I_smp1=simpson(fun1,a1,b1,N);
fprintf('\n\tSimpson integration for function 1 for N=%d is %f.\n',N,I_smp1)
fprintf('Error in Simpson integration for n=%d is %e.\n',N,abs(I_ext-I_smp1))

I_smp2=simpson(fun1,a1,b1,N/2);
fprintf('\n\tSimpson integration for function 1 for N=%d is %f.\n',N/2,I_smp2)
fprintf('Error in Simpson integration for n=%d is %e.\n',N/2,abs(I_ext-I_smp2))

%Integration value using Richardson extrapolation
I_rchrd=I_smp1+(I_smp1-I_smp2)/3;
fprintf('\n\tIntegration value using Richardson extrapolation is %f.\n',I_rchrd)
fprintf('Error in Richardson integration is %e.\n',abs(I_rchrd-I_ext))
%%Matlab function for Simpson integration
function val=simpson(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
    aa=func(xx(1));
    val=(dx/3)*(double(aa)+double(func(xx(end))));
    %loop for Riemann 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

%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%%%%