Question

2.4 2.

Modify the given code below which uses Simpson's 3/8 method with 6 subintervals to answer the question above.

if a==b
I=0
else
N=3;
h=(b-a)/N;

x=a:h:b;
y=Fun(x);
I=3*h/8*(y(1)+2*sum(y(4:3:(N-2)))+y(N+1));
N=2*N;
check=0;
% Calculating subsequent valus of I
while check==0;
h=(b-a)/N;
x=a:h:b;
y=Fun(x);
% Composite Simpson's 3/8 method, Eq. (9.22)
I_new=3*h/8*(y(1)+2*sum(y(4:3:(N-2)))+y(N+1));
I_new=I_new+3*h/8*3*(sum(y(2:3:(N-1)))+sum(y(3:3:N)));
% Compare solution with that calculated in the previous iteration
error=abs((I-I_new)/I)*100; % (*)
if error>0.1 % continue iteration
check=0;
N=N*2;
I=I_new;
elseif error<=0.1 % end iteration
check=1;
I=I_new;
end
end
end

Answer 1

------------------------simpson_three_eight_rule.m--------------------

function [simp_val] = simpson_three_eight_rule(y , h)

    % store the value of y3 + y6 + ... + y(n-3)

    sum3 = 0;

   

    % calculate the value of y3 + y6 + ... + y(n-3)

    for i = 4 : 3 : length(y) - 3

       

        sum3 = sum3 + y(i);

       

    end

   

    % store the value of y1 + y2 + y3 + ... + y(n-1)

    sum1 = 0;

    % calculate the value of y1 + y2 + y3 + ... + y(n-1)

    for i = 2 : length(y) - 1

       

        sum1 = sum1 + y(i);

       

    end

   

    % store the value of y1 + y2 + y4 + y5 + ... + y(n-1)

    sum2 = sum1 - sum3;

   

    simp_val = ( 3 * h / 8 ) * ( ( y(1) + y( length(y) ) ) + 3 * sum2 + 2 * sum3 );

   

end

----------------------main.m---------------------

%                               3h

% Simpson 1/3 rd Integration = -----[ ( y0 + yn ) + 3 * ( y1 + y2 + y4 + y5 + ... + y(n-1) ) + 2 * ( y3 + y6 + ... + y(n - 3) ) ) ]

%                                8

%               2x

% Integrate --------- in [0 , 2.4] with 6 intervals

%             1 + x^2

l_bound = 0;

r_bound = 2.4;

% store the number of intervals

n = 5;

% create a vector for x in [0 , 2.4] with 6 intervals

x = linspace( l_bound , r_bound , n + 1 );

% create a vector of size 1 x length(x) initilized to 0

y = zeros(1 , length(x));

% create a vector of values of ( sinx - logx + e^x )

for i = 1 : length(x)

   

   y(i) = 2 * x(i) / ( 1 + ( x(i) ^ 2 ) );

   

end

% the integration is fom 0 to 6 where no of interval is n

h = ( r_bound - l_bound ) / n;

simp_val = simpson_three_eight_rule( y , h );

fprintf('Simpson Integration : %f\n', simp_val);

Sample Output

Simpson Integration : 2.037146