Practice Set C:
Problems 1-3
1.
(a)
radiation = inline(vectorize('10000/(4*pi*((x - x0)^2 + (y - y0)^2 + 1))'), 'x', 'y', 'x0', 'y0') �
���� Inline function:
���� radiation(x,y,x0,y0) = 10000./(4.*pi.*((x - x0).^2 + (y - y0).^2 + 1)) �
x = zeros(1,
5); y = zeros(1, 5);�
for j = 1:5
� x(j) = 50*rand;
� y(j) = 50*rand;
end�
[X, Y] =
meshgrid(0:0.1:50, 0:0.1:50);�
�It is not so clear from the picture where to
hide, although it looks like the Captain has a pretty good chance of surviving
a small number of shots.� But 100 shots
may be enough to find him. Intuition says he ought to stay close to the
boundary.
(b)
Below
is a series of commands that place Picard at the center of the arena, fires the
death ray 100 times, and then determines the health of Picard.� It uses the following function file lifeordeath, which computes the fate of
the Captain after a single shot.
function r = lifeordeath(x1, y1, x0, y0)
%This file computes the number of illumatons.
%that arrive at the point (x1, y1), assuming the
death,
%ray strikes 1 meter above the point (x0, y0).
%If that number exceeds 50,a "1" is
returned in the
%variable "r"; otherwise a "0"
is returned for "r".
dosage = 10000/(4*pi*((x1 - x0)^2 + (y1 - y0)^2 +
1));
if dosage>50
��� r = 1;
else
��� r = 0;
end
Here
is the series of commands to test the Captain's survival possibilities:
for n = 1:100
� x1 = 50*rand;
y1 = 50*rand;
r = lifeordeath(x1, y1, x0, y0);
h = h + r;
end
if h>0
� disp('The Captain is
dead!')
else
� disp('Picard lives!')
end ��
In fact if you run this sequence of commands
multiple times, you will see the Captain die far more often than he lives.
(c)
So
let's do a Monte Carlo simulation to see what his odds are:
for k = 1:100
� h = 0;
for n = 1:100
x1 = 50*rand;
y1 = 50*rand;
r = lifeordeath(x1, y1, x0,
y0);
h = h + r;
end
if h>0
� c = c;
else
� c = c + 1;
end
end
disp(strcat('The chances of Picard surviving are = ', num2str(c/100))) �
The chances of Picard surviving are =0.12 �
We ran this a few times and saw survival chances
ranging from 9% to 16%.
(d)
for k = 1:100
� h = 0;
for n = 1:100
x1 = 50*rand;
y1 = 50*rand;
r = lifeordeath(x1, y1, x0,
y0);
h = h + r;
end
if h>0
� c = c;
else
� c = c + 1;�
end
end
disp(strcat('The chances of Picard surviving are = ',
num2str(c/100))) ��
The chances of Picard surviving are =0.16 �
This time the numbers were between 10 and
18%.� Let's keep moving him toward the
periphery.
(e)
for k = 1:100
� h = 0;
for n = 1:100
x1 = 50*rand;
y1 = 50*rand;
r = lifeordeath(x1, y1, x0,
y0);
h = h + r;
end
if h>0
� c = c;
else
� c = c + 1;�
end
end
disp(strcat('The chances of Picard surviving are = ',
num2str(c/100))) �
The chances of Picard surviving are =0.45 �
The numbers now hover between 36 and 47% upon
multiple runnings of this scenario; so finally, suppose he cowers in the
corner.
x0 = 50; y0 = 50; c = 0;
for k = 1:100
� h = 0;
for n = 1:100
x1 = 50*rand;
y1 = 50*rand;
r = lifeordeath(x1, y1, x0, y0);
h = h + r;
end
if h>0
� c = c;
else
� c = c + 1;�
end
end
disp(strcat('The chances of
Picard surviving are = ', num2str(c/100)))
The chances of Picard surviving are =0.64 �
We
saw numbers between 56 and 64%.�
They
say a brave man dies but a single time, but a coward dies a thousand
deaths.� But the person who said that
probably never encountered a Cardassian.�
Long live Picard!
2.
(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‰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 clearly 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‰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 + � + (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[((1 + J)n+1
- 1)/((1 + J) �� 1) � 1] = S[((1 + J)n+1
� 1)/J � 1].
(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 = .1/12.
months =
solve('(1 + (0.05)/12)^n = 2')�
years = months/12�
166.70165674865177999568182581405
years =
13.891804729054314999640152151171 �
months =
solve('(1 + (0.1)/12)^n = 2')�
years = months/12�
83.523755900375880189555262964714
years =
6.9603129916979900157962719137262 �
If
you double the interest rate, you halve the time to achieve the goal.
(e)
solve('1000000 = S*((((1 + (0.08/12))^(35*12 + 1) - 1)/(0.08/12)) - 1)') �
433.05508895308253849797798306477 �
You
need to deposit $433.06 every month.
(f)
solve('1000000 = 300*((((1 + (0.08/12))^(n + 1) - 1)/(0.08/12)) - 1)') �
472.38046393034711345013989217261 �
ans/12 �
39.365038660862259454178324347717 �
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
format bank; option1 = 65000*(1 + (.05/12))^(12*20) �
���� 176321.62 �
The
stash grows to about $176,322.� The
second option yields
������� 333.33 �
option2 = S*((1/(.05/12))*(((1 + (.05/12))^241) - 1) - 1) �
���� 137582.10 �
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]; �
for k = 0:4
� T = 50000;
for j = 1:5
� T = T*(1 + rates(k + j));
end
disp([k + 1,T])
end �
��������� 2.00����� 72794.74
��������� 3.00���� �72794.74
��������� 4.00����� 72794.74
��������� 5.00����� 72794.74 �
The
results are all the same, you wind up with $72,795 regardless of where you
enter in the cycle, because the product �1�j�5
(1 + rates(j)) is independent of the order in which
you place the factors.� If you put the
50K in a bank account paying 8%, you get
����� 73466.40 �
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 + r1); then at the end of the
second year (S(1 + r1) + S)(1 + r2),
where we have written rj
for rates(j).� So at the end of 5 years,
the amount in the account will be the product of S and the number
�j�1
(1 + rj) + �j�2
(1 + rj) + �j�3
(1 + rj) + �j�4
(1 + rj) + (1 + r5).
If
you enter at a different year in the business cycle the terms get cycled
appropriately.� So now we can compute
� 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�
��� 2.0000��� 6.4000
��� 3.0000��� 6.8358
��� 4.0000��� 6.1885
��� 5.0000��� 6.0192 �
Multiplying
each of these by $10,000 gives the portfolio amounts for the five
scenarios.�� Note all are less than what
one obtains by investing the original $50,000 all at once�not surprising.� 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 �
��� 6.3359 �
which
is better than some investment models and worse than others.
3.
(a)
First
we define an expression that computes whether Tony gets a hit or not during a
single at bat, based on a random number chosen between 0 and 1.� If the random number is less than or equal
to 0.339, Tony is credited with a hit, whereas if the number exceeds 0.339, he
is retired by the opposition.�
Here
is an M-file, called atbat.m, which computes the outcome of a single at bat:
%This file simulates a single at
bat.
%The variable r contains a
"1" if Tony gets a hit,
%that is, if rand <= 0.339; and
it contains a "0"
%if Tony fails to hit safely, that
is, if rand > 0.339.
s = rand;
if s <= 0.339
��� r = 1;
else
��� r = 0;
end
We
can simulate a year in Tony's career by evaluating the script M-file atbat 500 times.� The following program does exactly that,
then it computes his average by adding up the number of hits and dividing by
the number of at bats, that is 500.� We build
in a variable that allows for a varying number of at bats in a year, although
we shall only use 500.
function y = yearbattingaverage(n)
%This function file computes
Tony's batting average for
%a single year, by simulating n at
bats, adding up the
%number of hits, and then dividing
by n.
X = zeros(1, n);
for i = 1:n
��� atbat;
��� X(i) = r;
end
y = sum(X)/n;�
yearbattingaverage(500) �
ans =
��� 0.3200 �
(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.
function y = career(n,k)
%This function file computes the
batting average for each
%year in a k-year career, asuming
n at bats in each year.� %Then it lists
the maximum, minimum and lifetime average.
Y = zeros(1, k);
for j = 1:k
��� Y(j) =
yearbattingaverage(n);
end
disp(strcat('Best avg: ',
num2str(max(Y))))
disp(strcat('Worst average: ',
num2str(min(Y))))
disp(strcat('Lifetime avg: ',
num2str(sum(Y)/k)))
Worst average:0.29
Lifetime avg:0.341 �
(c)
Now
we run the simulation four more times:
Worst average:0.292
Lifetime avg:0.3342 �
Worst average:0.29
Lifetime avg:0.345 �
career(500, 20)�
Best avg:0.378
Worst average:0.312
Lifetime avg:0.3428 �
Worst average:0.324
Lifetime avg:0.3473 �
(d)
The
average for the five different 20-year careers is:
(.344 + .340 + .338 + .343 + .331)/5 �
ans =
��� 0.3392 �
How
about that!
If
we ran the simulation 100 times and took the average it would likely be very
close to .339.