Newton Cotes Numerical integration

%Newton Cotes Numerical Integration
%minster
clear all
close all
a = 0; % limits of integration
b = 1;
n=100;
h=10^-6;
esp1=10^-12/h^2
esp2=10^-12/h^4
h1=10^-1;
h2=10^-4;
%h=10^-6;%eps/h^2
n=100; %n points for simpson
f = @(x) (cos(x.^2));
x = [a:0.01:b];

d2 = (-4.*x.*x.*cos(x.^2)) -2.*sin(x.^2);%2nd deriv

for i =1:101
fx(i) = f(x(i));
p1x(i) = (1/(b-a))*(f(b)*(x(i)-a) – f(a)*(x(i)-b)); %p1(x) inter poly
end

traparea= 0.5*h*(f(a) + f(b))
errorNC = abs((1/12)*((h)^3)*max(d2))

dx= h/n;
dx1=esp1/n;
dx2=esp2/n;
sum = 0;
plot(x,fx)
title(‘Function and First Degree Interpolating Polynomial’)
hold on
plot(x,p1x,’k’)

for i =1:16

dx = dx/10;
dx1=esp1/n;
dx2=esp2/n;
sum = 0;
for i=0:n, % Loop over nodes.
x = a + i*h;
fx = f(x); % Evaluate integrand at x.
if i==0 | i==n,
sum = sum + fx/6; % Add 1/6 * Value of integrand at endpoints.
else
sum = sum + fx/3; % Add 1/3 * Value of integrand at interior points.
end;

if i < n, xmid = x + .5*dx; fxmid = f(xmid); sum = sum + 2*fxmid/3; % Add 2/3 * Value of integrand at midpoints. end; end; end; simpson = sum*dx % Multiply final result by dx. simpson1 = sum*dx1 simpson2=sum*dx2 disp(' Table 1. h1=10^-3\n') disp(' Trapezoid Area error ') disp(' =========================================================== ') disp(sprintf(' %20.12e %28.12e %20.3e ',traparea,errorNC)) sprintf('\n') disp(' Table 2. h1=10^-3') disp(' Simpson error (simpson-q) ') disp(' =========================================================== ') disp(sprintf(' %20.12e %20.12e', simpson,abs(simpson-0.9045))) sprintf('\n') disp(' Table 3. eps/h^2 and eps/h^4 ') disp(' eps/h^2 eps/h^4 ') disp(' =========================================================== ') disp(sprintf(' %20.12e %28.12e %20.3e ',abs(simpson1-0.9045),abs(simpson2-0.9045))) q=quad('cos(x.^2)',0,1,[1.e-12 1.e-12]) figure(2) plot(x,d2,'k+') title('Second Derivative') hold on plot(x,simpson, 'ro') hold all