```Two Incomplete Matlab Codes for Numerical Integration
Related to a Computer Project

Note:  You need to modify the codes to get them work.

% integral1.m: comparison of numerical integration rules

a = 0; b = pi/2; M = 10;

for n = 1:M
h = (b-a)/n;

rl = 0; rm = 0;
for i = 1:n
rl = rl + sin(a+(i-1)*h);         % left endpont rectangle rule
rm = rm + sin( ? );               % midpoint rectangle rule
end
rl = h*rl; erl(n) = log(abs(rl-1));
rm = h*rm; erm(n) = log(abs(rm-1));

t = sin(a) + sin(b);                % trapezoid rule
for i = 1:n-1
t = t + ?;
end
t = 0.5*h*t; et(n) = log(abs(t-1));

s = sin(a) + sin(b);                % Simpson's rule
for i = 1:n-1
s = s + 2*sin(a+i*h);
end
for i = 1:n
s = s + ?*sin( ? );
end
s = s*h/6; es(n) = log(abs(s-1));

g2 = 0;                             % Two-point Gaussian quadrature
for i = 1:n
g2 = g2 + sin(a+(i-0.5)*h-0.5*h/sqrt(3)) ...
+ sin(a+(i-0.5)*h+0.5*h/sqrt(3));
end
g2 = 0.5*h*g2; eg2(n) = log(abs(g2-1));

disp(sprintf('%3.0f %12.8f %12.8f %12.8f %12.8f %12.8f',n,rl,rm,t,s,g2))
end

x = [];
for n = 1:M
x(n) = log(n);
end
plot(x, erl, 'o', x, erm, '--', x, et, '-', x, es, '-.', x, eg2, ':')
print error1.ps

```

```% integral2.m : comparison of the trapezoid rule, the extrapolated
%               trapezoid rule, and the Simpson rule.

a = 0; b = pi/2; M = 8;

for n = 1:M+1
k = 2^(n-1);
h = (b-a)/k;

t = sin(a) + sin(b);                     % trapezoid rule
for i = 1:k-1
t = t + 2*sin(a+i*h);
end
t = 0.5*h*t; tra(n) = t;

s = sin(a) + sin(b);                     % Simpson's rule
for i = 1:k-1
s = s + 2*sin(a+i*h);
end
for i = 1:k
s = s + 4*sin(a+(i-0.5)*h);
end
s = s*h/6; sim(n) = s;
end

for n = 1:M                                % Richardson extrapolation
etr(n) = ? ;
end

for n = 1:M
k = 2^(n-1);
disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,tra(n),etr(n),sim(n)))
end
disp(sprintf('%3.0f %20.16f %41.16f',256,tra(9),sim(9)))
disp(sprintf(''))

for n = 1:M
k = 2^(n-1);
eetr(n) = log(abs(etr(n)-1));
etra(n) = log(abs(tra(n)-1));
esim(n) = log(abs(sim(n)-1));
disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,etra(n),eetr(n),esim(n)))
end

x = [];
for n = 1:M
x(n) = log(2^(n-1));
end
plot(x, etra, ':', x, eetr, '-', x, esim, 'o')
print error2.ps

```