Chapter 4 M-Files

M-File Home

Script File Index

 ShowNCIdea Displays the idea behind the Newton-Cotes rules. ShowQNC Computes Newton-Cotes estimates for specified integrands. ShowNCError Illustrates NCerror. ShowNCOpenError Illustrates NCOpenerror. ShowCompQNC Illustrates CompQNC on three examples. ShowAdapts Illustrates AdaptQNC. ShowQuads Illustrates quad and quad8. GLvsNC Compares Gauss-Legendre and Newton-Cotes rules. ShowSplineQ Illustrates SplineQ.

Function File Index

 NCWeights Constructs the Newton-Cotes weight vector. QNC The simple Newton-Cotes rule. NCError Error in the simple Newton-Cotes rule. NCOpenWeights Constructs the open Newton-Cotes weight vector. QNCOpen The simple open Newton-Cotes rule. NCOpenError Error in the simple open Newton-Cotes rule. CompQNC Equally-spaced, composite Newton-Cotes rule. AdaptQNC Adaptive Newton-Cotes quadrature. SpecHumps The humps function with function call counters. GLWeights Constructs the Gauss-Legendre weight vector. QGL The simple Gauss-Legendre rule. SplineQ Spline quadrature.

```% Script File: ShowNCIdea
% Illustrates the idea behind the Newton-Cotes approach

close all

% Generate a function and its interpolant.
x = linspace(0,2);
y = exp(-x).*sin(2*pi*x);
x0 = linspace(0,2,8);
y0 = exp(-x0).*sin(2*pi*x0);
c = InterpN(x0,y0);
pvals = HornerN(c,x0,x);

% Plot the function.
subplot(2,2,1)
plot(x,y,x,zeros(1,100),'.')
title('f(x)');
axis([0 2 -.5 1])

% Plot the interpolant.
subplot(2,2,2)
plot(x,y,x,pvals,x0,y0,'*',x,zeros(1,100),'.')
title('The Interpolant');
axis([0 2 -.5 1])

% Display the integral of the function.
subplot(2,2,3)
plot(x,y)
fill(x,y,[.5 .5 .5])
title('Integral of f(x)');
axis([0 2 -.5 1])

% Display the integral of the interpolant
subplot(2,2,4)
plot(x,pvals)
fill(x,pvals,[.5 .5 .5])
title('Integral of Interpolant');
axis([0 2 -.5 1])
```

```% Script File: ShowQNC
% Examines the closed Newton-Cotes rules.

while input('Another exammple? (1=yes, 0=no)')
fname = input('Enter within quotes the name of the integrand function:');
a = input('Enter left endpoint: ');
b = input('Enter right endpoint: ');
s = ['QNC(''' fname '''' sprintf(',%6.3f,%6.3f,m )',a,b)];
clc
disp(['  m      ' s])
disp(' ')
for m=2:11
numI = QNC(fname,a,b,m);
disp(sprintf(' %2.0f      %20.16f',m,numI))
end
end
```

```% Script File: ShowNCerror
% Examines the quality of the Newton-Cotes error bound.

clc
disp('Easy case: Integral from 0 to pi/2 of sin(x)')
disp(' ')
disp('Take DerBound = 1.')
disp(' ')
disp('   m           QNC(m)            Error      Error Bound')
disp(' ')
for m=2:11
numI = QNC('sin',0,pi/2,m);
err = abs(numI-1);
errBound = NCerror(0,pi/2,m,1);
s = sprintf('%20.16f   %10.3e    %10.3e',numI,err,errBound);
disp([ sprintf('  %2.0f   ',m) s])
end
```

```% Script File: ShowNCOpenError
% Examines the quality of the Newton-Cotes error bound.

clc
disp('Easy case: Integral from 0 to pi/2 of sin(x)')
disp(' ')
disp('Take DerBound = 1.')
disp(' ')
disp('   m           QNCOpen(m)        Error      Error Bound')
disp(' ')
for m=1:7
numI = QNCOpen('sin',0,pi/2,m);
err = abs(numI-1);
errBound = NCOpenError(0,pi/2,m,1);
s = sprintf('%20.16f   %10.3e    %10.3e',numI,err,errBound);
disp([ sprintf('  %2.0f   ',m) s])
end
```

```% Script File: ShowCompQNC
% Illustrates composite Newton-Cotes rules on three different
% integrands.

% Show QNC(m,n) errors for integral of sin from 0 to pi/2.
close all
figure
for m = [3 5 7]
% m-point rule used.
err = [];
for n = [1 2 4 8 16 32]
% n = number of subintervals.
err = [err  abs(compQNC('sin',0,pi/2,m,n) -1)+eps];
end
semilogy([1 2 4 8 16 32],err)
axis([0 40 10^(-17) 10^0])
text(33,err(6),sprintf('m = %1.0f',m))
hold on
end
title('Example 1.   QNC(m,n) error for integral of sin from 0 to pi/2')
xlabel('n = number of subintervals.')
ylabel('Error in QNC(m,n)')

% Show QNC(m,n) errors for integral of sqrt from 0 to 1.
figure
for m = [3 5 7]
% m-point rule used.
err = [];
for n = [1 2 4 8 16 32]
% n = number of subintervals.
err = [err abs(compQNC('sqrt',0,1,m,n) - (2/3))+eps];
end
semilogy([1 2 4 8 16 32],err)
axis([0 40 10^(-5) 10^(-1)])
text(33,err(6),sprintf('m = %1.0f',m))
hold on
end
title('Example 2.   QNC(m,n) error for integral of sqrt from 0 to 1')
xlabel('n = number of subintervals.')
ylabel('Error in QNC(m,n)')

% Show QNC(m,n) errors for integral of 1/(1+25x^2) from 0 to 1.
figure
for m = [3 5 7]
% m-point rule used.
err = [];
for n = [1 2 4 8 16 32]
% n = number of subintervals.
err = [err  abs(compQNC('Runge',0,1,m,n) - (atan(5)/5))+eps;];
end
semilogy([1 2 4 8 16 32],err)
axis([0 40 10^(-17) 10^0])
text(33,err(6),sprintf('m = %1.0f',m))
hold on
end
title('Example 3. QNC(m,n) error for integral of 1/(1+25x^2) from 0 to 1')
xlabel('n = number of subintervals.')
ylabel('Error in QNC(m,n)')
```

```% Script File: ShowAdapts
% Illustrates AdaptQNC for a range of tol values and a range of
% underlying NC rules. Uses the humps function for an integrand
% and displays information about the function evaluations.

global FunEvals VecFunEvals
close all
x = linspace(0,1,100);
y = humps(x);
plot(x,y)
u=[];
for tol = [.01  .001  .0001 .00001]
for m=3:2:9
figure
plot(x,y)
hold on
title(sprintf('m = %2.0f   tol = %10.6f',m,tol))
FunEvals = 0;
VecFunEvals = 0;
xLabel(sprintf('Scalar Evals = %3.0f    Vector Evals = %3.0f',FunEvals,VecFunEvals))
hold off
u = [u; FunEvals VecFunEvals];
end
end
```

```% Script File: ShowQuads
% humps function from 0 to 1.

clc
disp('------------------------------------------------------')
for tol = logspace(-2,-6,5)
s1 = sprintf('  %12.7f  %5.0f',Q1,count1);
s2 = sprintf('  %12.7f  %5.0f',Q2,count2);
disp([sprintf('%8.6f ',tol) s1 s2])
end

disp(sprintf('\n\n\n'))
disp('----------------------------')
for alpha = [-3 -2 -1];
for beta = 1:3
disp(sprintf('   %3.0f     %3.0f   %10.6f',alpha,beta,G))
end
end
```

```% Script File: GLvsNC
% Compares the m-point Newton-Cotes and Gauss-Legendre rules
% applied to the integral of sin(x) from 0 to pi/2.

clc
disp('Approximating the integral of sin from 0 to pi/2.')
disp(' ')
disp(' m           NC(m)                 GL(m) ')
disp('------------------------------------------------')
for m=2:6
NC = QNC('sin',0,pi/2,m);
GL = QGL('sin',0,pi/2,m);
disp(sprintf(' %1.0f  %20.16f   %20.16f',m,NC,GL))
end
```

ShowSplineQ

```% Script File: ShowSplineQ
% Examine spline quadrature on integral of sine from 0 to pi/2.

clc
disp('----------------------------')
for m=[5 50 500]
x    = linspace(0,pi/2,m);
y    = sin(x);
numI = SplineQ(x,y);
disp(sprintf(' %3.0f  %20.16f',m,numI))
end
```

NCWeights

```  function w = NCweights(m)
% w = NCweights(m)
%
% w is a column m-vector consisting of the weights for the m-point Newton-Cotes rule.
% m is an integer that satisfies 2 <= m <= 11.

if      m==2,   w=[1 1]'/2;
elseif  m==3,   w=[1 4 1]'/6;
elseif  m==4,   w=[1 3 3 1]'/8;
elseif  m==5,   w=[7 32 12 32 7]'/90;
elseif  m==6,   w=[19 75 50 50 75 19]'/288;
elseif  m==7,   w=[41 216 27 272 27 216 41]'/840;
elseif  m==8,   w=[751 3577 1323 2989 2989 1323 3577 751]'/17280;
elseif  m==9,   w=[989 5888 -928 10496 -4540 10496 -928 5888 989]'/28350;
elseif  m==10,  w=[2857 15741 1080 19344 5778 5778 19344 1080 15741 2857]'/89600;
else            w=[16067 106300 -48525 272400 -260550 427368 -260550 272400 -48525 106300 16067]'/598752;
end
```

QNC

```  function numI = QNC(fname,a,b,m)
% numI = QNC(fname,a,b,m)
%
% Integrates a function of the form f(x) named by the string fname from a to b.
% f must be defined on [a,b] and it must return a column vector if x is a column vector.
% m is an integer that satisfies 2 <= m <= 11.
% numI is the m-point Newton-Cotes approximation of the integral of f from a to b.

w = NCweights(m);
x = linspace(a,b,m)';
f = feval(fname,x);
numI = (b-a)*(w'*f);
```

NCError

```  function error = NCerror(a,b,m,M)
% error = NCerror(a,b,m,M)
%
% The error bound for the m-point Newton Cotes rule when applied to
% the integral from a to b of a function f(x). It is assumed that
% a<=b and 2<=m<=11. M is an upper bound for the (d+1)-st derivative of the
% function f(x) on [a,b] where d = m if m is odd, and m-1 if m is even.

if     m==2,   d=1;  c = -1/12;
elseif m==3,   d=3;  c = -1/90;
elseif m==4,   d=3;  c = -3/80;
elseif m==5,   d=5;  c = -8/945;
elseif m==6,   d=5;  c = -275/12096;
elseif m==7,   d=7;  c = -9/1400;
elseif m==8,   d=7;  c = -8183/518400;
elseif m==9,   d=9;  c = -2368/467775;
elseif m==10,  d=9;  c = -173/14620;
else           d=11; c = -1346350/326918592;
end
h = (b-a)/(m-1);
error = abs(c*M*h^(d+2));
```

NCOpenWeights

```  function w = NCOpenWeights(m)
% w = NCOpenWeights(m)
%
% w is a column m-vector consisting of the weights for the m-point open
% Newton-Cotes rule. m is an integer that satisfies 1 <= m <= 7.

if      m==1,   w = [1];
elseif  m==2,   w = [1 1]'/2;
elseif  m==3,   w = [2 -1 2]'/3;
elseif  m==4,   w = [11 1 1 11]'/24;
elseif  m==5,   w = [11 -14 26 -14 11]'/20;
elseif  m==6,   w = [611 -453 562 562 -453 611]'/1440;
else            w = [460 -954 2196 -2459 2196 -954 460]'/945;
end
```

QNCOpen

```  function numI = QNCOpen(fname,a,b,m)
% numI = QNC(fname,a,b,m,tol)
%
% Integrates a function of the form f(x) named by the string fname from a to b.
% f must be defined on [a,b] and it must return a column vector if x is a column vector.
% m is an integer that satisfies 1 <= m <= 9.
% numI is the m-point open Newton-Cotes approximation of the integral of f from a to b.

w = NCOpenWeights(m);
h = (b-a)/(m+1);
x = linspace(a+h,b-h,m)';
f = feval(fname,x);
numI = (b-a)*(w'*f);
```

NCOpenError

```  function error = NCOpenError(a,b,m,M)
% error = NCOpenError(a,b,m,M)
%
% The error bound for the m-point Newton Cotes rule when applied to
% the integral from a to b of a function f(x). It is assumed that
% a<=b and 1<=m<=7. M is an upper bound for the (d+1)-st derivative of the
% function f(x) on [a,b] where d = m if m is odd, and m-1 if m is even.

if     m==1,   d=1;  c = 1/3;
elseif m==2,   d=1;  c = 1/4;
elseif m==3,   d=3;  c = 28/90;
elseif m==4,   d=3;  c = 95/144;
elseif m==5,   d=5;  c = 41/140;
elseif m==6,   d=5;  c = 5257/8640;
else           d=7;  c = 3956/14175;
end
h = (b-a)/(m+1);
error = c*M*h^(d+2);
```

CompQNC

```  function numI = CompQNC(fname,a,b,m,n)
% numI = CompQNC(fname,a,b,m)
%
% Integrates a function of the form f(x) named by the string fname from a to b.
% f must be defined on [a,b] and it must return a column vector if x is a column vector.
% m is an integer that satisfies 2 <= m <= 11.
% numI is the composite m-point Newton-Cotes approximation of the integral of f
% from a to b with n equal length subintervals.

Delta = (b-a)/n;
h = Delta/(m-1);
x = a+h*(0:(n*(m-1)))';
w = NCWeights(m);
x = linspace(a,b,n*(m-1)+1)';
f = feval(fname,x);
numI = 0;
first = 1;
last = m;
for i=1:n
%Add in the inner product for the i-th subintegral.
numI = numI + w'*f(first:last);
first = last;
last = last+m-1;
end
numI = Delta*numI;
```

```  function numI = AdaptQNC(fname,a,b,m,tol)
%
% Integrates a function from a to b
% fname is a string that names an available function of the form f(x) that
% is defined on [a,b]. f should return a column vector if x is a column vector.
% a,b are real scalars, m is an integer that satisfies 2 <= m <=11, and
% tol is a positive real.
%
% numI is a composite m-point Newton-Cotes approximation of the
% integral of f(x) from a to b, with the abscissae chosen adaptively.

A1 = CompQNC(fname,a,b,m,1);
A2 = CompQNC(fname,a,b,m,2);
d = 2*floor((m-1)/2)+1;
error = (A2-A1)/(2^(d+1)-1);
if abs(error) <= tol
% A2 is acceptable
numI = A2+error;
else
% Sum suitably accurate left and right integrals
mid = (a+b)/2;
end
```

SpecHumps

```  function y = SpecHumps(x)
% y = SpecHumps(x)
%
% Yields humps(x) where x is an n-vector.
% FunEvals is an initialized global scalar that is increased by length(x).
% VecFunEvals is an initialized global scalar that is increased by 1.

global FunEvals VecFunEvals;
y = 1 ./ ((x-.3).^2 + .01) + 1 ./ ((x-.9).^2 + .04) - 6;
hold on
plot(x,y,'*')
hold off
FunEvals = FunEvals + length(x);
VecFunEvals = VecFunEvals + 1;
```

GLWeights

```  function [w,x] = GLWeights(m)
% [w,x] = GLWeights(m)
%
% w is a column m-vector consisting of the weights for the m-point Gauss-Legendre rule.
% x is a column m-vector consisting of the abscissae.
% m is an integer that satisfies 2 <= m <= 6.

w = ones(m,1);
x = ones(m,1);
if m==2
w(1) =  1.000000000000000; w(2) =  w(1);
x(1) = -0.577350269189626; x(2) = -x(1);
elseif m==3
w(1) =  0.555555555555558; w(3) =  w(1);
w(2) =  0.888888888888889;
x(1) = -0.774596669241483; x(3) = -x(1);
x(2) =  0.000000000000000;
elseif m==4
w(1) =  0.347854845137454; w(4) =  w(1);
w(2) =  0.652145154862546; w(3) =  w(2);
x(1) = -0.861136311594053; x(4) = -x(1);
x(2) = -0.339981043584856; x(3) = -x(2);
elseif m==5
w(1) =  0.236926885056189; w(5) =  w(1);
w(2) =  0.478628670499366; w(4) =  w(2);
w(3) =  0.568888888888889;
x(1) = -0.906179845938664; x(5) = -x(1);
x(2) = -0.538469310105683; x(4) = -x(2);
x(3) =  0.000000000000000;
else
w(1) =  0.171324492379170; w(6) =  w(1);
w(2) =  0.360761573048139; w(5) =  w(2);
w(3) =  0.467913934572691; w(4) =  w(3);
x(1) = -0.932469514203152; x(6) = -x(1);
x(2) = -0.661209386466265; x(5) = -x(2);
x(3) = -0.238619186083197; x(4) = -x(3);
end
```

```  function numI = QGL(fname,a,b,m)
% numI = QGL(fname,a,b,m,tol)
%
% Integrates a function from a to b
% fname is a string that names an available function of the form f(x) that
% is defined on [a,b]. f should return a column vector if x is a column vector.
% a,b are real scalars, m is an integer that satisfies 2 <= m <= 6.
%
% numI is the m-point Gauss-Legendre approximation of the
% integral of f(x) from a to b.

[w,x] = GLWeights(m);
fvals = feval(fname,((b-a)/2)*x + ((a+b)/2)*ones(m,1));
numI = ((b-a)/2)*w'*fvals;
```

SplineQ

```  function numI = SplineQ(x,y)
% numI = SplineQ(x,y)
%
% Integrates the spline interpolant of the data specified by the
% column n-vectors x and y. It is a assumed that x(1) < ... < x(n)
% and that the spline is produced by SPLINE. The integral is from
% x(1) to x(n).

S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
sum = 0;
for i=1:L
% Add in the integral from x(i) to x(i+1).
h = x(i+1)-x(i);
subI = h*(((rho(i,1)*h/4 + rho(i,2)/3)*h + rho(i,3)/2)*h + rho(i,4));
sum = sum + subI;
end
numI = sum;
```