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

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

P4.3.2

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

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

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

P4.5.2

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