Chapter 1: Problem Solutions
% 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)
% 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
% 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]))
% 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)
% 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
% 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')
% 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))
% 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))
% 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
% 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]')
% 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]')
% 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)
% 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)')
% Problem P1_3_2
%
% Roots of random quadratics.
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
% 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))
% 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))
% 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
% 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
% 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))
% 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))
% 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')
% 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)')
z
% 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))
Not Available.
Not Available
Not Available
Not Available
Not Available
Not Available
Not Available
% 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;
Not Available
Not Available
Not Available
Not Available
Not Available