Chapter 4 M-Files
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
% to quadrature.
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;
num0 = AdaptQNC('SpecHumps',0,1,m,tol);
xLabel(sprintf('Scalar Evals = %3.0f Vector Evals = %3.0f',FunEvals,VecFunEvals))
hold off
u = [u; FunEvals VecFunEvals];
end
end
% Script File: ShowQuads
% Uses quad and quad8 to estimate the integral of the
% humps function from 0 to 1.
clc
disp('Tolerance quad N-quad quad8 N-quad8')
disp('------------------------------------------------------')
for tol = logspace(-2,-6,5)
[Q1,count1] = quad('humps',0,1,tol);
s1 = sprintf(' %12.7f %5.0f',Q1,count1);
[Q2,count2] = quad8('humps',0,1,tol);
s2 = sprintf(' %12.7f %5.0f',Q2,count2);
disp([sprintf('%8.6f ',tol) s1 s2])
end
disp(sprintf('\n\n\n'))
disp(' alpha beta quad')
disp('----------------------------')
for alpha = [-3 -2 -1];
for beta = 1:3
G = quad('F4_3_2',0,1,.001,0,alpha,beta);
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
% Script File: ShowSplineQ
% Examine spline quadrature on integral of sine from 0 to pi/2.
clc
disp(' m Spline Quadrature')
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
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
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);
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));
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
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);
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);
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) % 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; numI = AdaptQNC(fname,a,mid,m,tol/2) + AdaptQNC(fname,mid,b,m,tol/2); end
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;
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;
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;