Chapter 4: Problem Solutions

Problems Home

 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

P4.1.1

```% 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;
```

P4.1.2

```% 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;
```

P4.1.3

```% 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
```

P4.1.4

```% 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
```

P4.1.5

```% 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
```

P4.2.1

```% 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);
```

P4.2.2

```% 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;
```

P4.2.3

```% 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;

FunEvals = 0;
VecFunEvals = 0;
end
end

%
% 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;
end
```

P4.2.4

```% 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;
```

P4.2.5

```z
```

P4.3.1

```% Problem P4_3_1
%
% humps function from 0 to 1.

clc

disp('-----------------------------------------------------------------')
for tol = logspace(-2,-5,4)

tic;
t1 = toc;

tic;
t2 = toc;

disp(sprintf('%7.5f     %10.8f   %10.8f          %8.3f',tol,Q1,Q2,t2/t1))
end
```

P4.3.2

```% Problem P4_3_2
%
% Using the linearity of the quadrature problem.

tol = .001;
a = 0;
b = 1;

numDiff = (numG - numH)/2;
numSum  = (numG + numH)/2;
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))
```

P4.3.3

```% 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;

FunEvals = 0;
VecFunEvals = 0;
end
end

%
% 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;
end
```

P4.3.4

```z
```

P4.3.5

```z
```

P4.3.6

```(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;

```

P4.3.7

```z
```

P4.3.8

```z
```

P4.4.1

```% 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
```

P4.4.2

```% 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));

```

P4.4.3

```% 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;
```

P4.4.4

``` z
```

P4.4.5

```z
```

P4.4.6

```z
```

P4.5.1

```% Problem P4_5_1
%
% Scaling tol when adding integrals

a = 0;
b = 2;
tol = .0001;
m=5;
clc
home
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
numI = numI +d;
end
disp(sprintf('  %2.0f               %15.9e',n-1,abs(numI-exact)))
end
```

P4.5.2

```% Problem P4_5_2
%
% Compare Static Scheduling and Pool-of-task scheduling.

p = 4;

close all

figure
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
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])

% 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

% Static assignment:

% How much work is each processor assigned?
proc_vec = zeros(p,1);
for k=1:p
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;

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