Chapter 8: Problem Solutions

Problems Home

P8.1.1 P8.2.1 P8.3.1 P8.4.1
P8.1.2 P8.2.2 P8.3.2 P8.4.2
P8.1.3 P8.2.3 P8.3.3 P8.4.3
P8.1.4 P8.2.4 P8.3.4 P8.4.4
P8.1.5 P8.2.5 P8.4.5
P8.1.6 P8.2.6 P8.4.6
P8.1.7 P8.2.7 P8.4.7
P8.1.8 P8.2.8 P8.4.8
P8.1.9 P8.2.9
P8.1.10
P8.1.11
P8.1.12
P8.1.13
P8.1.14
P8.1.15
P8.1.16

 


P8.1.1

% Problem P8_1_1
%
% Plots the error in the function Mysqrt with min max starting
% approximation.

close all

s = 17/48;
t = 2/3;
x = linspace(.25,1)';
y = sqrt(x);
z = s+t*x;
plot(x,y,x,z)
title('sqrt(x) and (17+32x)/48')
Avals = linspace(.25,1,300);
s = zeros(1,300);
for k=1:300
   s(k) = MySqrt1(Avals(k));
end
error = abs(s-sqrt(Avals))./sqrt(Avals);
figure
plot(Avals,error+eps*.01)
Title('Relative Error in the Function MySqrt1(A)')
xlabel('A');

  function x = Mysqrt1(A)
% 
% A is a non-negative real number.
% 
% x is the square root of A.

if A==0
   x = 0;
else
   TwoPower = 1;
   m = A;
   while m>=1
      m = m/4;
      TwoPower = 2*TwoPower;
   end
   while m < .25
      m = m*4;
      TwoPower = TwoPower/2;
   end;	  
   % sqrt(A) = sqrt(m)*TwoPower	  
   x = (17+32*m)/48;
   for k=1:4
      x = (x + (m/x))/2;
   end
   x = x*TwoPower;
end

P8.1.2

% Problem P8_1_2
%
% Plots the error in the function Mysqrt with least squares starting
% approximation.

close all
s = 10/27;
t = 88/135;
x = linspace(.25,1)';
y = sqrt(x);
z = s+t*x;
plot(x,y,x,z)
title('sqrt(x) and (50+88x)/135')
Avals = linspace(.25,1,300);
s = zeros(1,300);
for k=1:300
   s(k) = MySqrt2(Avals(k));
end
error = abs(s-sqrt(Avals))./sqrt(Avals);
figure
plot(Avals,error+eps*.01)
Title('Relative Error in the Function MySqrt2(A)')
xlabel('A')


  function x = Mysqrt2(A)
% 
% A is a non-negative real number.
% 
% x is the square root of A.

if A==0
   x = 0;
else
   TwoPower = 1;
   m = A;
   while m>=1
      m = m/4;
      TwoPower = 2*TwoPower;
   end
   while m < .25
      m = m*4;
      TwoPower = TwoPower/2;
   end;	  
   % sqrt(A) = sqrt(m)*TwoPower	  
   x = (50+88*m)/135;
   for k=1:4
      x = (x + (m/x))/2;
   end
   x = x*TwoPower;
end




P8.1.3

% Problem P8_1_3
%
% Confirms the error bound (8.1).

close all
x = linspace(.25,1,500)';
y = sqrt(x);
z = (1+2*x)/3;
e = abs(y-z);
plot(x,e)
title(sprintf('Max error = %8.5f',max(e)))

P8.1.4

% Problem P8_1_4
%
% Extension of MySqrt to matrices.

A = [ 0.00       100;...
      .01       900;...
      1.00         0;...
      .000001 40000];
   X = MySqrtMat(A);
clc 
disp('A =                           MySqrtMat(A) = ')
for i=1:4
   disp(sprintf('%12.6f  %12.6f  %22.3f    %12.3f',A(i,1),A(i,2),X(i,1),X(i,2)))
end
disp(' ')
disp('A random 10-by-7 example with rows 2 and 4 scaled.')
A = rand(10,7);
A(4,:) = .0001*A(4,:);
A(2,:) = 1000*A(2,:)
X = MySqrtMat(A);
Xexact = sqrt(A);
err = abs(X - Xexact)./Xexact

  function X = MysqrtMat(A)
% 
% A is an m-by-n matrix with  non-negative entries.
% 
% X(i,j) = sqrt(A(i,j)), i=1:m and j=1:n.

% Initialize arrays. (Make A a vector.)

[r,c] = size(A);
X = zeros(r,c);
xvec = zeros(r*c,1);
avec = zeros(r*c,1);
avec = A(:);
I = find(avec>0);
m = avec(I);
TwoPower = ones(size(m));

% Normalize nonzero entries to [.25,1).
idx = find(m>=1);
while ~isempty(idx)
   m(idx) = m(idx)/4;
   TwoPower(idx) = 2*TwoPower(idx);
   idx = find(m>=1);
end
idx = find(m<.25);
while ~isempty(idx)
   m(idx) = 4*m(idx);
   TwoPower(idx) = TwoPower(idx)/2;
   idx = find(m<.25);
end
% Newton iteration.
xvec(I) = (1 + 2*m)/3;
for k=1:4
   xvec(I) = (xvec(I) + m./xvec(I))/2;
end

% Scale the positive square roots and establish the final X array.
xvec(I) = xvec(I).*TwoPower;
X(:) = xvec;


P8.1.5

% Problem P8_1_5
%
% Illustrate a bad bisection implementation
% The midpoint of the initial interval is a root.

clc
x = Bisection1('atan',-1,1,.0000001);
disp(sprintf('Computed root of atan(x) in [-1,1] is %5.3f',x))
disp(sprintf('Exact    root of atan(x) in [-1,1] is %5.3f',0))


  function root = Bisection1(fname,a,b,delta)
% 
% fname is a string that names a continuous function  f(x) of 
%       a single variable.
% a,b   define an interval [a,b]
% f is continuous, f(a)f(b) < 0
% delta is a non-negative real.
%
% root is the midpoint of an interval [alpha,beta]
%        with the property that f(alpha)f(beta)<=0 and
%        |beta-alpha| <= delta+eps*max(|alpha|,|beta|)

fa = feval(fname,a); 
fb = feval(fname,b);
if fa*fb > 0
   disp('Initial interval is not bracketing.')
   return
end
if nargin==3
   delta = 0;
end
while abs(a-b) > delta+eps*max(abs(a),abs(b))
   mid = (a+b)/2;
   fmid = feval(fname,mid);    
   if fa*fmid<0
      % There is a root in [a,mid].
      b  = mid; 
      fb = fmid;
   else
      % There is a root in [mid,b].
      a  = mid; 
      fa = fmid;
   end  
end
root = (a+b)/2;

P8.1.6

% Problem P8_1_6
%
% Illustrate recursive bisection.

clc
x = rBisection('cos',0,pi,.000001);
disp(sprintf('pi = %9.6f',2*x))


  function root = rBisection(fname,a,b,delta,fa,fb)
% 
% fname is a string that names a continuous function  f(x) of 
%         a single variable.
% a,b define an interval [a,b].
% f is continuous, f(a)f(b) < 0.
% delta   non-negative real.
% fa,fb   the value of f at a and b (optional).
%
% root is the midpoint of an interval [alpha,beta]
%         with the property that f(alpha)f(beta)<=0 and
%         |beta-alpha| <= delta+eps*max(|alpha|,|beta|)

if nargin==4
   fa = feval(fname,a);
   fb = feval(fname,b);
end
if abs(a-b) > delta+eps*max(abs(a),abs(b))
   mid = (a+b)/2;
   fmid = feval(fname,mid);    
   if fa*fmid<=0
      % There is a root in [a,mid].
      b  = mid; 
      fb = fmid;
   else
      % There is a root in [mid,b].
      a  = mid; 
      fa = fmid;    end 
   root = rBisection(fname,a,b,delta,fa,fb);	   
else
   root = (a+b)/2;
end

P8.1.7

% Problem P8_1_7
%
% Roots of real cubics

clc 
disp('  b    c    d         Middle Root')
disp('-----------------------------------')
for b = -2:2
   for c = -2:2
      for d = -2:2
         r = MiddleRoot(1,b,c,d);
         if ~isempty(r)
            disp(sprintf('%3.0f  %3.0f  %3.0f    %20.16f',b,c,d,r))		   
         end
      end
   end
end

  function r = MiddleRoot(a,b,c,d)
%
% a,b,c,d are real scalars.
%
% If ax^3 + bx^2 + cx + d has 3 distinct real roots then
% r is the  the middle root. Otherwise, r is the empty matrix.

a1 = 3*a; b1 = 2*b; c1 = c;   % Derivative = a1*x^2 + b1*x + c1
discrim = b1^2 - 4*a1*c1;
if (discrim<=0 | a1 == 0)
   % Does not have 3 real roots.
   r = [];
   return
else
   L = (-b1 - sqrt(discrim))/(2*a1);
   R = (-b1 + sqrt(discrim))/(2*a1);	   
   fL = ((a*L + b)*L + c)*L + d;
   fR = ((a*R + b)*R + c)*R + d;
   if fL*fR > 0
      % Does not have 3 real roots because the critical points are
      % on the same side of the x-axis.
      r = [];
      return
   else 
      % The middle root is in between the two critical x-values.
      while abs(L-R) > eps*max(abs(L),abs(R))
         mid = (L+R)/2;
         fmid = ((a*mid + b)*mid + c)*mid + d;
         
         if fL*fmid<=0
            % There is a root in [L,mid].
            R  = mid; 
            fR = fmid;
         else
            % There is a root in [mid,R].
            L  = mid; 
            fL = fmid;
         end  
      end
      r = (L+R)/2;
   end
end

P8.1.8

% Problem P8_1_8
%
% Bisection with an function that requires a quadrature.

tol = 10^(-14);
a = 0;  fa = 0;
x = .5; fx = quad8('MyF818',0,x,tol);
b = 1;  fb = 1;
while abs(fx-.5)>tol   
   if fx > .5
      % There is a root in [a,x].
      b  = x; 
      fb = fx;			 
   else
      % There is a root in [x,b].
      a  = x; 
      fa = fx;
   end 	  
   newx = (a+b)/2; 
   fx = fx + quad8('MyF818',x,newx,tol);
   x = newx;
end
clc
disp(sprintf('x = %16.14f',x))


  function y = MyF818(x)
y = pi*sin(pi*x/2)/2;



P8.1.9

% Problem P8_1_9
%
% A bisection application with a function that involves
% solving a linear system.

clc

% Playing with MyF819 shows that if alfa = 1 the 2-norm of the
% solution is bigger than 1 and if alfa = 3 then the 2-norm of
% the solution is less than 1.

alfa = Bisection('MyF819',1,3,.000000000001);
disp(sprintf('alfa = %15.12f',alfa))

  function nx = MyF819(alfa)
%
% alfa is a real scalar
%
% nx equals x'*x -1  where x satisfies (A+alfaI)x = ones(8,1) and 
%             A = toeplitz([6 -4  1 0 0 0 0 0]).

x = toeplitz([6+alfa -4 1 0 0 0 0 0])\ones(8,1);
% (There are more efficient ways to solve this linear system.)
nx = x'*x-1;

P8.1.10

Not Available

P8.1.11

Not Available



P8.1.12

Not Available.

P8.1.13

Not Available.

P8.1.14

Not Available.

P8.1.15

Not Available.

P8.1.16

Not Available.

P8.2.1

% Problem P8_2_1
%
% Modified GraphSearch Environment
 MyGraphSearch('DistMercEarth',0,1000,50)
 
 
   function MyGraphSearch(fname,a,b,nfevals)
%
% fname is a string that names a function f(x) that is defined 
%       on the interval [a,b]. 
%
% Produces a set of plots of the function f(x). The user
% specifies the x-ranges by mouseclicks. The fourth argument
% s used to determine how the plots are saved. If Save is
% nonzero, then each plot is saved in a separate figure window.
% If Save is zero or if GraphSearch is called with just three
% arguments, then only the final plot is saved. [L,R] is the
% x-range of the final plot. The plots are based on nfevals
% function evaluations. 

close all

%  Click in the search intervals.
x = linspace(a,b,nfevals);
y = feval(fname,x);
plot(x,y)
AnotherInterval=1;
n=0;
L = min(x);
R = max(x);
d = (R-L)/10;
hold on
while AnotherInterval
   xlabel('Enter a point of interest. (Click off x-range to terminate.)')
   [x1,y1] = ginput(1);
   if (L<=x1) & (x1<=R)
      n = n+1;
      a(n) = x1-d;
      b(n) = x1+d;
   else
      AnotherInterval = 0;
   end       
end
hold off

close all
for k=1:n
   % Search the kth interval
   figure
   AnotherPlot = 1;
   L = a(k);
   R = b(k);
   while AnotherPlot
      x = linspace(L,R,nfevals);
      y = feval(fname,x);
      ymin = min(y);
      plot(x,y,[L R],[ymin ymin])
      title(['The Function ' fname])
      v = axis;
      v(1) = L;
      v(2) = R;
      v(3) = v(3)-(v(4)-v(3))/10;
      axis(v)
      xlabel('Enter New Left Endpoint. (Click off x-range to terminate.)')
      text(R,ymin,['  ' num2str(ymin)])
      [x1,y1] = ginput(1);
      if (x1<L) | (R<x1) 	  
         AnotherPlot=0;
      else
         xlabel('Enter New Right Endpoint. (Click off x-range to terminate.)')
         [x2,y2] = ginput(1);
         if (x2<L) | (R<x2) 
            AnotherPlot=0; 
         else
            L = x1;
            R = x2;
         end
      end
      xlabel(' ')
   end 
end 



P8.2.2

% Problem P8_2_2
%
% Max area triangle in ellipse

close all

% The optimum triangle is isosceles. Solution with vertical orientation:

tstar = Golden('MyF822A',pi/2,1.5*pi);
A = -MyF822A(tstar);
tvals = linspace(pi/2,1.5*pi);
Avals = -MyF822A(tvals);
plot(tvals,Avals,tstar,A,'*')
title(sprintf('tstar = %6.3f  Area = %6.3f',tstar,A))
figure
theta = linspace(0,2*pi);
x = 2*cos(theta);
y = 3*sin(theta);
plot(x,y,[0 2*cos(tstar) -2*cos(tstar) 0],[3 3*sin(tstar) 3*sin(tstar) 3])
axis square equal

% The optimum triangle is isosceles. Solution with horizontal orientation:

tstar = Golden('MyF822B',pi/2,1.5*pi);
A = -MyF822B(tstar);
tvals = linspace(pi/2,1.5*pi);
Avals = -MyF822B(tvals);
figure
plot(tvals,Avals,tstar,A,'*')
title(sprintf('tstar = %6.3f  Area = %6.3f',tstar,A))

figure
plot(x,y,[2 2*cos(tstar) 2*cos(tstar) 2],[0 3*sin(tstar) -3*sin(tstar) 0])
axis square equal




  function Area = MyF822A(t)
% t is a real scalar
%  
% Area is the negative of the area of the triangle with vertices
%
%              (0,3), (2cos(t),3sin(t)) , (-2cos(t),3sin(t))

x = 2*cos(t);
y = 3*sin(t);
base   = 2*abs(x);
height = 3-y;
Area   = -base.*height/2;


  function Area = MyF822B(t)
% t is a real scalar
%
% Area is the negative of the area of the triangle with vertices
%
%              (2,0), (2cos(t),3sin(t)) , (2cos(t),-3sin(t))

x = 2*cos(t);
y = 3*sin(t);
base   = 2*abs(y);
height = 2-x;
Area   = -base.*height/2;



P8.2.3

Not Available.

P8.2.4

% Problem 8_2_4
%
% Fitting a circle to nearly circular data.


% Generate some nearly circular data

randn('state',0);    % .This way we all get the same answer.
n = 100;
r_true = 10;
delta = .5;
theta = linspace(0,2*pi,n)';
x = r_true*cos(theta); x_noisey = x + delta*randn(n,1);
y = r_true*sin(theta); y_noisey = y + delta*randn(n,1);

% Compute the "best" fitting circle's radius.

r = fmin('F824',0,max(sqrt(x_noisey.^2+y_noisey.^2)),[1],x_noisey,y_noisey);

xfit = r*cos(linspace(0,2*pi,200))';
yfit = r*sin(linspace(0,2*pi,200))';
plot(xfit,yfit,'r',x_noisey,y_noisey,'o')
axis equal
title(sprintf('Optimum r = %6.3f',F824(r,x_noisey,y_noisey)))


  function d = F824(r,x,y)
% r is a scalar and x and y are column n-vectors
%
% d = |sqrt(x(1)^2 + y(1)^2)) - r | + ... + |sqrt(x(n)^2 + y(n)^2)) - r |
%
% Each term is the distance from (x(i),y(i)) to the circle centered at (0,0)
% with radius r.

d = sum(abs(sqrt(x.^2 + y.^2)-r));


P8.2.5

Not Available.

P8.2.6

Not Available.

P8.2.7

Not Available.

P8.2.8

Not Available.

P8.2.9

Not Available.

P8.3.1

% Change the linesearch fragment to this:

% Line Search
% Try to get L<=Lmax so either f(xc+(L/2)*sc) < f(xc) or
% f(xc+L*sc) < f(xc).

L = Lmax;
Halvings = 0;
fL2 = feval(fName,xc+L*sc,plist);
fL1 = feval(fName,xc+(L/2)*sc,plist);
fLmax = fL2;
while (fL1>=fc) & (fL2>=fc) & Halvings<=20
   L = L/2;
   Halvings = Halvings+1;
   fL2 = fL1;
   fL1 = feval(fName,xc+(L/2)*sc,plist);
end 

% Find the quadratic interpolant 
%             q(lambda) = a(1) + a(2)lambda + a(3)lambda^2
% of (0,fc), (L/2,fL1), (L,fL2). Use Vandermonde approach.

a = [1 0 0 ; 1 (L/2) (L/2)^2 ; 1 L L^2]\[fc; fL1; fL2];
% Find the lambda that minimizes  q on [0,Lmax].
Lcrit = -.5*a(2)/a(3);

if (a(3)>0) & (0<=Lcrit) & (Lcrit<=Lmax)
   % There is a local minimum in the interval
   lambda = Lcrit;      
else
   % Check endpoints
   if fc < fLmax
      lambda=0;
   else
      lambda=Lmax;
   end
end


P8.3.2

Not Available.

P8.3.3

Not Available.

P8.3.4

Not Available.

P8.4.1

Not Available.

P8.4.2

Not Available.

P8.4.3

Not Available.

P8.4.4

Not Available.

P8.4.5

Not Available.

P8.4.6

% Problem 8_4_6
% A nonlinear least squares fitting problem.

clc
randn('seed',10);
disp('    L      R        a(1)        a(2)     lambda(1)   lambda(2)')
disp('----------------------------------------------------------------')
for s = [0 1 2 20; 2 3 4 21]
   L = s(1); R = s(2);
   [t,fvals] = GenProb(10,L,R,.001);
   lambda_best = FMINS('F846',[0;-1],[],[],t,fvals);
   a_best = [exp(lambda_best(1)*t) exp(lambda_best(2)*t)]\fvals;
   disp(sprintf('%5.0f  %5.0f  %10.4f  %10.4f  %10.4f   %10.4f',L,R,a_best,lambda_best))
end

function [t,fvals] = GenProb(m,L,R,noise)
t = linspace(L,R,m)';
fvals = 3*exp(-.2*t) + 5*exp(-4*t) + noise*randn(m,1);


function d = F846(lambda,t,fvals)
% lambda is a column 2-vector and t and fvals are column m-vectors.

% Get the best column 2-vector a given lambda, i.e., fit the
% data as best you can given lambda.

E_lambda = [exp(lambda(1)*t) exp(lambda(2)*t)];
a_best = E_lambda\fvals;

% Compute the norm of the residual.
d = norm(E_lambda*a_best - fvals);

P8.4.7

Not Available.

P8.4.8

Not Avialable.