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

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

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

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