%% Solutions to Practice Set B: Math, Graphics, and Programming

%%
clear all
close all
format short
set(0, 'DefaultAxesFontSize', 16)

%% 1.
%% (a) 

[X, Y] = meshgrid(-1:0.1:1, -1:0.1:1);
contour(X, Y, 3*Y + Y.^3 - X.^3, 'k')

%%
[X, Y] = meshgrid(-10:0.1:10, -10:0.1:10); 
contour(X, Y, 3*Y + Y.^3 - X.^3, 'k')

%%
% Here is a plot with more level curves.  (Note that while *|ezcountour|*
% can make the previous plot, it lacks the capability to adjust the number
% of curves plotted.)
[X, Y] = meshgrid(-10:0.1:10, -10:0.1:10); 
contour(X, Y, 3*Y + Y.^3 - X.^3, 30, 'k')  

%% (b) 
% Now we plot the level curve with value $5$.
[X, Y] = meshgrid(-5:0.1:5, -5:0.1:5); 
contour(X, Y, 3.*Y + Y.^3 - X.^3, [5 5], 'k')  

%% (c) 
% We note that $f(1, 1) = 0$, so the appropriate commands 
% to plot the level curve of $f$ through the point $(1, 1)$
% are:

[X, Y] = meshgrid(0:0.1:2, 0:0.1:2); 
contour(X, Y, Y.*log(X) + X.*log(Y), [0 0], 'k') 

%% 2.
% We find the derivatives of the given functions:  

syms x r   

%% (a)
diff((2*x - 1)/(x^2 + 1), x)

%%
simplify(ans)  

%% (b)
diff(sin(3*x^2 + 2), x)  

%% (c)
diff(x^r, x)   

%% 3.
% We compute the following integrals:

%% (a)
int(cos(x), x, 0, pi/2)    

%% (b)
int(sin(3*x)*sqrt(1 - cos(3*x)), x)

%%
diff(ans, x)  

%% (c)
int(x^2*sqrt(x + 3), x)

%%
diff(ans, x)

%%
simplify(ans)  

%% (d)
int(exp(-x^2), x, -Inf, Inf)  

%% 4.
%% (a)
format long; integral(@(x) exp(sin(x)), 0, pi)

%%
% This is an integral MATLAB can't do symbolically:
syms x; int(exp(sin(x)), x, 0, pi)

%%
% By asking the Symbolic Toolbox to evaluate the integral numerically,
% we get the same answer as above.
double(int(exp(sin(x)), x, 0, pi))

%% (b)
integral(@(x) sqrt(x.^3 + 1), 0, 1)

%%
% MATLAB can do this integral symbolically, though the exact answer is
% rather long and cryptic.  Here we convert it to a decimal:
vpa(int(sqrt(x^3 + 1), x, 0, 1), 20)

%%
% The numerical approximation is accurate to at least 15 decimal places.

%% (c)  
integral(@(x) exp(-x.^2), -Inf, Inf)

%%
% In problem 3(d), we saw that the exact integral is $\sqrt{\pi}$.
sqrt(pi)

%%
% The numerical approximation using *|integral|* is accurate to 9 decimal
% places.  We can get more accuracy by decreasing the relative error
% tolerance from its default value of *|1e-6|*.
integral(@(x) exp(-x.^2), -Inf, Inf, 'RelTol', 1e-13)

%% 5. 

%% (a)
limit((1 + cos(x))/(x + pi), x, -pi)  

%% (b)
limit(x^2*exp(-x), x, Inf)  

%% (c)
limit(1/(x - 1), x, 1, 'left')  

%% (d)
limit(sin(1/x), x, 0, 'right')

%%
% The answer *|NaN|* (for ``not a number'') means there is no single limiting
% value.  In fact, every real number in the interval between $-1$ and $+1$
% is a ``limit point'' of $\sin(1/x)$ as $x$ tends to zero.  You can see
% why from the graph of $\sin(1/x)$ on the interval $(0,1]$.
X = 0:0.001:1; plot(X, sin(1./X), 'k')

%% 6.

%% (a) 
syms k n; symsum(k^2, k, 0, n)   

%% (b) 
syms r; symsum(r^k, k, 0, n)  

%%
% The answer is a little easier to read in ``pretty" form.
pretty(ans)  

%% (c) 
syms x; symsum(x^k/factorial(k), k, 0, Inf) 

%% 7.

%% (a)
taylor(exp(x), 'Order', 7)  

%% (b)
pretty(taylor(sin(x), 'ExpansionPoint', 2))  

%% (c)
taylor(log(x), 'ExpansionPoint', 1, 'Order', 5)   

%% 8.

%% (a) 
syms x y; ezsurf(sin(x)*sin(y), [-3*pi 3*pi -3*pi 3*pi])  

%% (b) 
syms x y; ezsurf((x^2 + y^2)*cos(x^2 + y^2), [-1 1 -1 1])  

%% 9.
T = 0:0.01:1;  
for j = 0:16
    fill(4*cos(j*pi/8)+(1/2)*cos(2*pi*T), ...
        4*sin(j*pi/8)+(1/2)*sin(2*pi*T), 'r');
    axis equal; axis([-5 5 -5 5]);
    M(j+1) = getframe;
end  
figure(gcf), movie(M)  

%%
% Obviously, we can't show the movie in this book, but you see the last
% frame above.

%% 10.

%% (a) 
A1 = [3 4 5; 2 -3 7; 1 -6 1]; b = [2; -1; 3];  
format short; x = A1\b

%%
A1*x  

%% (b) 
A2 = [3 -9 8; 2 -3 7; 1 -6 1]; b = [2; -1; 3];
x = A2\b  

%% 
% The matrix *|A2|* is indeed singular, as we will see in the next problem.
% The answer given is a valid solution, though as we saw in Problem 4 of
% Practice Set A, it is not the only solution.
A2*x

%% (c) 
A3 = [1 3 -2 4; -2 3 4 -1; -4 -3 1 2; 2 3 -4 1];
b3 = [1; 1; 1; 1]; x = A3\b3

%%
A3*x 

%% (d) 
syms a b c d u v
A4 = [a b; c d]; A4\[u; v] 

%%
det(A4)

%%
% The determinant of the coefficient matrix is the denominator 
% in the solution.  So we get a valid solution only if the coefficient 
% matrix is non-singular.

%% 11.

%% (a) 
[rank(A1) rank(A2) rank(A3) rank(A4)]

%%
% For the fourth matrix, MATLAB implicitly assumes that $ad - bc$ is not 0.

%% (b) 
% Only the second matrix is singular.

%% (c) 
det(A1)
%%
inv(A1)
%%
det(A2)  
%%
inv(A2)
%%
% MATLAB gives an answer, but notice that the size of the entries is quite
% large, roughly the inverse of the size of the determinant.  The actual
% determinant of *|A2|* is 0, and *|A2|* does not in fact have an
% inverse.
det(sym(A2))
%%
inv(sym(A2))
%%
det(A3)
%%
inv(A3)
%%
det(A4)
%%
inv(A4)  

%% 12.

%% (a)
[U1, R1] = eig(A1)  
%%
A1*U1 - U1*R1  
%%
% This is essentially zero.  Notice the *|e-14|*.

[U2, R2] = eig(A2)
%%
A2*U2 - U2*R2  
%%
% Same comment as in (a).

[U3, R3] = eig(A3);
%%
max(abs(A3*U3 - U3*R3))
%%
% Ditto yet again. (Note that we suppressed the output of the
% first command,
% and took the maximum of the absolute values of the entries in
% each column of *|A3*U3 - U3*R3|*, in order
% to keep the data from scrolling off the screen.)

[U4, R4] = eig(A4);
pretty(U4)
%%
pretty(R4)
%%
simplify(A4*U4 - U4*R4)

%% (b)
A = [1 0 2; -1 0 4; -1 -1 5]; 
[U1, R1] = eig(A)  
%%
B = [5 2 -8; 3 6 -10; 3 3 -7];  
[U2, R2] = eig(B)  
%%
% We observe that the first and second columns of *|U1|* are negatives 
% of the corresponding columns of *|U2|*, and the third columns are
% identical.  Finally:

A*B - B*A  

%% 13.

%% (a)
% If we set $X_n$ to be the column matrix with entries $x_n$, $y_n$, and $z_n$,
% and $M$ the square matrix with entries
% $1, 1/4, 0; 0, 1/2, 0; 0, 1/4, 1$,
% then $X_{n+1} = MX_n$.

%% (b)
% We have $X_n = MX_{n-1} = M^2X_{n-2} = ... = M^nX_0$.

%% (c)
M = [1, 1/4, 0; 0, 1/2, 0; 0, 1/4, 1];
[U,R] = eig(M) 

%% (d)
% $M$ should be $URU^{-1}$.  Let's check:

M - U*R*inv(U) 
%%
% Since $R^n$ is the diagonal matrix with entries $1, 1, 1/2^n$,
% we know that $R_\infty$ is the diagonal matrix 
% with entries $1, 1, 0$.  Therefore, $M_\infty = UR_\infty U^{-1}$. So:

Minf = U*diag([1, 1, 0])*inv(U) 

%% (e)
syms x0 y0 z0; X0 = [x0; y0; z0]; Minf*X0 
%%
% Half of the mixed genotype migrates to the dominant 
% genotype and the other half of the mixed genotype 
% migrates to the recessive genotype.  These are added 
% to the two original pure types, whose proportions are preserved.

%% (f)
M^5*X0
%%
M^10*X0  

%% (g)
% With the suggested alternate model, 
% only the first three columns of the table are relevant, 
% the transition matrix *|M|* becomes *|M = [1 1/2 0; 0 1/2 1; 0 0 0]|*.
% You can compute that the eventual population 
% distribution is *|[1 0 0]|*, and is independent of the initial population.

%% 14.
% The following commands create a JPEG file of the French flag.
% The *|set|* command has no effect on the JPEG file; we include
% it only to turn off the axis ticks in the figure window.

A = zeros(200, 300, 3);
A(:, 1:100, 3) = ones(200, 100);
A(:, 101:200, 1) = ones(200, 100);
A(:, 101:200, 2) = ones(200, 100);
A(:, 101:200, 3) = ones(200, 100);
A(:, 201:300, 1) = ones(200, 100);
image(A); axis equal off
imwrite(A, 'tricolore.jpg')

%%
% We now redefine the color of the leftmost third of the array and
% create a JPEG file of the Italian flag.

A(:, 1:100, 3) = zeros(200, 100);
A(:, 1:100, 2) = ones(200, 100);
image(A); axis equal off
imwrite(A, 'italia.jpg')

%%
% Of course you will have to run the commands above yourself
% to see the flags in color.  Then you can run the commands
% below to see a movie transforming
% the French flag into the Italian flag.

clear M
for k = 0:10
    A(:, 1:100, 3) = (10-k)*ones(200, 100)/10;
    A(:, 1:100, 2) = k*ones(200, 100)/10;
    M(k+1) = im2frame(A);
end
figure(gcf), movie(M)

%% 15.
% Our solution and its output is below.  First, we set *|n|* to $500$ in order 
% to save typing in the following lines and make it easier to change this 
% value later.  Then we initialize a zero matrix *|A|* of 
% the appropriate size and begin a loop that successively defines each 
% row of the matrix.  Finally, we extract the maximum value from 
% the list of eigenvalues of *|A|*.

n = 500;
A = zeros(n);
for k = 1:n
    A(k,:) = 1./(k:(k + n - 1));
end
max(eig(A))

%%
% There is also a one-line solution using the MATLAB function *|hankel|*,
% which creates a matrix constant along anti-diagonals.
max(eig(hankel((1:500).^(-1), (500:999).^(-1))))

%% 16.
% Again we display below our solution and its output.  First, we define a 
% vector $t$ of values between $0$ and $2\pi$, in order to be able later
% to represent 
% circles parametrically as $x = r \cos t, y = r \sin t$.  Then, we clear 
% any previous figure that might exist and prepare to create the figure in 
% several steps.  Let's say that the red circle will have radius $1$; then the 
% first black ring should have inner radius $2$ and outer radius $3$, and thus 
% the tenth black ring should have inner radius $20$ and outer radius $21$.  
% We start drawing from the outside in because the idea is to fill the 
% largest circle in black, then fill the next largest circle in white 
% leaving only a ring of black, then fill the next largest circle in black 
% leaving a ring of white, etc.  The *|if|* statement tests true when *|r|* is 
% odd and false when it is even.  We stop the alternation of black and 
% white at a radius of $2$ in order to make the last circle red instead of 
% black, then we adjust the axes to make the circles appear round.

t = linspace(0, 2*pi, 100);
cla reset; hold on
for r = 21:-1:2
    if mod(r,2)
        fill(r*cos(t), r*sin(t), 'k')
    else
        fill(r*cos(t), r*sin(t), 'w')
    end
end
fill(cos(t), sin(t), 'r')
axis equal; hold off  


%% 17.
% Here are the contents of our solution M-file.

type mylcm

%%
% Here are some examples. First, cases in which it should work:

mylcm(4, 5, 6)

%%

mylcm([6 7 12 15])

%%
% Next, a case in which we expect an error:

%%
try
    mylcm(4.5, 6)
catch err
    disp(err.message)
end

%% 18.
% Here are the contents of our solution M-file.

type letcount

%%
% Here is the output when we run this M-file on itself.

letcount('letcount.m')  

%% 19.
% The following M-file performs the desired task:

type editrecent

%%
% By assigning the output of *|dir|* to a variable, we get a structure
% array that contains among other things the names of all the M-files
% in the current folder and their modification dates.  If the
% structure array has $0$ elements, we print an error message.
% Otherwise, we form the dates into a cell array of strings, and use
% *|datenum|* to convert the strings into numbers.  We find the index of
% the largest number, indicating the most recent date, and finally
% pass the name of the corresponding file to *|edit|*.