Question

10 points (done by MATLAB:code submission required) It is known that -1 (a) Write a MATLAB code using composite Simpson 1/3 rule to compute this integral; (b) If a relative error of 1% is required, what should be the maximum step size h (or minimum number of subintervals to be used)? Run your code to find the answer for the maximum step size, and the numerical value of the integral.

Answer 1

Please find the code below with detailed inline comments.

CODE

=======================

function I = simpson3_f ( f )
f = @(x) sqrt(1 - x*x)
% for calculating integrals using Simpson's 1/3 rule when function is known

%asking for the range and desired accuracy

error = 1;
a = -1, b = 1;
%intial h and n
n = 100;
h = (b -a )/100;
%for calculating maximum of f''''(x) in the given region

for k = 0:100
x( 1, k+1 ) = a + k*h ;
y4( 1, k+1) = feval ( f, x(1,k+1) + 4*h ) - 4*feval( f, x(1,k+1) + 3*h )...
+ 6*feval( f, x(k+1) + 2*h ) - 4* feval ( f, x(1,k+1) + h ) ...
+ feval( f, x(k+1) );% fourth order difference
end
[ y i ] = max( y4);
x_opt = x(1,i);

% for calculating the desired value of h for 1% error
m=0;
ddf = feval ( f, x_opt + 4*h ) - 4*feval( f, x_opt + 3*h )...
+ 6*feval( f, x_opt + 2*h ) - 4* feval ( f, x_opt + h ) ...
+ feval( f, x_opt );% fourth order difference
% dff defined outside bracket just for convinence
while ddf * ( b -a )/180 > error % error for Simpson's 1/3 rule
m = m +1;
h = h*10^-m;
ddf = feval ( f, x_opt + 4*h ) - 4*feval( f, x_opt + 3*h )...
+ 6*feval( f, x_opt + 2*h ) - 4* feval ( f, x_opt + h ) ...
+ feval( f, x_opt );% defined here for looping
end

%calculating n
n = ceil( (b-a)/h );
if rem( n,2) == 0
n = n+1;
end
h = ( b-a )/n;
% calculating matrix X

for k = 1:(n+1)
X(k,1) = a + (k-1)*h;
X(k,2) =feval ( f, X(k,1));
end
%calculating integration

i= 1; I1 = 0;
while ( 2*i ) < (n+1)
I1 = I1 + X ( ( 2*i) , 2 );
i = i +1;
end
j=1; I2 =0;
while (2*j + 1) < (n+1)
I2 = I2 + X ( ( 2*j + 1) , 2);
j = j + 1;
end
I = h/3 * ( X( 1,2) + 4*I1 + 2*I2 + X(n,2) );% Simpson's 1/3 rule
% Display final result

h
n
disp(sprintf(' Using this integration for the given function in the range ( %10f , %10f ) is %10.6f .',a,b,I))

1 function 1= simpsons-f(f) 2 f e(x) sqrt(1 -x*x) % for calculating integrals using Simpson's 1/3 rule when function is known 4 5 %asking for the range and desired accuracy error = 1; 9 %intial h and n 10 n 100; 11 h-(b -a/100; 12 %for calculating maximum of f++"(x) in the given region 13 14 for k 0:100 15 16 17 18 19 end x( 1, k+1 ) = a + k*h ; + 6*feval( f, x(k+1) + 2*h4* feval ( f, x(1,k+1) + h) + feva ( f , x(k+1) );% fourth order difference 21 x opt x(1,i); 23 % for calculating the desired value of h for 1% error 24 m=0; 25 ddf feval Cf, x.opt 4*h -4*feval f, x opt + 3*h... 26 27 28 29 while ddf * ( b-a )/180 > error % error for Simpson's 1/3 rule 30 + fevol( f, X.opt );% fourth order difference % dff defined outside bracket just for convinence h-h*10A-m; 32 ddf = feval ( f, X-Opt + 4*h ) -4*feval( f, x-opt + 3*h ). .. + 6*feval( f, x.opt 2*h -4* feval ( f, x_opt + h + fevol( f, x. opt );% defined here for looping 34 35 end

36 37 %calculating n 38 n=ceil( (ba)/h); 39 if rem n,2)0 40 41 end 42 h=( b-a)/n; 43 % calculating matrix X 45 for k 1:(n+1) 46 47 48 end 49 %calculating integration 50 51 i- 1; 11 = 0; 52 while( 2#1 )< (n+1) 53 54 55 end 56 j-1; 12-0; 57 while (2j 1) < (n+1) 58 59 60 end 61 l=h/3*( X( 1,2) + 4*11 + 2*12 + xcn,2) );% Simpson's 1/3 rule 62 % Display final result 63 64 h 65 n 66 dispCsprintFC Using this integration for the given function in the range ( %10f , X10f ) is %10.6f . ',a,b,1)) 67 x(k,1) = a + (k-1)*h; X(k,2) -feval (f, X(k,1)); I2 12 X CC 2*j +1) , 2);

OUTPUT

======================

th Result Soctave -qf -no-window-system demo.m (x) sqrt (1 -xx) a-1 h= 0.019802 n= 101 Using this integration for the given function in the range ( -1.000000 , 1.000000 ) is ї.570466