Chapter 1: Problem Solutions

Problems Home

 P1.1.1 P1.2.1 P1.3.1 P1.4.1 P1.5.1 P1.6.1 P1.7.1 P1.1.2 P1.2.2 P1.3.2 P1.4.2 P1.5.2 P1.6.2 P1.2.3 P1.3.3 P1.4.3 P1.5.3 P1.6.3 P1.2.4 P1.3.4 P1.4.4 P1.5.4 P1.2.5 P1.4.5 P1.5.5 P1.2.6 P1.4.6 P1.5.6 P1.2.7 P1.4.7 P1.5.7 P1.2.8 P1.4.8 P1.2.9 P1.4.9 P1.2.10 P1.4.10

P1.1.1

% Problem P1_1_1
%
% Plots various built-in functions

close all

x = linspace(-1,1);
plot(x,abs(x))
title('The Absolute Value Function')
pause(2)

figure
x = linspace(0,10);
plot(x,sqrt(x))
title('The Square Root Function')
pause(2)

figure
x = linspace(-2,2);
plot(x,exp(x))
title('The Exponential')
pause(2)

figure
x = linspace(1,100);
plot(x,log(x))
title('The Natural Log')
pause(2)

figure
x = linspace(0,4*pi);
plot(x,sin(x))
title('The Sine Function')
pause(2)

figure
x = linspace(0,4*pi);
plot(x,cos(x))
title('The Cosine Function')
pause(2)

figure
x = linspace(-1,1);
plot(x,asin(x))
title('The Inverse Sine Function')
pause(2)

figure
x = linspace(-1,1);
plot(x,acos(x))
title('The Inverse Cosine Function')
pause(2)

figure
x = linspace(-5,5);
plot(x,atan(x))
title('The Inverse Tangent Function')
pause(2)


P1.1.2

% Problem P1_1_2
%
% Plotting a periodic function.

close all

z1 = linspace(0,2,51);
y2 = sqrt(1 -(z1-1).^2);

% Note, plot(z1,y2) displays the first of the four semicircles.
z = [z1 z1(2:51)+2 z1(2:51)+4 z1(2:51)+6];
y = [y2 y2(2:51) y2(2:51) y2(2:51)];
plot(z,y)
axis equal


P1.2.1

% Problem P1_2_1
%
% Subvector illustration

clc
z = [10 40 20 80 30 70 60 90]
disp(' ')
disp('z(1:2:7) = ')
disp(' ')
disp(z(1:2:7))
disp(' ')
disp('z(7:-2:1) = ')
disp(' ')
disp(z(7:-2:1))

disp(' ')
disp('z([3 1 4 8 1]) = ')
disp(' ')
disp(z([3 1 4 8 1]))


P1.2.2

% Problem P1_2_2
%
% Subvector assignment illustration

clc
disp(' ')
disp('If z =')
z = [10 40 20 80 30 70 60 90];
disp(' ')
disp(z)
disp(' ')
disp('then after z(1:2:7) = zeros(1,4) it looks like')
disp(' ')
z(1:2:7) = zeros(1,4);
disp(z)

disp(' ')
disp('If z =')
z = [10 40 20 80 30 70 60 90];
disp(' ')
disp(z)
disp(' ')
disp('then after z(7:-2:1) = zeros(1,4) it looks like')
disp(' ')
z(7:-2:1) = zeros(1,4);
disp(z)

disp(' ')
disp('If z =')
z = [10 40 20 80 30 70 60 90];
disp(' ')
disp(z)
disp(' ')
disp('then after z([3 4 8 1]) = zeros(1,4) it looks like')
disp(' ')
z([3 4 8 1]) = zeros(1,4);
disp(z)


P1.2.3

% Problem P1_2_3
%
% Efficient circle plotting

close all

x = linspace(0,1,200);
y = sqrt(1 - x.^2);
plot(x,y)
axis square equal
hold on
plot(x,-y)
plot(-x,y)
plot(-x,-y)
hold off


P1.2.4

% Problem P1_2_4
%
% Superpositioning

close all

% Method 1:
figure
x = linspace(0,2*pi);
plot(x,sin(x),x,sin(2*x),x,sin(3*x),x,sin(4*x),x,sin(5*x))
title('Method 1')

% Method 2:
figure
plot(x,sin(x));
hold on
for k=2:5
plot(x,sin(k*x),'-')
end
hold off
title('Method 2')


P1.2.5

% Problem P1_2_5
%
% Plotting x, x^2, ... ,x^m where m is a solicited integer.

close all
m = input('Enter positive integer m: ');
x = linspace(0,1); xk = x;
plot(x,x);
hold on
for k=2:m
xk = xk.*x;
plot(x,xk)
end
hold off
title(sprintf(' Plots of x^k across [0,1] for k=1:%1.0f',m))


P1.2.6

% Problem P1_2_6
%
% Partial sum of the exponential power series.

close all

x = linspace(0,1);   % or any other vector
m = input('Enter m: ');

y = ones(size(x));
term = x;
for k=1:m
term = x.*term/k;
y = y+term;
end

plot(x,y)
title(sprintf('Partial sum of exp(x) using %2.0f terms',m))


P1.2.7

% Problem P1_2_7
%
% Ellipse plots.

close all

x = linspace(0,2*pi);
c = cos(x);
s = sin(x);
plot(3*6*c,-2+9*s,7+2*c,8+6*s)
axis square equal


P1.2.8

% Problem P1_2_8
%
% Sine superpositioning

close all

x = linspace(0,2*pi);
y = sin(x);
plot(x/2,y)
hold on
for k=1:3
plot((k*pi)+x/2,y)
end
hold off
title('sin(x) across [0,4*pi]')


P1.2.9

% Problem P1_2_9
%
% Plotting exp(-x)*sin(x)

close all

x = linspace(0,2*pi);
y = sin(x);
z = exp(-x);

plot(x,y.*z)
hold on
plot(x+2*pi,z(100)*(y.*z))
hold off
title('e^{-x}sin(x) across [0,4\pi]')


P1.2.10

% Problem P1_2_10
%
% Plot a function and its derivative

close all

subplot(2,1,1)
x = linspace(-10,10,200)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y = A*[2;3;7;5];
plot(x,y)
title('f(x) = 2sin(x) + 3sin(2x) + 7sin(3x)  + 5sin(4x)')

subplot(2,1,2)
Ap = [cos(x) 2*cos(2*x) 3*cos(3*x) 4*sin(4*x)];
yp = Ap*[2;3;7;5];
plot(x,yp)
title('f''(x) = 2cos(x) + 6cos(2x) + 21cos(3x)  + 20sin(4x)


P1.3.1

% Problem P1_3_1
%
% Explore the Updown Sequence

n = 200;
g = zeros(n,1);
for m=1:n
%  How long does it take for the UpDown sequence to converge
%  with starting value = m?
x = m;
k = 1;
while (x > 1)
if rem(x,2) == 0;
x = x/2;
else
x = 3*x+1;
end
k = k+1;
end
g(m) = k;
end
plot(g)
title('g(m) = is the index of the first iterate that equals 1')
xlabel('m = the starting value')
ylabel('g(m)')


P1.3.2

% Problem P1_3_2
%

clc
disp('       Prob Complex Roots  Prob Complex Roots')
disp('  n       (Uniform Coeff)     (Normal Coeff) ')
disp('--------------------------------------------------')
for n=100:100:800

% Uniformly distributed coefficients:
a = rand(n,3);
% The ith quadratic is a(i,1)x^2 + a(i,2)x + a(i,3). It has complex
% roots if the discriminant a(i,2)^2 - 4*a(i,1)*a(i,3)<0
pcmplx_U = sum(a(:,2).^2 < 4*a(:,1).*a(:,3))/n;

% Normally distributed coefficients:
a = randn(n,3);
% The ith quadratic is a(i,1)x^2 + a(i,2)x + a(i,3). It has complex
% roots if the discriminant a(i,2)^2 - 4*a(i,1)*a(i,3)<0
pcmplx_N = sum(a(:,2).^2 < 4*a(:,1).*a(:,3))/n;
disp(sprintf('%5d        %6.4f              %6.4f ',n,pcmplx_U,pcmplx_N))
end



P1.3.3

% Problem P1_3_3
%
% Volume of a 4-dimensional sphere

n = input('Enter number of trials: ');

% Throw n darts at the 4-cube with center (0,0,0,0) and edge length 2.
% This cube has volume 16.

A = -1 + 2*rand(n,4);

% The location of the ith dart is (A(i,1),A(i,2), A(i,3), A(i,4)).
% If the sum of squares of the coordinates is less than 1, it is inside
% the unit sphere.

hits = sum(A(:,1).^2 + A(:,2).^2 + A(:,3).^2 + A(:,4).^2 <= 1);
volume = (hits/n)*16;
clc
disp(sprintf('Volume of unit sphere in 4-space = %6.4f (based on %5d trials)',volume,n))


P1.3.4

% Problem P1_3_4
%
% Monte Carlo area estimate.

x0 = .2;
y0 = .4;

n = input('Enter number of trials: ');

x = -1+2*rand(n,1);
y = -1+2*rand(n,1);

% The points (x(i),y(i)) are uniformly distributed in S.

d = sqrt((x-x0).^2+(y-y0).^2);             % Distances to (x0,y0).
d_edge = min(min(abs(x+1),abs(x-1)),min(abs(y+1),abs(y-1)));   % Distances to nearest edge.

inside_S0 = sum(d<d_edge);
area_S0 = (inside_S0/n)*4;

clc
disp(sprintf('Area of S0 = %6.4f (based on %5d trials)',area_S0,n))


P1.4.1

% Problem P1_4_1
%
% Prints a table showing error in Stirling estimation of binomial coefficients.

close all
clc

disp( '  ')
disp('                           Stirling          Absolute    Relative')
disp('   k    52-choose-k      Approximation         Error       Error')
disp('----------------------------------------------------------------')
e=exp(1);
n = 52;
sn = sqrt(2*pi*n)*((n/e)^n);
true = n;
for k=2:13
sk   = sqrt(2*pi*k)*((k/e)^k);
snmk = sqrt(2*pi*(n-k))*(((n-k)/e)^(n-k));
bnk = (sn/sk)/snmk;
true = true*(n-k+1)/k;
abserror = abs(true-bnk);
relerror = abserror/true;
s1 = sprintf('  %2.0f %14.0f   %16.2f',k,true,bnk);
s2 = sprintf('   %13.2e    %5.2e',abserror,relerror);
disp([s1 s2])
end


P1.4.2

% Problem P1_4_2
%
% Plots, as a function of n, the relative error in the
% Taylor sin(x) approximation
%
%            x + x^2/2! +...+ x^n/n!

close all
nTerms = 50;
for x=[1 5 10 20 40]
figure
error = sin(x)*ones(nTerms,1);
s = x;
term = x;
x2 = x^2;
for k=1:nTerms
term = -x2*term/(2*k*(2*k+1));
s = s+ term;
error(k) = abs(error(k) - s);
end;
relerr = error/abs(sin(x));
semilogy(1:nTerms,relerr)
ylabel('Relative Error in Partial Sum.')
xlabel('Order of Partial Sum.')
title(sprintf('x = %5.2f',x))
pause(3)
end


P1.4.3

% Problem P1_4_3
%
% Taylor Approximation to sin(x)

close all
n = input('Enter number of terms:');
x = linspace(0,2*pi);
s = x;
term = x;
s2 = x.^2;
for k=1:n
term = -(term.*s2)/(2*k*(2*k+1));
s = s + term;
end
plot(x,sin(x),x,s)
title(sprintf('Approximation of sin(x) with %2.0f-term Taylor Sum',n))


P1.4.4

% Problem P1_4_4
%
% Exact representation of n!
% It is safe to assume that mantissa length will be the
% limiting factor because a factorial in the neighborhood
% of 2^100 will involve far more than 24 significant bits
% Compute MaxInt = the largest integer that can be stored in 24 bits.

clc
b = 24;
MaxInt = 2^b - 1;
n = 0;
nfact = 1;
s = MaxInt-1;

disp(' n              n!        x             n!/x        MaxInt-(n!/x)')
disp('------------------------------------------------------------------')
while(s>=0)
n = n+1;
nfact = nfact*n;
% Compute x = largest power of 2 that divides n!.

TwoPower = 1;
while floor((nfact/TwoPower))*TwoPower == nfact
TwoPower = 2*TwoPower;
end
x = TwoPower/2;

% n! is exactly representable if n!/x <= MaxInt

s = MaxInt-(nfact/x);
disp(sprintf('%2.0f %15.0f %8.0f  %15.0f    %15.0f',n,nfact,x,nfact/x,s))
end
disp(' ')
disp(sprintf('The largest n so n! is exactly representable is n = %1d.',n-1))


P1.4.5

% Problem P1_4_5
%
% Spacing of the floating point numbers.

clc
L = 4;
disp(' ')
disp('  Interval     Spacing of F.P. numbers')
disp('----------------------------------------')
for spacing=-12:-6
R = 2*L;
disp(sprintf('   [%3.0f,%3.0f)            2^(%3.0f)',L,R,spacing))
L = R;
end
disp(' ')
disp('So the next f.p. number after 70 is 70 + 2^-8')


P1.4.6

% Problem P1_4_6
%
% Spacing of the floating point numbers.

clc
L = 1;
disp(' ')
disp('  Interval     Spacing of F.P. numbers')
disp('----------------------------------------')
for k=1:6
R = 2*L;
disp(sprintf('   [%2.0f,%2.0f)            2^(-t+%1.0f)',L,R,k))
L = R;
end
disp(' ')
disp('So the minimum  y-x = 2^(-t+4) + 2^(-t+3)')


P1.4.7

z


P1.4.8

% Problem 1_4_8
%
% Nearest floating point number neighbor to a positive integer.

x = input('Enter a positive integer:');
t = input('Enter t, the length of the base-2 mantissa:');

% Compute m (.5<=m<1) and e so x = m * 2^e

m=x;e=0;
while (m>=1)
m=m/2; e=e+1;
end
while (m<1/2)
m=2*m; e=e-1;
end

clc
disp(sprintf('x = %1d',x))
disp(sprintf('t = %1d',t))
if (m==1/2)
disp(sprintf('The previous floating point number is %1d - 2^%1d',x,-t+e-1))
else
disp(sprintf('The previous floating point number is %1d - 2^%1d',x,-t+e))
end
disp(sprintf(   'The  next  floating  point  number is %1d + 2^%1d',x,-t+e))



P1.4.9

Not Available.


P1.4.10

Not Available


P1.5.1

Not Available


P1.5.2

Not Available


P1.5.3

Not Available


P1.5.4

Not Available


P1.5.5

Not Available


P1.5.6

% Problem P1_5_6
%
% Vectorized sine-cosine computation.

% Suppose you have a sine and cosine table with entries for 0,1,...,90 degrees.
% How could you "expand" these tables so that they reported the sines and
% cosines for 0, 0.5, 1.0, 1.5, 2.0,..., 89.0, 89.5, and 90.0 degrees?
% This problem is about that.

clc
a = input('Enter a:');
b = input('Enter b with b>a:');
n = input('Enter integer n with n>=1:');

c = cos(linspace(a,b,n+1));
s = sin(linspace(a,b,n+1));

[cnew,snew] = F(c,s,a,b);

err_cnew = sum(abs(cnew - cos(linspace(a,b,2*n+1)')))
err_snew = sum(abs(snew - sin(linspace(a,b,2*n+1)')))

function [cnew,snew] = F(c,s,a,b)
% a and b are scalars with a<b. c and s are row (n+1)-vectors with the
% property that c = cos(linspace(a,b,n+1)) and s = sin(linspace(a,b,n+1)).
%
% cnew and snew are column (2n+1)-vectors with the property that
% c = cos(linspace(a,b,2n+1)) and s = sin(linspace(a,b,2n+1)).

% Establish cnew and snew:

n = length(c)-1;
cnew = zeros(2*n+1,1);
snew = zeros(2*n+1,1);

% If theta = linspace(a,b,2*n+1), then
%
%             c = cos(theta(1:2:2*n+1))
%             s = sin(theta(1:2:2*n+1)).
%
% Thus, the odd-indexed components of cnew and snew are "free":

cnew(1:2:2*n+1) = c';
snew(1:2:2*n+1) = s';

% If delta = (b-a)/(2*n), the spacing between the theta(k), then using the given
% trig identities we have
%
%       cos(theta(2*k)) = cos(theta(2*k-1))cos(delta) - sin(theta(2*k-1))*sin(delta)
%       sin(theta(2*k)) = sin(theta(2*k-1))cos(delta) + cos(theta(2*k-1))*sin(delta)
%
% for k=1:n.

delta = (b-a)/(2*n);
c_delta = cos(delta);
s_delta = sin(delta);

cnew(2:2:2*n) = cnew(1:2:2*n-1)*c_delta - snew(1:2:2*n-1)*s_delta;
snew(2:2:2*n) = snew(1:2:2*n-1)*c_delta + cnew(1:2:2*n-1)*s_delta;


P1.5.7

Not Available


P1.6.1

Not Available


P1.6.2

Not Available


P1.6.3

Not Available


P1.7.1

Not Available