Scientific Computing

Assignment 1

Maastricht Science Program 
Coordinator: Frans Oliehoek <frans.oliehoek@maastrichtuniversity.nl>

Student

Name: <ENTER YOUR NAME>
Email: <ENTER YOUR EMAIL>

Working with Mathematica Notebooks

Currently, you are looking at a mathematica ‘notebook’, which is the document type used by mathematica. It allows you to combine text and computation in one document. This is facilitated by using ‘cells’ which are indicated by the ‘]’ brackets in the right margin. 
For instance, the text that you are currently reading is a ‘text’ cell...

...and this is another...

etc.
You can verify this by positioning the cursor on this text, or clicking on the left (inner) most bracket, and clicking on the “Format->Style” dropdown menu”. This will also allow you to change the cell types and show you the shortcut keys to do this.
Try inserting a new cell here: position the cursur on the last line of this paragraph, a line should appear, then press <ALT>+7 and start typing.

If you find these cells confusing, you can have a quick look at this video:
http://www.wolfram.com/broadcast/screencasts/handsonstart/
(but please use a headphone!)

Interleaving Code and Text

Of course, having just text does not help us solve interesting problems, the interesting part of Mathematica is that is also allows you to perform (very sophisticated) computations!
For instance, put the cursor on the cell below, and ‘execute’ or ‘evaluate’ it by pressing <SHIFT>+<ENTER>.

2+3

These executable cells are called ‘input’ cells. As you may have noticed, Mathematica automatically creates an output cell when you evaluate an input cell. You can also use the numpad <ENTER> key

7*8

As you can see this is quite versatile. It is also possible to have multiple commands in one cell as this:

4-7-2
7/8
(1+1)/3
4^2

To create an input cell your self, just position the cursur in a non-cell area such that the ‘insertion line’ appears and start typing.

Assignment:

Compute the number of seconds in a year. To do this, you have to insert a new ‘input cell’ right below this cell.

You will be working in this notebook, and I want to see the results. Therefore, it is good to regularly save your work (as usual the save option is under the file menu).

Assignment:

Save this notebook.

Using Functions

One of the nice things in mathematica is that it comes with many build in functions. For instance try running the next cells and understanding what they do.  
If you are uncertain, try consulting the mathematica help function: right click a function name and select “get help”,  or open the documentation center and search for the term you want help about.

Pi

N[Pi]

N[Pi, 20]

%

2*%

Plot[Sin[x^2],{x,0,2Pi}]

Integrate[1/(x^3+1),{x,0,1}]
TraditionalForm[%]

myList={4,1,2,3}
Sort[myList]

Note that this does not change the list:

myList

myList[[3]]

Permutations[myList]

---------------------------------------------------------
IMPORTANT 

note that: 
-Function names start with capitals
-Function arguments go in square brackets: [ ]
-Curly braces { } are used for lists
-Round parens () are used for grouping toghether computations.

---------------------------------------------------------

Assignment:

Call me  over and explain the above examples.

Don’t Panic!

Sometimes, evaluation does not seem to end. To abort evaluation, press: ALT+.
(see also the options in the ‘Evaluation’ pull down menu)

Try to run this command, and abort it:

Integrate[1/(x^1.3+1),{x,0,1}]

(Integrate is a command for integration, one of the powerful functions of Mathematica that we will revisit in a later lab session. If you are already curious, test the following commands:

Integrate[1/(x^2),x]
Integrate[1/(x^2),{x,0.1,1}]

)

Using Variables

Variables are useful to better understand your computations. For instance, if we want to figure out the number of minutes in 15 weeks, this:

HoursPD=24;
MinutesPH=60;
numWeeks=15;
numDays=numWeeks*7;
numMinutes =numDays*HoursPD*MinutesPH

is clearer than

24*60*15*7

Of course, there are some rules to using variables:

only use letters

capitalization matters!

as a rule, you should start variable names with a lowercase letter. (names starting with a capital are ‘reserved’ for mathematica)

For more information, see “tutorial/DefiningVariables” in help.

When using variables, it may happen that you assign a value to a variable, which you later want cleared. For this, there is a special ‘Clear’ command. It is best demonstrated by an example.

temp=6

temp

Clear[temp]

temp

(Did you notice how the color of the variable ‘temp’ changes to black if it is defined, and blue when it is not)

Assignment

Compute the number of inches in one light year. Use variables to make the computation understandable.

Trajectory of a Ball: Defining Functions & Plotting

In this assignment, we explore defining functions, and plotting.

Using functions is convenient, but of course Mathematica does not define all the functions you might want to use. Therefore, it is convenient to define your own functions. For instance if you want to plot them.

Here we will define a simple linear function:

y[x_]:=5*x+3

y[4]

There are a couple of things to note about the definition of the function that we called ‘y’

the underscore ‘_’ after ‘x’ indicates that ‘x’ is an argument.

we use ‘:=’, called the SetDelayed operator,  instead of ‘=’, the normal assignment operator. The difference is that a function defined with ‘:=’ will only try and figure out what value the right hand side (rhs) has when it is called, while ‘=’ tries to figure out this value when the function is defined.

The difference between immediate and delayed assignment are actually a bit subtle (see “tutorial/ImmediateAndDelayedDefinitions” in help, if you are interested. Usually, however, it is better to use the delayed version, as illustrated in the following example.

Clear[x]

immFunction[x_]=2*x

delFunction[x_]:=2*x

Usually this works fine:

immFunction[10]
delFunction[10]

Howevever, look what happens if we happened to have defined ‘x’ before:

x=4
immFunction[x_]=2*x

immFunction[10]

In contrast,  the following still works:

delFunction[x_]:=2*x
delFunction[10]

Clear[x]

Now, returning to our function ‘y’, we can use this function for plotting, using the following syntax:

Plot[y[x],{x,-1,8}]

Of course, we often want to plot legends and titles, these can be added as options:

Plot[y[x],{x,-1,8},AxesLabel->{"x-as", "y-as"}, PlotLabel->"Plot of function y" ]

Assignment:

Try right-clicking the figure and saving it as ‘eps’ (encapsulated postscript - this is great for use in latex!) or ‘pdf’ (great for pdflatex!)

Now, we are in the position to consider the problem of the trajectory of a ball that is thrown. In particular, we are interested in the height of the ball versus the horizontal distance traveled.

We will use the following symbols
    g: the gravitational acceleration—usually taken to be 9.81 m/s2 near the Earth’s surface
    theta: the angle at which the projectile is launched
    v0: the velocity at which the projectile is launched
    y0: the initial height of the projectile
    d: the total horizontal distance traveled by the projectile
(See also: http://en.wikipedia.org/wiki/Trajectory_of_a_projectile)

Clear[x,y]
g=9.81;
theta=.2*Pi;
v0=10;
y0=1.8;

The height of the ball at a particular distance is given by the following formula:

y(x)=y0+x tan(theta)-((g (x^2))/(2 (((v0 cos(theta)))^2)))

Assignment

Define a function ‘y’ and plot it over the interval x∈[0,10]. Add labels to the axis, as well as a meaningful title.

Assignment

Define a function ‘yNoG’ for the trajectory of the same ball if there is no gravity. Plot it together with ‘y’ in one plot.
(Hint: `plot[ {y1[x], y2[x]}, {x,0,10}]’ plots two functions ‘y1’ and ‘y2’ in the same plot.)

Growth of Bacteria

(This assignment is heavily inspired by Allman&Rhodes, 2004, Mathematical Models in Biology.)

In this assignment we investigate the so-called Maltusian model of growth of bacteria. This (simplistic) model assumes that we can consider the size of the population at discrete time steps, and that at each time step, the growth is given by:

(Subscript[P, t])= (Subscript[P, t-1])+ΔP

ΔP=b*P-d*P=(b-d)*P

with b  the birth rate, d  the death rate, and P the population size. As you might know, such an equation, where the next value depends on the previous value is known as a difference equation.

Since we are interested in the size of the population over time, we can express the above equations as

(Subscript[P, t])= (Subscript[P, t-1])+ΔP=(Subscript[P, t-1])+(b-d)*(Subscript[P, t-1])=r*(Subscript[P, t-1])

where we implicitly defined r=(1+b-d). 
Subsequently, we can express the population at time t as

(Subscript[P, t])=r*(Subscript[P, t-1])=r*(r*(Subscript[P, t-2]))= ... = (r^t)*(Subscript[P, 0])

Assignment

Define a function ‘pop[t]’ to express the population of bacteria, assuming a Maltusian model, using:

b=1.47;
d=1.42;
P0=7;

Subsequently, plot the population.

Let Mathematica do the work.

Now, we consider how we can  let mathematica do the hard work for us. That is, rather than doing all the math above, ourselves, we can let Mathematica deal with this.

Clear[pop,n];
deltaP[n_]:=(b-d)*pop[n];

RSolve[{pop[n+1]==pop[n]+deltaP[n],pop[0]==17}, pop[n],n]

At this point, mathematica has found an expression of the function P[n]. However, as you can see, the output is in a ‘strange’ form. It is a list of lists.  Let us save the first (and only!) element of this list in a new variable ‘popFunc’

popFunc=First[%]

This ‘popFunc’ variable now contains a ‘rule’, which can be understood as a substitution command. To work with the result returned by RSolve, we need to use a ‘ReplaceAll’ command, indicated by ‘./’
For instance, we can now replace the expression ‘pop[n]’  with 17.` ((0.9345794392523364`)^(-1.` n)) using:

pop[n]/.popFunc

This, in turn, we can use to create a table :

t=Table[pop[n]/.popFunc,{n,0,15}]

The points in such a table can be plotted using the ‘ListPlot’ command

ListPlot[t]

Assignment

Create a list of the first 50 points of ‘pop[n]’ and use ListPlot to plot them. Compare the result to the plot you made earlier.

More Complex models.

However, when we make the model a bit more complicated, even Mathematica cannot find closed form solutions... (they don’t exist!) We will illustrate this by now using a so-called finite capacity model of population growth. In this model, there are only enough resources to maintain K individual (called the ‘carrying capacity’). The formula for this model is given by

ΔP= (Subscript[P, t])(b-d)(1-((Subscript[P, t])/K))

Note that this leads to a growth of nearly r when (Subscript[P, t]) is small, and a growth of 0 when (Subscript[P, t])=K. 
This leads to

(Subscript[P, t+1])=(Subscript[P, t])+(Subscript[P, t]) (b-d)(1-((Subscript[P, t])/K))=(Subscript[P, t])(1+(b-d)(1-((Subscript[P, t])/K)))

Putting this in Mathematica as before:

Clear[K,n];
K=400;
RSolve[{popk[n+1]==popk[n]*(1+(b-d)*(1-((popk[n])/K))),popk[0]==P0}, popk[n],n]

You will notice that, in contrast to before, mathematica did not produce a nice formula. This means that it could not find a closed form solution (which typically means that no closed-form solution exists).
Not all hope is lost, however. We can still evaluate the problem numerically, by creating a table:

Clear[popk,n]

t=RecurrenceTable[{popk[n+1]==popk[n]*(1+(b-d)*(1-((popk[n])/K))),popk[0]==P0}, popk,{n,0,200}]

Remark:

You may have noticed that in the above examples, we are using the function ‘popk’ which has not been defined at all. In fact, ‘popk’ (and also ‘pop’ before) are not functions of the type we defined before. Instead, these are parts of the specification of a (system of) ‘difference equation(s)’, also called recurrence equations. That is, in the last example
{popk[n+1]==popk[n]*(1+(b-d)*(1-((popk[n])/K))),popk[0]==P0}
is a set of two equations. Mathematica can interpret these equations without us having to define actual functions  (as we did in the trajectory of a ball assignment.)

Assignment

Plot the growth according to the finite capacity growth model.

Example: Gene Frequency in Populations

(Again, this is heavily inspired by the treatment by Allman&Rhodes, 2004, Mathematical Models in Biology.)

We end this session with a slightly more complex example: a model gene frequencies. In particular, we will model how for a particular gene (e.g., encoding eye color), the frequency of the alleles (variants of that gene, e.g. ‘blue’, ‘green’, ‘brown’) changes over time.

We will consider a hypothetical gene, which has two alleles P, Q. We will write p to denote the frequency of P, and q for the frequency of Q.
Since every individual has two copies of each gene, an individual can have genotype PP, PQ, or QQ. From the frequency of genotypes, we will determine the allele frequencies. Then we will examine how they change (or not!) under random mating. Finally, we consider the scenario where certain genotypes lead to a reproductive advantage.

Determining Allele Frequencies

Suppose that at time 0, we have that the frequency of genotypes is as follows:

PP=.09;
PQ=.42;
QQ=.49;

from this we can derive that the total number of P alleles is

totalPs=2*PP*popsize+PQ*popsize

(can you explain the ‘2*’ part?!)
and similarly,

totalQs=2*QQ*popsize+PQ*popsize

From these, we can now derive the frequencies of the alleles:

p0=totalPs/(totalPs+totalQs)
q0=totalQs/(totalPs+totalQs)

(notice that ‘popsize’ nicely drops out, without us ever having defined it!)
Finally, we can check that our results make sense:

{p0^2,2*p0*q0,q0^2}

Random Mating and Weinberg Equilibrium

Suppose that we consider a ‘random mating’ model of reproduction: every individual in the population is just as likely to produce offspring with any other individual. Also, for convenience, we will ignore sex completely.
Under such a model, the probability that a newborn receives a P for its first chromosome is

singleP=((PP+.5PQ)/(PP+PQ+QQ))

singleQ=((QQ+.5PQ)/(PP+PQ+QQ))

PPnew=singleP^2
PQnew=2*singleP*singleQ
QQnew=singleQ^2

That is, the frequency of genotypes seems unaltered in the new generation! This would impy that the state of the population is ‘in equilibrium’ (in balance) according to the model of random mating. This is indeed the case, and this type of equilibrium is called Weinberg Equilibrium. We will further investigate why this happens next.

Let us start by rewriting the derivation above in a somewhat more general way:

Clear[p,q]
Clear[PP,PQ,QQ]

Given frequencies p,q, we can derive the frequencies of genotypes

PP=p^2
PQ=2*p*q
QQ=q^2

In turn, these give us allele frequencies for the next generation:

pnew=PP+(1/2)PQ
qnew=QQ+(1/2)PQ

In other words, we can express the new allele frequencies as a function of the previous ones:

p(n+1)==((p(n) q(n)+((p(n))^2))/( p(n) q(n)+((p(n))^2)+((q(n))^2)))

q(n+1)==((p(n) q(n)+((q(n))^2))/( p(n) q(n)+((p(n))^2)+((q(n))^2)))

Note that these are a system of (coupled) difference equations! As such, we might hope that RSolve can again shed some light on what will be the behavior of these sequences over time. 
We do this via:

ClearAll[q,p]
RSolve[{
p[n+1]==(p[n]^2+p[n]q[n])/(p[n]^2+2*p[n]q[n]+q[n]^2),
q[n+1]==(q[n]^2+p[n]q[n])/(p[n]^2+2*p[n]q[n]+q[n]^2),
p[0]==0.3,
q[0]==0.7
}
 , {p[n],q[n]},n]

This confirms that the system is at a steady state, a Hardy-Weinberg Equilibrium, of p=0.3, q=0.7.

Incorporating Relative Fittness

One aspect of the theory of evolution is that certain genotypes may lead to an advantage, while others may cause a disadvantage; leading to differences in the amount of offspring produced. 

Here, we examine how such considerations may influence our model, by introducing ‘relative fittness’ factors for the different genotypes. In particular, we will assume that PP leads to the highest reproduction rate, PQ has a bit of a disadvantage, while QQ has a more severe disadvantage. We model this with the following relative fitness levels:

Fpp=1;
Fpq=.98;
Fqq=.92;

We can now incorporate these fitness factors in the reproduction equations above, and again use RSolve:

RSolve[{
p[n+1]==(Fpp*p[n]^2+Fpq*p[n]q[n])/(Fpp*p[n]^2+Fpq*2*p[n]q[n]+Fqq*q[n]^2),
q[n+1]==(Fqq*q[n]^2+Fpq*p[n]q[n])/(Fpp*p[n]^2+Fpq*2*p[n]q[n]+Fqq*q[n]^2),
p[0]==0.3,
q[0]==0.7
}
 , {p[n],q[n]},n]

We notice that in this more complex problem, mathematica is not able to find a closed-form solution. Still, as before, we can still tabulate the values:

t=RecurrenceTable[
{
p[n+1]==(Fpp*p[n]^2+Fpq*p[n]q[n])/(Fpp*p[n]^2+Fpq*2*p[n]q[n]+Fqq*q[n]^2),
q[n+1]==(Fqq*q[n]^2+Fpq*p[n]q[n])/(Fpp*p[n]^2+Fpq*2*p[n]q[n]+Fqq*q[n]^2),
p[0]==0.3,
q[0]==0.7
},{p,q},{n,1,100}]

Note that this output is a list of {p,q} pairs. We can not directly plot these, but by taking the transpose:

tt=Transpose[t]

we get a list that contains 2 lists: one of p’s and one of q’s. The first of these (the list of p’s) we can extract using tt[[1]], and the second using tt[[2]]. 
(the double brackets denote an operation that extracts an element from a list)

These can now be put in listplot:

Needs["PlotLegends`"]
ListPlot[{tt[[1]],tt[[2]]},PlotLegend->{p,q},LegendPosition->{1.1,-0.4}]

As expected, we see that since P leads to an advantage, its frequency increases over time.

Hand-in Assignment: Predator-Prey Dynamics

NOTE: this is a hand-in assignment, which means that you will have to hand in a mini-report describing the assignment. Instructions for handing in etc. are listed on the course website. For Mathematica assignments, like the current assignment, it will be the easiest to just submit your entire notebook. Make sure however, that you write text to actually make it like a small report.

In this assignment, you are asked to use Mathematica to plot the development of the populations of two linked species. The first, the pray species, will be denoted by P, and in absence of any predators, its population follows a finite capacity growth as above:

ΔP= (Subscript[P, t])r(1-((Subscript[P, t])/K)).

However, when the second, predator, species (denoted by Q) are present, these will eat a fraction of the prey population. The number of pray killed depends on the probability that predators find prays, which we can model with a product ‘interaction’ term: s *(Subscript[P, t])*(Subscript[Q, t ]), with s  a positive constant. As such the pray population evolves as

ΔP= (Subscript[P, t])r(1-((Subscript[P, t])/K))-s (Subscript[P, t]) (Subscript[Q, t ]).

The predator population evolves according to the following model:

ΔQ= (1-u)(Subscript[Q, t])+v (Subscript[P, t]) (Subscript[Q, t ]).

Here, u denotes the fraction of predators that will die at every iteration when no pray are present. When pray are present, however, the predators can feed and reproduce, as modeled by the interaction term v (Subscript[P, t]) (Subscript[Q, t ]). Again v  is a positive constant.

Assignment

Plot the evolution of the populations of predator and prey species for 50 iterations. Make a nice looking plot!
Use the following constants: K=1, r=1.3, s=.5, u=.7, v=1.6, (Subscript[P, 0])=1.1, (Subscript[Q, 0])=0.4.

Everything required for this assignment is treated in earlier examples.
Hints:

First express the population sizes at the next stage in terms of the previous-stage population sizes.

Use some of the recurrence solving functionality!