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
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
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.