Chapter 4: Problem Solutions
| P4.1.1 | P4.2.1 | P4.3.1 | P4.4.1 | P4.5.1 |
| P4.1.2 | P4.2.2 | P4.3.2 | P4.4.2 | P4.5.2 |
| P4.1.3 | P4.2.3 | P4.3.3 | P4.4.3 | |
| P4.1.4 | P4.2.4 | P4.3.4 | P4.4.4 | |
| P4.1.5 | P4.2.5 | P4.3.5 | P4.4.5 | |
| P4.3.6 | P4.4.6 | |||
| P4.3.7 | ||||
| P4.3.8 |
% Problem P4_1_1
%
% Error in the Corrected trapezoidal rule.
% Apply to the function f(x) = x^4.
clc
numI = CorrTrap('MyF411','DMyF411',0,1);
exact = 1/5;
err = abs(numI - exact);
% Use the equation
%
% error = c * (b-a)^5 * 24
%
% to solve for c:
c = err/24;
disp(sprintf('Computed constant = %18.15f',c))
disp(' Actual constant = 1/720')
function numI = CorrTrap(fname,fpname,a,b)
% fname is a string that names a function f(x)
% fpname is a string that names the derivative of f.
% a,b are real numbers
% numI is an approximation to the integral of f(x) from a to b.
fa = feval(fname,a);
fb = feval(fname,b);
Dfa = feval(fpname,a);
Dfb = feval(fpname,b);
h = (b-a);
numI = (h/2)*(fa + fb) + (h^2/12)*(Dfa - Dfb);
function y = MyF411(x)
y = x^4;
function y = DMyF411(x)
y = 4*x^3;
% Problem P4_1_2
%
% Computing the NC weights via a linear system solve
clc
disp(' m norm(NCweights(m)-MyNCweights(m)')
disp('-------------------------------------')
for m=2:11
disp(sprintf('%2.0f %20.16f',m,norm(NCweights(m)-MyNCweights(m))))
end
function w = MyNCweights(m);
% m a positive integer >= 2.
% w a column m-vector of weights associated with the
% m-point closed Newton-Cotes rule.
A = ones(m,m);
A(2,:) = linspace(0,1,m);
for i=3:m
A(i,:) = A(i-1,:).*A(2,:);
end
rhs = 1./(1:m)';
w = A\rhs;
% Problem P4_1_3
%
% NCweights(m) is exact for polynomials of degree m.
clc
disp(' m Error in NCweights(m) applied to x^m on [0,1]')
disp('----------------------------------------------')
for m=3:2:11
x = (linspace(0,1,m)').^m;
I = 1/(m+1);
numI = NCweights(m)'*x;
disp(sprintf('%2.0f %20.16f',m,abs(I-numI)))
end
% Problem P4_1_4
%
% Examines the quality of the Newton-Cotes error bound
global AVAL
AVAL = input('Integrand = 1/(1+a*x). Enter a = ');
clc
disp(sprintf('Integral from 0 to 1 of 1/(1+%2.0f*x)',AVAL))
disp(' ')
disp(' m QNC(m) Error Error Bound')
disp(' ')
for m=2:11
numI = QNC('NCTest',0,1,m);
err = abs(numI-log(AVAL+1)/AVAL);
if 2*floor(m/2) == m;
d = m-1;
else
d = m;
end
DerBound=1;
for k=1:d
DerBound = AVAL*k*DerBound;
end
errBound = NCerror(0,1,m,DerBound);
s = sprintf('%20.16f %10.3e %10.3e',numI,err,errBound);
disp([ sprintf(' %2.0f ',m) s])
end
% Problem P4_1_5
%
% Examines the quality of the Open 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=2: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
% Problem P4_2_1
%
% Error bound for composite NC rules applied to the integral
% of sin(x) from 0 to pi/2.
clc
disp(' m n=1 n=10 n=100')
disp('--------------------------------------')
for m=2:11
e1 = CompNCerror(0,pi/2,m,1,1);
e2 = CompNCerror(0,pi/2,m,1,10);
e3 = CompNCerror(0,pi/2,m,1,100);
disp(sprintf('%3.0f %8.3e %8.3e %8.3e',m,e1,e2,e3))
end
function error = CompNCerror(a,b,m,DerBound,n)
% a,b are real scalars
% m are positive integer with 2<=m<=11
% DerBound is a positive real
% n is apositive integer
%
% error is an upper bound for the error when the composite
% NC(m) rule is applied to integral of f(x) from a to b.
% Assumes that the d+1st derivative of f is bounded by DerBound
% where d = m-1 if m is even and d = m if m is odd.
if 2*floor(m/2)==m
d = m-1;
else
d = m;
end
error = NCerror(a,b,m,DerBound)/n^(d+1);
% Problem P4_2_2
%
% A more vectorized CompQNC
A = zeros(10,3);
for i=1:10
m = i+1;
for n=1:3
numI = CompQNC('sin',0,pi/2,m,n);
numI1 = CompQNC1('sin',0,pi/2,m,n);
A(i,n) = abs(numI-numI1);
end
end
clc
disp('Difference between CompQNC and CompQNC1.')
disp('Tested on integral from 0 to pi/2 of sin(x).')
disp(' ');
disp(' m n=1 n=2 n=3')
disp('-------------------------------------------')
for i=1:10
disp(sprintf('%2.0f %8.3e %8.3e %8.3e',i+1,A(i,:)))
end
function numI = CompQNC1(fname,a,b,m,n)
% 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
% n is a positive integer
%
% numI is the the composite m-point Newton-Cotes approximation of the
% integral of f(x) from a to b. The rule is applied on
% each of n equal subintervals of [a,b].
w = CompNCweights(m,n);
x = linspace(a,b,n*(m-1)+1)';
f = feval(fname,x);
numI = (b-a)*w'*f;
function w = CompNCweights(m,n)
% m and n are positive integers with 2<=m<=11.
% The weight vector associated with composite NC(m)
% rule with n equal length subintervals.
w = zeros(n*(m-1)+1,1);
for k=1:n
i=(1+(k-1)*(m-1)):(1+k*(m-1));
w(i) = w(i) + NCweights(m);
end
w = w/n;
% Problem P4_3_3
%
% 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
clc
home
for tol = [.01 .001 .0001 .00001]
for m=3:2:9
disp(' ')
disp(sprintf('tol = %8.5f m = %1.0f',tol,m ))
FunEvals = 0;
VecFunEvals = 0;
num = AdaptQNC('SpecHumps',0,1,m,tol);
disp(sprintf(' AdaptQNC: %4.0f %4.0f',FunEvals,VecFunEvals))
FunEvals = 0;
VecFunEvals = 0;
num0 = AdaptQNC1('SpecHumps',0,1,m,tol);
disp(sprintf(' AdaptQNC1: %4.0f %4.0f',FunEvals,VecFunEvals))
end
end
function numI = AdaptQNC1(fname,a,b,m,tol,OldEvals)
%
% 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
% tol is a positive real
% OldEvals is the value of f at linspace(a,b,m) (optional)
%
% numI is the composite m-point Newton-Cotes approximation of the
% integral of f(x) from a to b, with the abscissae chosen adaptively.
% Compute the integral two ways and estimate the error.
xvals = linspace(a,b,(2*m-1))';
if nargin == 5
OldEvals = feval(fname,xvals(1:2:(2*m-1)));
end
NewEvals = feval(fname,xvals(2:2:(2*m-2)));
fEvals = zeros(2*m-1,1);
fEvals(1:2:(2*m-1)) = OldEvals;
fEvals(2:2:(2*m-2)) = NewEvals;
A1 = (b-a)*NCweights(m)'*fEvals(1:2:(2*m-1));
A2 = ((b-a)/2)*NCweights(m)'*(fEvals(1:m) + fEvals(m:(2*m-1)));
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 = AdaptQNC1(fname,a,mid,m,tol/2,fEvals(1:m)) + ...
AdaptQNC1(fname,mid,b,m,tol/2,fEvals(m:2*m-1));
end
% Problem P4_2_4
%
% Illustrate the corrected trapezoid rule.
a = 0;
b = pi/2;
clc
disp('Integral of sin(x) from 0 to pi/2')
disp(' ')
disp(' n Trapezoidal Corrected Trapezoidal')
disp(' Rule Rule')
disp('----------------------------------------------')
for n=[1 2 4 8 16 32 64];
numI = CompQNC('sin',a,b,2,n);
numI1 = CompCorrTrap('sin','cos',a,b,n);
disp(sprintf('%3.0f %15.12f %15.12f',n,numI,numI1))
end
function numI = CompCorrTrap(fname,fpname,a,b,n)
% fname,fpname are strings that name a function f(x) and its derivative.
% a,b are reals
% n is a positive integer
%
% numI is an approximation of the integral of f from a to b.
h = (b-a)/n;
corr = (h^2/12)*(feval(fpname,a) - feval(fpname,b));
numI = CompQNC(fname,a,b,2,n) + corr;
z
% Problem P4_3_1
%
% Uses quad and Adapt to estimate the integral of the
% humps function from 0 to 1.
clc
disp(' Tol quad Adapt Quad Time / Adapt Time')
disp('-----------------------------------------------------------------')
for tol = logspace(-2,-5,4)
tic;
Q1 = quad('humps',0,1,tol);
t1 = toc;
tic;
Q2 = Adapt('humps',0,1,tol);
t2 = toc;
disp(sprintf('%7.5f %10.8f %10.8f %8.3f',tol,Q1,Q2,t2/t1))
end
% Problem P4_3_2
%
% Using the linearity of the quadrature problem.
tol = .001;
a = 0;
b = 1;
[numG,countG] = quad('humps',a,b,tol/2);
[numH,countH] = quad('sin',a,b,tol/2);
numDiff = (numG - numH)/2;
numSum = (numG + numH)/2;
[numGmH,countGmH] = quad('GmH',a,b,tol);
[numGpH,countGpH] = quad('GpH',a,b,tol);
count = countGmH + countGpH;
clc
home
disp('Method 1 computes the integrals of g and h separately and then combines.')
disp('Method 2 computes the integrals of (g-h)/2 and (g+h)/2 directly.')
disp(' ')
disp('g(x) = humps(x) is hard, h(x) = sin(x) = easy.')
disp(' ')
disp(' Idiff Isum g-evals h-evals')
disp('---------------------------------------------------')
disp(sprintf('Method 1: %8.4f %8.4f %4.0f %4.0f',numDiff,numSum,countG,countH))
disp(sprintf('Method 2: %8.4f %8.4f %4.0f %4.0f',numGmH,numGpH,count,count))
% Problem P4_3_3
%
% 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
clc
home
for tol = [.01 .001 .0001 .00001]
for m=3:2:9
disp(' ')
disp(sprintf('tol = %8.5f m = %1.0f',tol,m ))
FunEvals = 0;
VecFunEvals = 0;
num = AdaptQNC('SpecHumps',0,1,m,tol);
disp(sprintf(' AdaptQNC: %4.0f %4.0f',FunEvals,VecFunEvals))
FunEvals = 0;
VecFunEvals = 0;
num0 = AdaptQNC1('SpecHumps',0,1,m,tol);
disp(sprintf(' AdaptQNC1: %4.0f %4.0f',FunEvals,VecFunEvals))
end
end
function numI = AdaptQNC1(fname,a,b,m,tol,OldEvals)
%
% 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
% tol is a positive real
% OldEvals is the value of f at linspace(a,b,m) (optional)
%
% numI is the composite m-point Newton-Cotes approximation of the
% integral of f(x) from a to b, with the abscissae chosen adaptively.
% Compute the integral two ways and estimate the error.
xvals = linspace(a,b,(2*m-1))';
if nargin == 5
OldEvals = feval(fname,xvals(1:2:(2*m-1)));
end
NewEvals = feval(fname,xvals(2:2:(2*m-2)));
fEvals = zeros(2*m-1,1);
fEvals(1:2:(2*m-1)) = OldEvals;
fEvals(2:2:(2*m-2)) = NewEvals;
A1 = (b-a)*NCweights(m)'*fEvals(1:2:(2*m-1));
A2 = ((b-a)/2)*NCweights(m)'*(fEvals(1:m) + fEvals(m:(2*m-1)));
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 = AdaptQNC1(fname,a,mid,m,tol/2,fEvals(1:m)) + ...
AdaptQNC1(fname,mid,b,m,tol/2,fEvals(m:2*m-1));
end
z
z
(a) Solution:
1. I = Q_n+Ch^2; (by assumption)
2. So, I = Q_2n+Ch^2/4, which can be seen as follows.
Error composite rule = error simple rule/(r^(d+1)), where
r = ratio of points adding in, i.e., r=(2*n)/n=2,
d = m-1 (if m is even), m (if m is odd),
m = 2 for the trapezoid rule, so d=1.
So, r^(d+1)=2^(1+1)=4.
3. Subtracting 2 from 1 yields: abs(Q_n-Q_2n) = 3/4*Ch^2.
4. Solving 3 for Ch^2 yields Ch^2 = 4/3*abs(Q_n-Q_2n).
5. Therefore, abs(I-Q_2n) = Ch^2/4 (from 2)
= 4/3*abs(Q_n-Q_2n)*1/4 (from 4)
= 1/3*abs(Q_n-Q_2n) (simplifying).
% This script computes Q_2^(k+1) efficiently, where k is the
% smallest positive integer so that abs(I-Q_(2^(k+1))) < = tol.
% Establish some values for testing purposes.
a = 0; % left endpoint.
b = 2*pi; % right endpoint.
f = 'sin'; % function name.
tol = 10^(-3); % tolerance.
n=1; % number of points.
% Note: This script assumes the existence of a trapezoid function that
% can integrate a function f from a to b using the composite trapezoid
% rule.
% Compute Q_n and Q_2n initially.
Q_n = trapezoid(f,a,b,n);
Q_2n = trapezoid(f,a,b,2*n);
n=2;
while(1/2*abs(Q_n-Q_2n) > tol) % while abs(I-Q_2n) > tol
n = 2*n;
Q_n = Q_2n; % Take advantage of powers of 2 to avoid an extra trapezoid call.
Q_2n = trapezoid(f,a,b,n);
end;
z
z
% Problem P4_4_1
%
% Gauss Quadrature error bound applied to the integral of
% exp(cx) from a to b. (c>0)
% Em = ((b-a)^(2m+1))/(2m+1) * (m!)^4/((2m)!)^3
a = 0;
c = 1;
Mvals = zeros(1,3);
clc
disp(' b tol=10^-3 tol=10^-6 tol=10^-9 ')
disp('-----------------------------------------')
for b = 1:10
for j = 1:3
tol = .001^j;
m = 1;
E = (b-a)^3*b*exp(c*b)/24;
while E > tol
m = m+1;
E = E*(b-a)^2 * ((2*m-1)/(2*m+1)) * m * (m/((2*m-1)*2*m))^3;
end
Mvals(j) = m;
end
disp(sprintf('%2.0f %2.0f %2.0f %2.0f',b,Mvals))
end
% Problem P4_4_2
%
% Illustrates composite Gauss-Legendre rules on three different
% integrands.
%
% Show QNC(m,n) errors for integral of sin from 0 to pi/2.
close all
figure
for m = 2:6
% m-point rule used.
err = [];
for n = [1 2 4 8 16]
% n = number of subintervals.
err = [err abs(compQGL('sin',0,pi/2,m,n) -1)+eps];
end
semilogy([1 2 4 8 16],err)
axis([0 20 10^(-17) 10^0])
text(16,err(5),sprintf('m = %1.0f',m))
hold on
end
title('Example 1. QGL(m,n) error for integral of sin from 0 to pi/2')
xlabel('n = number of subintervals.')
ylabel('Error in QGL(m,n)')
% Show QGL(m,n) errors for integral of sqrt from 0 to 1.
figure
for m = 2:6
% m-point rule used.
err = [];
for n = [1 2 4 8 16 ]
% n = number of subintervals.
err = [err abs(compQGL('sqrt',0,1,m,n) - (2/3))+eps];
end
semilogy([1 2 4 8 16],err)
axis([0 20 10^(-6) 10^(-2)])
text(16,err(5),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 QGL(m,n)')
% Show QGL(m,n) errors for integral of 1/(1+25x^2) from 0 to 1.
figure
for m=2:6
% m-point rule used.
err = [];
for n = [1 2 4 8 16]
% n = number of subintervals.
err = [err abs(compQGL('Runge',0,1,m,n) - (atan(5)/5))+eps;];
end
semilogy([1 2 4 8 16],err)
axis([0 20 10^(-17) 10^0])
text(16,err(5),sprintf('m = %1.0f',m))
hold on
end
title('Example 3. QGL(m,n) error for integral of 1/(1+25x^2) from 0 to 1')
xlabel('n = number of subintervals.')
ylabel('Error in QGL(m,n)')
function numI = CompQGL(fname,a,b,m,n) % numI = CompQGL(fname,a,b,m,n) % % 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 % n is a positive integer % % numI is the composite m-point Gauss-Legendre approximation of the % integral of f(x) from a to b. The rule is applied on % each of n equal subintervals of [a,b]. z = linspace(a,b,n+1); % The easy approach which involves n calls to f numI=0; for i=1:n numI = numI + QGL(fname,z(i),z(i+1),m); end % An approach that involves a single call to f. Uses the fact that % % w'*f1 + ... + w'*f2 = [w;...;w]'*[f1;f2;...;fn] where w and the fi are m-vectors. [w,x] = GLweights(m); wtilde = []; xtilde = []; mdpt = (z(1:n)+z(2:n+1))/2; % The midpoints of the intervals. h_half = .5*(b-a)/n; % Half the subinterval length. for i=1:n wtilde = [wtilde;w]; xtilde = [xtilde;mdpt(i)+h_half*x]; end numI = h_half*(wtilde'*feval(fname,xtilde));
% Problem P4_4_3
%
% Test generalized SplineQ
clc
disp('Test on integral of x^3 - 3x^2 + 4x - 20 from a to b with spline')
disp('based on linspace(0,5,6).')
disp(' ')
disp(' ')
disp(' a b SplineQ1 Exact')
disp('-----------------------------------')
x = linspace(0,5,6)';
y = x.^3 - 3*x.^2 + 4*x - 20;
a = 0;
b = 5;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 1.5;
b = 5;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 1.5;
b = 3.5;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 2.3;
b = 2.6;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 5;
b = 5;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 3;
b = 4;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 0;
b = .3;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 0;
b = 0;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
a = 2;
b = 2;
numI = SplineQ1(x,y,a,b);
exact = (b^4/4 - b^3 + 2*b^2 -20*b ) - (a^4/4 - a^3 + 2*a^2 -20*a );
disp(sprintf('%4.2f %4.2f %10.6f %10.6f',a,b,numI,exact))
function numI = SplineQ1(x,y,a,b)
% x,y are n-vectors with x(1) < ... < x(n)
% a,b are such that x(1) <= a <= b <= x(n)
%
% numI is the integral from x(1) to x(n) of the not-a-knot spline
% interpolant of (x(i),y(i)), i=1:n
S = spline(x,y);
[x,rho,L,k] = unmkpp(S);
if a==x(L+1)
numI = 0;
return
end
ia = locate(x,a); %x(ia) <= a <= x(ia+1)
ib = locate(x,b); %x(ib) <= b <= x(ib+1)
if ia==ib
% same interval
hb = (b-x(ia));
fb = hb*(((rho(ia,1)*hb/4 + rho(ia,2)/3)*hb + rho(ia,3)/2)*hb + rho(ia,4));
ha = (a-x(ia));
fa = ha*(((rho(ia,1)*ha/4 + rho(ia,2)/3)*ha + rho(ia,3)/2)*ha + rho(ia,4));
numI = fb-fa;
return
end
% Start in the subinterval that houses a:
h1 = (x(ia+1)-x(ia));
f1 = h1*(((rho(ia,1)*h1/4 + rho(ia,2)/3)*h1 + rho(ia,3)/2)*h1 + rho(ia,4));
ha = (a-x(ia));
fa = ha*(((rho(ia,1)*ha/4 + rho(ia,2)/3)*ha + rho(ia,3)/2)*ha + rho(ia,4));
sum = f1-fa;
% Add in the contributions of the intervals in between a and b but not
% including either of those intervals.
for i=(ia+1):(ib-1)
% 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
% Add in the contribution of the rightmost interval that houses b.
hb = b-x(ib);
subI = hb*(((rho(ib,1)*hb/4 + rho(ib,2)/3)*hb + rho(ib,3)/2)*hb + rho(ib,4));
numI = sum + subI;
z
z
z
% Problem P4_5_1
%
% Scaling tol when adding integrals
a = 0;
b = 2;
tol = .0001;
m=5;
exact = quad8('humps',a,b,tol,.0001);
clc
home
disp('Error bounds are additive with AdaptQNC. Thus, if')
disp('an integral is expressed as the sum of n integrals and')
disp('each of those is computed to within tol, then the best we can')
disp('say about the sum is that it is correct to within n*tol.')
disp(' ')
disp('Illustrate with integral of humps from 0 to 2, tol = .0001 and')
disp('using AdaptQNC with m = 5. Total integral broken into n equal-length')
disp('subintegrals.')
disp(' ')
disp(' Subintegrals Error')
disp('------------------------------------')
for n=2:10
numI = 0;
x = linspace(a,b,n);
for i=1:n-1
d = AdaptQNC('humps',x(i),x(i+1),m,tol);
numI = numI +d;
end
disp(sprintf(' %2.0f %15.9e',n-1,abs(numI-exact)))
end
% Problem P4_5_2
%
% Compare Static Scheduling and Pool-of-task scheduling.
p = 4;
close all
figure
task_vec = [1000; 100*ones(39,1)];
[tr1,tr2] = Xtime(task_vec,p);
subplot(2,1,1)
plot(tr1)
title(sprintf('Static Assignment. Total Time = %5.0f',length(tr1)))
ylabel('Procs Busy')
xlabel('Time Step')
axis([0 length(tr1) 0 p+1])
subplot(2,1,2)
plot(tr2)
title(sprintf('Pool-of-Task Assignment. Total Time = %5.0f',length(tr2)))
ylabel('Procs Busy')
xlabel('Time Step')
axis([0 length(tr1) 0 p+1])
figure
task_vec = [100*ones(39,1);1000];
[tr1,tr2] = Xtime(task_vec,p);
subplot(2,1,1)
plot(tr1)
title(sprintf('Static Assignment. Total Time = %5.0f',length(tr1)))
ylabel('Procs Busy')
xlabel('Time Step')
axis([0 length(tr1) 0 p+1])
subplot(2,1,2)
plot(tr2)
title(sprintf('Pool-of-Task Assignment. Total Time = %5.0f',length(tr2)))
ylabel('Procs Busy')
xlabel('Time Step')
axis([0 length(tr1) 0 p+1])
function [trace_static,trace_pool] = XTime(task_vec,p)
% ask_vec is a column n-vector of positive integers
% p is the number of processors, a positive integer.
%
% trace_static is a column vector where trace_static(k) is the number of
% processors busy during time step k, where
% k = length(trace_static). Based on static
% wrap mapping assignment of tasks.
%
% trace_pool is a column vector where trace_pool(k) is the number of
% processors busy during time step k, where
% k = length(trace_pool). Based on the pool-of-task
% paradigm.
% Static assignment:
n = length(task_vec);
% How much work is each processor assigned?
proc_vec = zeros(p,1);
for k=1:p
proc_vec(k) = sum(task_vec(k:p:n));
end
t = zeros(max(proc_vec),1);
for k=1:p
t(1:proc_vec(k)) = t(1:proc_vec(k)) + ones(proc_vec(k),1);
end
trace_static = t;
% Pool-of-task assignment:
next = 1;
busy = zeros(p,1);
t = [];
while (next<=n) | any(busy>0)
% Make sure all free processors are given a task if possible.
[z,free] = min(busy);
while (z==0) & (next<=n)
% The first free processor is given a task
busy(free) = task_vec(next);
next = next + 1;
[z,free] = min(busy);
end
t = [t;sum(busy>0)];
busy = busy-1;
end
trace_pool = t;