Chapter 8: Problem Solutions
% 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
% 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
% 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)))
% 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;
% 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;
% 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
% 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
% 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;
% 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;
Not Available
Not Available
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
% 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
% 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;
Not Available.
% 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));
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
% 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
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
Not Available.
% 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);
Not Available.
Not Avialable.