Module 3 - Dynamic models


This is a hand-in written by Group 60 (FIRE) for the weekly module #3 in the mathematical modelling course (DAT026) at Chalmers University of Technology.

Group 60 consists of:

  • Mazdak Farrokhzad

    • 901011-0279


    • Program: IT

    • Time spent: x hours, y minutes

  • Niclas Alexandersson

    • 920203-0111


    • Program: IT

    • Time spent: 18 hours, 48 minutes

We hereby declare that we have both actively participated in solving every exercise, and all solutions are entirely our own work.

Disclaimer: Please do note that the numbering of the sections in the TOC and elsewhere are not intended to follow the numbering of the problem. The problem referred to is always exlicitly stated in the title itself (e.g: Problem 1...)

Problem 1 - An epic tale of Krill and Whale

We are given two differential equations modelling the population dynamics of whales \((w)\) and krill (\(k\)), with the constants \(a, b, m, n\):

\[\begin{aligned} k' &= (a-bw) k\\ w' &= (-m+nk) w\end{aligned}\]

a) Interpretation of differential equations and constants

The system of non-linear differential equations create a system of co-dependently bounded exponential growth.

In this system, the constants used should have the follow meanings:

  • \(a =\) how good krill are at procreating as a factor how the current krill population size.

  • \(b =\) how much harder survival for krill becomes in relation to the current whale population size.

  • \(n =\) an inverse factor of how many krill it takes to feed a whale, and affects how the whale population grows with the current krill population size - the smaller the constant, the more krill per whale is needed to grow the population of whales. In this system of differential equations, the term \(nkw\) does not directly entail that if the krill population decreases, so will the whale population - but when \(m > nk\), the whale population will decrease.

  • \(m =\) death rate as a factor of current whale population size, may be due to “natural” causes such as dying due to age or whale hunting done by humans.

It is important to note that the equations are only defined for \(k(t) \geq 0, w(t) \geq 0, \text{where } t \geq 0\). A negative population size has no meaning.

b) Points of equilibrium

For an equilibrium to occur, the rate of change for both the krill and the whale population needs to be \(0\). I.e. we have \(k'=0\) and \(w'=0\). Looking at the equations, we can quickly see that for \(k'\) to equal \(0\), either \(a - bw\) or \(k\) must be \(0\). For \(w'\), \(-m+nk\) or \(w\) must be \(0\).

In an equilibrium where \(k=0\), we either have that \(m=0\) or \(w=0\). Since \(m\) is a constant which will likely not be \(0\) in practice, the only realistic case here is that \(w=0\). The same is true for \(a\) when \(w=0\). This would correspond to the case where both populations have gone extinct and cannot reproduce any more. When both species are extinct the equilibrium is fully stable.

More interesting perhaps is the case where \(a - bw = 0 \iff w = \frac{a}{b}\) and \(-m + nk = 0 \iff k=\frac{m}{n}\). This would mean that the populations have reached a state where the krill population is at a level exactly as large as needed to sustain a whale population of a size exactly as large as needed to keep the population level of krill constant.

c) Manual simulation with Eulers method

Without knowing the functions \(k(t), w(t)\) we are going to find values for them using numerical integration. The basic idea behind numerically finding the values is, for some function \(y(t)\):

\[y_{n+1} = y_n + \Delta t\frac{\Delta y}{\Delta t}\]

For differential equations, the Euler method is a usable method:

\[y_{n+1} = y_n + \Delta t y'(t_n, y_n)\]

The smaller the step size \(\Delta t\) is, the more accurate the value of \(y_{n+1}\) becomes.

Replacing the constants in the diff equations with the values: \(a = 0.2, b = 10^{-4}, m = 0.5, n = 10^{-6}\), we get:

\[\begin{aligned} k' &= (0.2 - 10^{-4} w) k\\ w' &= (-0.5 + 10^{-6} k) w\end{aligned}\]

Given the starting points, and the time step:

\[\begin{aligned} k(0) &= 7 \cdot 10^5\\ w(0) &= 3 \cdot 10^3\\ \Delta t &= 0.3\end{aligned}\]

we manually simulate \(k(t), w(t)\) for \(t_1 = 0.3, t_2 = 0.6\).

\(n\) \(t_n\) \(\Delta t\) \(k_n\) \(w_n\) \(k'(t_n, k_n)\) \(w'(t_n, w_n)\) \(k_{n+1}\) \(w_{n+1}\)
\(0\) \(0\) \(0.3\) \(700000\) \(3000\) \(-70000\) \(600\) \(679000\) \(3180\)
\(1\) \(0.3\) \(0.3\) \(679000\) \(3180\) \(-80122\) \(569.22\) \(654963\) \(3350.77\)
\(2\) \(0.6\) \(0.3\) \(654963\) \(3350.77\)        

For ease of computing, the following functions were defined and used in mathematica: kp[k_, w_] := (0.2 - 10^(-4) *w)* k wp[k_, w_] := (-0.5 + 10^(-6) * k) * w prims[k_, w_] := {kp[k, w], wp[k, w]} y[k_, w_] := {k + 0.3 kp[k, w], w + 0.3 wp[k, w]} all[k_, w_] := {prims[k, w], y[k, w]}

If we plot the data table, the graphs like this:

[ xlabel=\(t\), ylabel=\(k\), height=9cm, width=12cm, grid=major]

+[smooth] table[x=t, y=k] ;


[ xlabel=\(t\), ylabel=\(w\), legend style=at=(0.9,0.1),anchor=south east, height=9cm, width=12cm, grid=major]

+[smooth] table[x=t, y=w] ;