%% Solutions to Practice Set C: Developing Your MATLAB Skills

%%

clear all

%% 1.

%% (a)
%%
% Consider the status of the account on the last day of each month.  
% At the end of the first month, the account has $M + M\times J = M(1 + J)$
% dollars.
% Then at the end of the second month the account contains 
% $[M(1 + J)](1 + J) = M(1 + J)^2$ dollars.  Similarly, at the end of $n$ months,
% the account will hold $M(1 + J)^n$ dollars.  So our formula is
%%
% 
% $$T = M(1 + J)^n.$$
% 

%% (b)
%% 
% Now we take $M = 0$ and $S$ dollars deposited monthly.  At the end of the
% first month, the account has $S + S\times J = S(1 + J)$ dollars.  $S$ dollars are 
% added to that sum the next day, and then at the end of the second month,
% the account contains $[S(1 + J) + S](1 + J) = S[(1 + J)^2 + (1 + J)]$ dollars.
% Similarly, at the end of $n$ months, the account will hold     
%%
% 
% $$S[(1 + J)^n + \cdots + (1 + J)]$$
% 
%%
% dollars.  We recognize the geometric series---with the leading term missing,
% so the amount $T$ in the account after $n$ months will equal
%%
% 
% $$T = S\left[\frac{(1 + J)^{n+1} - 1}{(1 + J) - 1} - 1\right] = S\bigl[((1 + J)^{n+1} - 1)/J - 1\bigr].$$
% 

%% (c)
% By combining the two models, it is clear that in an account with an 
% initial balance $M$ and monthly deposits $S$, the amount of money $T$ after $n$ 
% months is given by
%%
% 
% $$T = M(1 + J)^n + S[((1 + J)^{n+1} - 1)/J - 1].$$
%

%% (d)
% We are asked to solve the equation
% $(1 + J)^n  = 2$
% with the values $J = 0.05/12$ and $J = 0.1/12$.
syms n
months = double(solve((1 + (0.05)/12)^n == 2)); 
years = months/12  

%% 
months = double(solve((1 + (0.1)/12)^n == 2));
years = months/12  

%%
% If you double the interest rate, you roughly halve the time to achieve the goal.

%% (e)
format bank
1000000/((((1 + 0.08/12)^(35*12 + 1) - 1)/(0.08/12)) - 1)

%%
% You need to deposit $433.06 every month.

%% (f)
syms n
months = solve(1000000 == ...
    300*(((1 + 0.08/12)^(n+1) - 1)/(0.08/12) - 1)); 
years = double(months)/12 

%%
% You have to work nearly 5 more years.

%% (g)
% First, taking the whole bundle at once, after 20 years the $65,000 left 
% after taxes generates

option1 = 65000*(1 + (.05/12))^(12*20)  

%%
% The stash grows to about $176,322.  The second option yields

S = .8*(100000/240);  
option2 = S*((1/(.05/12))*(((1 + (.05/12))^241) - 1) - 1)  

%%
% You only accumulate $137,582 this way.  Taking the lump sum up front 
% is clearly the better strategy.

%% (h)
rates = [.13, .15, -.03, .05, .10, .13, .15, -.03, .05];  

clear T
for k = 0:4
    T = 50000;
    for j = 1:5 
        T = T*(1 + rates(k + j));
    end
    disp([k + 1,T])
end  

%%
% The results are all the same; you wind up with $72,795 regardless of 
% where you enter in the cycle, because the product
% $\prod_{1\leq j\leq 5} (1 + \textbf{\texttt{rates}}(j))$ 
% is independent of the order in which you place the factors.  If you put 
% the $50,000 in a bank account paying 8%, you get

50000*(1.08)^5  

%%
% that is $73,466---better than the market.  The market's volatility hurts 
% you compared to the bank's stability.  But of course that assumes you 
% can find a bank that will pay 8%.  Now let's see what happens with no 
% stash, but an annual investment instead.  The analysis is more subtle 
% here.  Set $S = 10,000$.  At the end of one year, the account contains 
% $S(1 + r_1)$; then at the end of the second year $(S(1 + r_1) + S)(1 + r_2)$, 
% where we have written $r_j$ for $\textbf{\texttt{rates}}(j)$.
% So at the end of five years, the 
% amount in the account will be the product of S and the number
%%
%
% $$\prod_{j\geq 1}(1 + r_j) + \prod_{j\geq 2}(1 + r_j) + \prod_{j\geq 3}
% (1 + r_j) + \prod_{j\geq 4}(1 + r_j) + (1 + r_5).$$

%%
% If you enter at a different year in the business cycle, the terms get 
% cycled appropriately.  So now we can compute

format short  
clear T TT  
for k = 0:4
    T = ones(1, 5);
    for j = 1:5 
        TT = 1;
        for l = j:5
            TT = TT*(1 + rates(k + l));
        end
        T(j) = TT;
    end
    disp([k + 1,sum(T)])
end  

%%
% Multiplying each of these by $10,000 gives the portfolio amounts for the 
% five scenarios.  Not surprisingly, all are less than what one obtains by 
% investing the original $50,000 all at once. But in this model it matters 
% where you enter the business cycle.  It's clearly best to start your 
% investment program when a recession is in force and end in a boom.  
% Incidentally, the bank model yields in this case

(1/.08)*(((1.08)^6) - 1) - 1  

%%
% which is better than some investment models and worse than others.

%% 2.

%% (a)
% We can simulate a year in Tony's career by giving him credit for a
% hit each time a random number between 0 and 1 is $\le 0.338$.
% The following program does exactly that for $n$ at bats.

type yearbattingaverage

%% 
yearbattingaverage(500)  

%% (b)
% Now let's write a function M-file that simulates a 20-year career.  
% As with the number of at bats in a year, we'll allow for a varying 
% length career.

type career

%% 
career(500, 20)  

%% (c)
%
% If we ran the simulation 100 times and took the average, 
% it would be very close to .338.

%% 3. 
%
% We let $w, x, y, z$ denote the number of residences canvassed 
% in the four cities Gotham, Metropolis, Oz, and River City, respectively.  
% Then the linear inequalities specified by the given data are as follows:
% Non-negative data: $w,x,y,z \ge 0$;
% Pamphlets: $w + x + y + z \le 50,000$;
% Travel cost: $0.5w + 0.5x + y + 2z \le 40,000$;
% Time available: $2w + 3x + y + 4z \le 18,000$;
% Preferences: $w \le x$, $x + y \le z$;
% Contributions: $w + 0.25x + 0.5y + 3z \ge 10,000$.
% 
% The quantity to be maximized is:
% Voter support: $0.6w + 0.6x + 0.5y + 0.3z$.

%% (a)
% This enables us to set up and solve the linear programming problem in 
% MATLAB as follows: 

f = [-0.6 -0.6 -0.5 -0.3];  
A = [1 1 1 1; 0.5 0.5 1 2; 2 3 1 4; 1 -1 0 0; 0 1 1 -1; ...
    -1 -0.25 -0.5 -3; -1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1]; 
b = [50000; 40000; 18000; 0; 0; -10000; 0; 0; 0; 0];
simlp(f, A, b)

%%
% Jane should canvass 1268 residences in each of Gotham and Metropolis, 
% 1317 residences in Oz, and 2585 residences in River City.

%% (b)
% If the allotment for time doubles, then:

b = [50000; 40000; 36000; 0; 0; -10000; 0; 0; 0; 0];  
simlp(f, A, b)  
%%
% Jane should canvass 4000 residences in each of Gotham, Metropolis, and 
% River City, and ignore Oz. 

%% (c)
% Finally, if in addition she needs to raise $20,000 in contributions, then

b = [50000; 40000; 36000; 0; 0; -20000; 0; 0; 0; 0];  
simlp(f, A, b)  

%%
% Jane needs to canvass 2537 residences in each of Gotham and Metropolis, 
% 2634 residences in Oz, and 5171 in River City.

%% 4. 
%
% We let $w, x, y, z$ denote the number of hours that Joe spends 
% with the quarterback, the running backs, the receivers, and the 
% linemen, respectively.  Then the linear inequalities specified by 
% the given data are as follows:
% Non-negative data: $w , x, y, z \ge 0$;
% Time available: $w + x + y + z \le 50$;
% Point production: $0.5w + 0.3x + 0.4y + 0.1z \ge 20$;
% Criticisms: $w + 2x + 3y + 0.5z \le 75$;
% Prima Donna status: $x = y$, $w \ge x + y$, $x \ge z$;
% 
% The quantity to be maximized is:
% Personal satisfaction: $0.2w + 0.4x + 0.3y + 0.6z$.

%% (a)
% This enables us to set up and solve the linear programming problem 
% in MATLAB as follows: 

f = [-0.2 -0.4 -0.3 -0.6];  
A = [1 1 1 1; -0.5 -0.3 -0.4 -0.1; 1 2 3 0.5; 0 -1 1 0; ...
    0 1 -1 0; -1 1 1 0; 0 -1 0 1; -1 0 0 0; 0 -1 0 0; ...
    0 0 -1 0; 0 0 0 -1];  
b = [50; -20; 75; 0; 0; 0; 0; 0; 0; 0; 0];
simlp(f, A, b)  

%%
% Joe should spend 9.3 hours each with the running backs and receivers; 
% 5.6 hours with the linemen; and the majority of his time, 25.9 hours, 
% with the quarterback.

%% (b)
% If the team only needs 15 points to win, then:

b = [50; -15; 75; 0; 0; 0; 0; 0; 0; 0; 0];  
simlp(f, A, b)  

%%
% Joe can spread his time more evenly, 10 hours each with the running 
% backs, receivers, and linemen; but still the biggest chunk of his time, 
% 20 hours, with the quarterback.

%% 5.
% We let $f$ be the quantity to be set to zero.  It is a function
% of the diode voltage, $V_D$, and the parameters $V_0, R, I_0$, and $V_T$.
% Note that $f$ should be zero for some $x$ between 0 and $V_0$.
f = @(V, V0, R, I0, VT) V - V0 + R*I0*exp(V/VT)  

%% (a)

f1 = @(V) f(V, 1.5, 1000, 10^(-5), .0025)
VD = fzero(f1, [0, 1.5])  
%%
% That's the voltage; the current is therefore 
I = (1.5 - VD)/1000  

%% (b)
f2 = @(V) f(V, 1.5, 1000, 10^(-5)/2, .0025);
fzero(f2, [0, 1.5])  
%%
% Not surprisingly, the voltage goes up slightly.

%% (c)
% This time, we have to be careful; cutting $V_T$ in half means
% squaring $\exp(V/V_T)$. So we have to specify an 
% interval big enough to catch the solution, but small enough to 
% prevent the exponential from blowing up too drastically at the right 
% end-point. This will be the case even more dramatically in part (e)
f2 = @(V) f(V, 1.5, 1000, 10^(-5), .0025/2);
fzero(f2, [0, 0.5])  

%%
% This time the voltage goes down. 

%% (d)
% Next we halve both:

f2 = @(V) f(V, 1.5, 1000, 10^(-5)/2, .0025/2);
fzero(f2, [0, 0.5])  
%%
% The voltage is less than in part (b) but more than in part (c).

%% (e)
X = 10.^(0:-1:-5);
f3 = @(V, x) f(V, 1.5, 1000, 10^(-5)*x, .0025*x);
VD = 0:5; % initialize vector of VD-values
for j = 1:6
    f4 = @(V) f3(V, X(j));
    VD(j) = fzero(f4, [0, X(j)/10]);
end
loglog(10^(-5)*X, VD, 'x-')  
xlabel 'I_0'
ylabel 'V_D'

%%
% The *|loglog|* plot appears linear. That suggests $V_D$ is roughly
% a constant times a power of $I_0$. 

%% 6.

%% (a)
syms x(t) x0 t
dsolve(diff(x) == x - x^2)  

%%

sol = dsolve(diff(x) == x - x^2, x(0)== x0)  

%%
% When simplified, this includes the constant solutions
% 0 and 1; indeed:
bettersol = simplify(sol)  

%%
subs(bettersol, x0, [0 1])  

%%
% Let's also check the initial condition.
subs(bettersol, t, 0)  

%% b)

T = 0:0.1:5;   
 
clf reset, hold on
solcurves = @(t, x0) eval(vectorize(bettersol));
for initval = 0:0.25:2.0
    plot(T, solcurves(T, initval))
end  
axis tight
title( ...
'Solutions of dx/dt = x - x^2, with x(0) = 0, 0.25, ..., 2')
xlabel 't'
ylabel 'x'
hold off  

%%
% The graphical evidence suggests that the solution that starts at
% zero stays there, and all the others tend toward the constant solution
% $1$.
 
%% (c)

clf reset, hold on  
f = @(t,x) [x(1) - x(1)^2 - 0.5*x(1)*x(2); ...
    x(2) - x(2)^2 - 0.5*x(1)*x(2)];
for a = (0:13)/12, for b = (0:13)/12
        [t, xa] = ode45(f, [0 3], [a,b]);
        plot(xa(:, 1), xa(:, 2))
    end, end  
axis([0 13/12 0 13/12])  

%% (d)
% The endpoints on the curves are the start points.  So clearly any
% curve that starts out inside the first quadrant -- that is one that
% corresponds to a situation in which both populations are present at
% the outset -- tends toward a unique point, which from the graph appears
% to be about $(2/3,2/3)$.  In fact, if $x = y = 2/3$, then the right
% side of both equations in (C.4) vanishes, so the derivatives are zero
% and the values of $x(t)$ and $y(t)$ remain constant -- they don't
% depend on $t$.  If only one species is present at the outset, that
% is you start out on one of the axes, then the solution tends toward
% either $(1,0)$ or $(0,1)$, depending on whether $x$ or $y$ is the
% species present.  That is precisely the behavior we saw in part (b).  

%% (e)

clf reset, hold on
f = @(t,x) [x(1) - x(1)^2 - 2*x(1)*x(2); ...
    x(2) - x(2)^2 - 2*x(1)*x(2)];
for a = (0:13)/12, for b = (0:13)/12
        [t, xa] = ode45(f, [0 3], [a,b]);
        plot(xa(:, 1), xa(:, 2))
    end, end
axis([0 13/12 0 13/12])  

%%  
% This time most of the curves seem to be tending toward one of the
% points $(1,0)$ or $(0,1)$ -- in particular, any solution curve that
% starts on one of the axis (corresponding to no initial population for
% the other species) does so.  It seems that whichever species has a
% greater population at the outset will eventually take over all the
% population -- the other will die out.  But there is a delicate balance
% in the middle -- it appears that, if the two populations are about equal
% at the outset, then they tend to the unique population distribution
% at which, if you start there, nothing happens.  That value looks
% like $(1/3,1/3)$.  In fact, that is the value that renders both sides
% of (C.5) zero --- analogous to the role $(2/3,2/3)$ had in part (d).

%% (f)
% It makes sense to refer to the model (C.4) as ``peaceful coexistence,"
% since whatever initial populations you have -- provided both are
% present -- you wind up with equal populations eventually.  ``Doomsday" is
% an appropriate name for model (C.5), since if you start out with unequal
% populations, then the smaller group becomes extinct.  The lower
% coefficient 0.5 means relatively small interaction between the
% species, allowing for coexistence.  The larger coefficient 2 means
% stronger interaction and competition precluding the survival of both.

%% 7.
% Here is a Simulink model for redoing the pendulum application:

open_system pendulum

%%
% With the initial conditions $x(0) = 0$ and $\dot{x}(0) = 10$, the XY
% Graph block shows the following phase portrait of $x$ versus
% $\dot{x}$.
pout = sim('pendulum','SaveTime','on','TimeSaveName','t', ...
       'SaveState','on','StateSaveName','x');
x = pout.find('x'); t = pout.find('t');

%%
% Meanwhile, the Scope block or the following plot
% gives the following graph of $x$ as a
% function of $t$. Note the convergence to $4\pi$ (the
% dotted horizontal line) as $t\to\infty$.

plot(t, x(:,1)), hold on
plot([0 30], [4*pi, 4*pi], ':'), hold off
title 'Plot of x as a function of t'
axis([0 30 0 15])

%% 8.
% Here is a Simulink model for studying the equation of motion of a
% baseball.

open_system baseball

%%
% The way this works is as follows.  The Integrator block
% in the upper left integrates the acceleration (a vector quantity) to
% get the velocity (also a vector).
% This block requires the initial value of the velocity as
% an initial condition; we define it in the ``initial velocity''
% Constant block.  Output from the first Integrator goes into the
% second Integrator, which integrates the velocity to get the position
% (also a vector).  The initial condition for the position, *|[0, 4]|*,
% is stored in the parameters of this second Integrator.  The position
% vector is fed into a Demux block, which splits off the horizontal
% and vertical components of the position.  These are fed into the XY
% Graph block, and also the vertical component is fed into a scope
% block so that we can see the height of the ball as a function of
% time.  The hardest part is the computation of the right-hand side
% of (C.7).
% This is computed by adding the two terms on the right with the Sum
% block near the lower left.  The value of *|[0, -g]|* is stored in the
% ``gravity'' Constant block.  The second term on the right is
% computed in the Product block labeled ``Compute acceleration due to
% drag,'' which multiplies the velocity (a vector) by $-c$ times the
% speed (a scalar).  We compute the speed by taking the dot product of
% the velocity with itself and then taking the square root; then we
% multiply by $-c$ in the Gain block in the middle bottom of the
% model.  The Scope block in the lower right plots the ball's speed as
% a function of time.

%% (a)
% With $c$ set to $0$ (no air resistance) and the initial velocity set
% to $[80, 80]$, the ball follows a familiar parabolic trajectory, as
% seen in the following picture. 
bout = sim('baseball','SaveTime','on','TimeSaveName','t', ...
       'SaveState','on','StateSaveName','x');
x = bout.find('x'); t = bout.find('t');
%%
% Note that the ball travels about $400$ feet before hitting the
% ground, so the trajectory is just about what is required for a home
% run in most ballparks.  We can read off the flight time and final
% speed from the other two scopes or plot them as follows:

plot(t, x(:,2))
axis([0 10 0 150])
title 'Height versus Time'

%%

plot(t, sqrt(x(:,3).^2 + x(:,4).^2))
axis([0 10 50 300])
title 'Speed versus Time'

%%
% Thus the ball stays in the air about $5.0$ seconds and is traveling at
% about $115$ ft/sec when it hits the ground.  (Notice that the
% simulation continues beyond this time; since the model doesn't take
% into account hitting the ground, the height continues to decrease
% and the velocity continues to increase as if the ball had been
% batted off a cliff.)
%
% Now let's see what happens when we factor in air resistance, again
% with the initial velocity set to $[80, 80]$.  First, we take $c =
% 0.0017$.  The trajectory now looks like this:

set_param('baseball/Gain', 'Gain', '-0.0017')
[t, x] = sim('baseball');

%%
% Note the enormous difference air resistance makes; the ball
% travels only about $270$ feet.  We can also investigate the flight time
% and speed with the other two scopes or these plots:

plot(t, x(:,2))
axis([0 10 0 150])
title 'Height versus Time'

%%

plot(t, sqrt(x(:,3).^2 + x(:,4).^2))
axis([0 10 50 300])
title 'Speed versus Time'

%%
% So the ball is about $80$ feet up in the air at its peak, and hits the ground
% in about $4.5$ seconds.

%% (b) 
% Let's now redo exactly the same calculation with $c = 0.0014$
% (corresponding to playing in Denver).  The ball's trajectory is now

set_param('baseball/Gain', 'Gain', '-0.0014')
[t, x] = sim('baseball');

%%
% The ball goes about $285$ feet, or about $15$ feet further than when
% playing at sea level.  This particular ball is probably an easy play,
% but with some hard-hit balls, those extra $15$ feet could mean the
% difference between an out and a home run.  If we look at the height
% vs. time for the Denver calculation, we see: 

plot(t, x(:,2))
axis([0 10 0 150])
title 'Height versus Time'

%%
% so there is a very small increase in the flight time.  Similarly, if
% we look at the speed vs. time for the Denver calculation, we see:

plot(t, sqrt(x(:,3).^2 + x(:,4).^2))
axis([0 10 50 300])
title 'Speed versus Time'

%%
% So the final speed is a bit faster, about $83$ ft/sec.  

%% (c) 
% One would expect that batting averages and
% earned run averages would both be higher in Denver, as
% indeed is the case according to Major League Baseball statistics.

%% 9.
% We construct a Simulink model of the stochastic differential equation.

open_system BS
%%
% First we run the system with the starting parameters, and plot
% the result.
set_param('BS/Constant', 'Value', '0.1');
set_param('BS/Band-Limited White Noise', 'Cov', '0.2');
BSout = sim('BS','SaveTime','on','TimeSaveName','t', ...
       'SaveState','on','StateSaveName','S');
S = BSout.find('S'); t = BSout.find('t');
plot(t, S)
title 'plot of price vs. time for simulation with noise = 0.2'
%%
% Here is what happens if we increase the constant from $0.1$
% to $0.15$, keeping the noise the same.
set_param('BS/Constant', 'Value', '0.15');
BSout = sim('BS','SaveTime','on','TimeSaveName','t', ...
       'SaveState','on','StateSaveName','S');
S = BSout.find('S'); t = BSout.find('t');
plot(t, S)
title 'plot of price vs. time with bigger growth parameter'
%%
% Here is what happens if we increase the noise parameter, going
% back to the original growth parameter.
set_param('BS/Constant', 'Value', '0.1');
set_param('BS/Band-Limited White Noise', 'Cov', '0.3');
BSout = sim('BS','SaveTime','on','TimeSaveName','t', ...
       'SaveState','on','StateSaveName','S');
S = BSout.find('S'); t = BSout.find('t');
plot(t, S)
title 'plot of price vs. time with bigger variance'

%%
% The graphs in this solution look qualitatively similar to
% some of those in the Financial Data application, though
% the time scales are different.

%% 10.
% Here are the contents of our MuPAD notebook.

%% 11.
% As suggested in the hint, we start by defining an array *|Z|* of complex
% numbers whose real parts are chosen from the vector *|xvals|*
% and whose imaginary parts are chosen from the vector *|yvals|*.

xvals = linspace(-2, 2, 400);
yvals = linspace(-1.5, 1.5, 300);
[X, Y] = meshgrid(xvals, yvals);
Z = X + i*Y;

%%
% Next we then apply the given function $100$ times to each number
% in the array.  By that time, most of the sequences that are going
% to diverge to infinity have already become larger than the largest
% floating-point number that MATLAB can store.

for k = 1:100
    Z = Z.^2 - 0.75;
end

%%
% We now use *|isfinite|* to assign a value of $1$ to the values
% in the array that are still finite floating-point numbers and
% $0$ to those that have diverged.  We then use *|imagesc|* to
% color the points with ones differently from those with zeros,
% and select a gray colormap to make the colors white and black.
% Finally, we use *|set|* to reverse the vertical axis; without
% this command, the lower numbers would appear at the top of the
% graph.  This is the default for *|image|* and *|imagesc|*, but
% is not the way mathematical graphs are normally drawn.

clf reset
set(gcf, 'Color', 'White')
imagesc(xvals, yvals, isfinite(Z))
colormap(gray)
set(gca, 'YDir', 'normal')
axis equal tight

%%
% The Julia set consists (approximately) of the boundary between the
% black and white regions.  This boundary is an example of a
% ``fractal,'' and you can draw a variety of other fractals by changing
% the expression $z^2 - 0.75$.  For example, try changing $0.75$
% to $0.75 + 0.1i$ and see what a difference it makes!
