2.9 exercises.ipynb

Kernel: SageMath 7.6

Remember the M&Ms project? This project involved the difference equation $y_{n+1}=f(n,y_n)=\dfrac{1}{2}y_n+B$ where the value of $B$ was to be chosen, within limits.

Let's answer the question of this project as an introduction to section 2.9.

a) State your assumptions about the physical activity.
The usual simplifying assumptions are that each M&M has exactly on "m" on one face or the other (not both) and that, when the M&M's are gently shaken onto the plate, each "m" is equally likely to appear face up as it is to appear face down. I will also assume that most pigs can fly.

b) Offer up a description of what should happen. Care to make a prediction?
About half the M&M's will be removed from the population after each iteration. This will, I predict, lead to all the M&M's being removed after about $log_2(N) iterations.

c) In Table 1 record what happened and compare with what you thought would happen.
I will simulate this using the tabulate HTML library, like so:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML, display
import tabulate
table = [["Iteration","Number of M&M's at start of iteration"],
         [0,50],
         [1,29],
         [2,15],
         [3,7],
         [4,5],
         [5,2],
         [6,0]]
display(HTML(tabulate.tabulate(table, tablefmt='html')))

Iteration	Number of M&M's at start of iteration
0	50
1	29
2	15
3	7
4	5
5	2
6	0

np.log2(50)

5.6438561897747244

d) Compare your description/prediction with what actually happened.
$log_2(50)\approx 5.6, so the behavior I experienced was very close to what I expected.

e) Indicate which of your assumptions were reasonable and played a role in the experiment.
All were reasonable...except the bit about pigs.

f) Based on the observations and your assumptions, produce a reasonable formula for $a(n)$, i.e. offer up a discrete function $a(n)$, in the one variable, $n$, for $n = 0, 1, 2, \ldots$. Discuss the reasonableness of your function model for $a(n)$.
Based on our assumption, the best estimate of how many M&Ms survive is 1/2, so we'd expect $a(n)=a(n-1)/2$

g) How will you measure your "success" as a modeler in this situation? Check with others in the class and see if you are on target with your model. To measure success I would compare the absolute errors with the impirical data:

table = [["Iteration","M&M's","a(n)","a(n)-M&Ms"],
         [0,50,50,0],
         [1,29,25,-4],
         [2,15,12.5,-2.5],
         [3,7,6.25,-0.75],
         [4,5,3.125,-1.875],
         [5,2,1.5625,-0.4375],
         [6,0,0.78125,0.78125]]
display(HTML(tabulate.tabulate(table, tablefmt='html')))

Iteration	M&M's	a(n)	a(n)-M&Ms
0	50	50	0
1	29	25	-4
2	15	12.5000000000000	-2.50000000000000
3	7	6.25000000000000	-0.750000000000000
4	5	3.12500000000000	-1.87500000000000
5	2	1.56250000000000	-0.437500000000000
6	0	0.781250000000000	0.781250000000000

h) Check how good your model is at predicting your experiment and defend your rationale.
The errors are pretty small...so it's a pretty good model.

i) State your assumptions about the physical activity.
The assumptions are the same, except I'll be adding 14 M&M's at each iteration.

j) Offer up a description of what will happen.
The M&M population will decrease by half and then increase by 14, a process modeled by
$a(n+1)=a(n)/2+14$

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import HTML, display
import tabulate
table = [["Iteration","Number of M&M's at start of iteration"],
         [0,50],
         [1,37],
         [2,34],
         [3,22],
         [4,24],
         [5,24],
         [6,30],
         [7,26],
         [8,28],
         [9,27],
         [10,26]]
display(HTML(tabulate.tabulate(table, tablefmt='html')))

Iteration	Number of M&M's at start of iteration
0	50
1	37
2	34
3	22
4	24
5	24
6	30
7	26
8	28
9	27
10	26

k) In the above table I recorded what happened.
Huh. It seems to be reaching a steady state...kind of, hovering around 28.

l) To verify that, consider the
model in a steady state: $a(n+1)=a(n)=a(n)/2+14$ and solve for $a(n)$ to get $a(n)=28$. Note that if $a(n+1)<28$
then $\Delta a=a(n+1)-a(n)>0$ and if $a(n+1)>28$ then $\Delta a=a(n+1)-a(n)<0$ so 28 is a stable equilibrium.

n) If we assume the population of M&Ms is a continuous (rather than discrete) variable, then we can solve $\dfrac{da}{dt}=-\dfrac{a}{2}+B$
by using the integrating factor $\mu(t)=e^{t/2}$ to get $\dfrac{d}{dt}(e^{t/2}a)=Be^{t/2}\Leftrightarrow e^{t/2}a=2Be^{t/2}+c$
So $a(t)=2B+c\cdot e{-t/2}\rightarrow 2B$ as $t\rightarrow\infty$. That confirms our experimental results.

#Non-linear Difference Equations.

Nonlinear difference equations are much more complicated and have
much more varied solutions than linear equations. We will restrict our attention to a single
equation, the logistic difference equation
$y_{n+1}=\rho y_n\left(1-\dfrac{y_n}{k}\right)$ (19)
which is analogous to the logistic differential equation
$\dfrac{dy}{dt}=ry\left(1-\dfrac{y}{K}\right)$ (20)
that was discussed in Section 2.5. Note that if the derivative $dy/dt$
in equation (20) is replaced by the difference quotient $\dfrac{y_{n+1}−y_n}{h}$,
then equation (20) becomes $\dfrac{y_{n+1}−y_n}{h}=ry_n\left(1-\dfrac{y_n}{K}\right)\Leftrightarrow y_{n+1}=y_n+hry_n\left(1-\dfrac{y_n}{K}\right)$
$=y_n+y_n\left(hr-\dfrac{hry_n}{K}\right)=y_n\left(1+hr-\dfrac{hry_n}{K}\right)=y_n\left[(1+hr)\left(1+\dfrac{hry_n}{K(1+hr)}\right)\right]$
$=(1+hr)y_n\left(1+\dfrac{y_n}{K(1+hr)/hr}\right)$ reduces to equation (19) with $\rho=1+hr$ and $k=(1+hr)K/hr$. Phew!

0.0075*100000/(1-1.0075^(-240))

899.725955850170