Runge Kutta Equations 4. Order¶

an adaptive Simulation of the Lotka-Volterra Model¶

created by Sophie Schmeißner, summer term 2020

Usage¶

You can scale the Reproduction- and Mortality-Koefficient of Predator & Prey. The first Plot shows the time development, the second one the Phase Diagram.

Left Side: The orange Graph shows the Population of Predators , the blue one shows the Population of Preys. You can see the typical Lotka Volterra course: The more Prey, the more Predators can be saturated and start to multiply. The more Predators, the less Preys until the Preys start to die out, so the Predators are not longer saturated. The less Predators, the more Preys will survive and start to multiply. These dynamics continue in a cycle of growth and decline

Right Side: Population of Preys and Predators are plotted against each other. The green dot shows the Starting point, the blue dot shows the Fixpoint. Because of the Periodicity you can't see any start or end, it oscillates with the same values for $s=0$ in every time

Constant of motion¶

Here you can see the known phase diagram with different initial Values. Every initial value leads to periodic solutions, but for identical Koefficients they will have the same Fixpoints.

The Period is independent from the initial Values.¶

mathematical approach:¶

with equating $f_x$ and $f_y$ one gets a new Definition, what can be integrated for getting the Potential.
$n \cdot f_x + b \cdot f_y = a \cdot n \cdot x - b \cdot m \cdot y$

$\frac {m} {x} \cdot f_x + \frac {a} {y} \cdot f_y = a \cdot n \cdot x- b \cdot m \cdot y$

$ n \cdot f_y - \frac {m} {y} \cdot f_y + b \cdot f_x - \frac {a} {y} \cdot f_x = 0$

$V(x,y):= n \cdot y-m \cdot ln(Y)+ b \cdot y - a \cdot ln(x)=const$

$ \frac {dV(x,y)} {dt} = 0$

=> V(x,y) is an inviariant of motion

Limiting Factor¶

For a more real Model you can initialize a Koefficient s, which is used in

$ \frac {dx} {dt} = f_x(x,y) = a \cdot x -b \cdot x\cdot y - x^2 \cdot s $

and works as a limiting factor or a Damping , what means, that you can't have as many Preys as you want.

A limiting factor is a variable of a system that causes a noticeable change in output or another measure of a type of system, here used as Resource Limitation via the derivation of Preys. With $s > 0 $ the Oszillation begins to be damped, the Phase diagram decreases, the Fixpoint won't be reached.

Theory¶

the Runge-Kutta Methods:¶

$ \frac {dx} {dt} = f_x(x,y) = a \cdot x -b \cdot x\cdot y - x^2 \cdot s $

$ \frac {dy} {dt} = f_y(x,y) = -m \cdot y +n * x\cdot y $

$ y(t_0) = y_0 $

$ x_{n+1} = x_n + \frac{1} {6} \cdot h \cdot (k_1 + 2 k_2 + 2 k_3 + k_4) $

So $x_{t+1}$ is the sum of $x_n$ and the weighted average of $k_{1-4}$, where

$ k_{1,x} = f_x \left(x_n, y_n \right) $ is the slope at the beginning of the interval, using $y$

$ k_{2,x} = f_x \left(x_n + \frac {h} {2} \cdot k_{1,x}, y_n + h \frac {k_{1,y}} {2}\right)$ is the slope at the midpoint of the interval, using $y$ and $k_2$

$ k_{3,x} = f_x\left(x_n + \frac {h} {2} \cdot k_{2,x}, y_n + h \frac {k_{2,y}} {2}\right) $ is again the slope of the midpoint, but using $y$ and $k_2$

$ k_{4,x} = f_x\left(x_n +h \cdot k_{3,x}, y_n + h \cdot k_{3,y}\right) $ is the slope of the end of the interval, using $y$ and $k_3$

we use $ h = 0.1 $ and $n = 800 $

The System has the Fixpoints

$ x = \frac {m} {g} $ $ y = \frac {a} {b} $

where¶

$x$ is the number of Prey (for example, rabbits)

$y$ is the number of Predator (for example, foxes)

$a$ is the Reproduction Koefficient for Prey

$b$ ist the Mortality Koefficient for Prey

$m$ is the Reproduction Koefficient for Predator

$g$ is the Mortality Koefficient for Predator

$h$ is the increment between $x_n$ and $x_{n+1}$

$s$ is a disorder for Resource Limitation

Where does the Method come from?:¶

Derivation of the second order¶

$k_1 = h\cdot f\left(y_n, t_n\right)$

$k_2 = h\cdot f\left(y_n + ß\cdot k_1, t_n + a\cdot h\right)$

$y_{n+1}= y_n + a\cdot k_1 + b\cdot k_2 $

The $ Taylor$ $series$ $expansion $ of $y$ in the neighborhood of $t_n$ is

$y(t_{n+1})= y(t_n) + h\cdot \frac {dy} {dt}|_{t_n} + \frac {h^2} {2} \cdot \frac {d^2y} {dt^2} |_{t_n} + O(h^3) $

with $ \frac {dy} {dt} = f(y,t) $ we can say

$ \frac {df(y,t)} {dt} = \frac {df} {dt} + \frac {df} {dy} \frac {dy} {dt} = \frac {df} {dt} + f\cdot \frac {df} {dy}$

so we get

$y_{n+1} = y_n + h\cdot f(y_n,t_n) + \frac {h^2} {2} \left[ \frac {df} {dt} + f\cdot \frac {df} {dy}\right] (y_n, t_n) + O(h^3)$

$k_2$ can be expandet correct to $O(h^3)$ as

$k_2 = h\cdot f(y_n + ß\cdot k_1, t_n + gh)$

$ = h \cdot \left(f \left( y_n,t_n \right)+g \cdot h \cdot \frac {df} {dy} \left(y_n, t_n\right) + ß\cdot k_1 \frac {df} {dy}\left(y_n, t_n \right) \right) + O(h^3) $

Now substituting for $k_2$ gives

$ y_{n+1} = y_n + (a+b)\cdot h\cdot f\left(y_n, t_n\right) + b\cdot h^2 \cdot \left(gV \frac {df} {dy} +b f \cdot \frac {df} {dy}\right)\left(y_n, t_n\right) + O(h^3)$

what gives us

$a+b=1$

$g*b= 0.5$

$ß \cdot b=0.5$

With the Choice of $g=ß=1$ , $a=b= 0.5$ we have the classical second order Runge Kutta Method:

$k_1 = h \cdot f(y_n, t_n)$

$k_2 = h \cdot f(y_n + k_1, t_n +h)$

$y_{n+1}= y_n + \frac {k_1+k_2}{ 2} $ , Second Order Runge-Kutta Method

Higher Orders can be developed in a similar fashion¶

References¶

[1] Runge Kutta Methods https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node5.html

[2] Runge Kutta Methods, Uni Münster https://www.uni-muenster.de/imperia/md/content/physik_tp/lectures/ss2017/numerische_Methoden_fuer_komplexe_Systeme_II/rkm-1.pdf

[3] 'Predator vs Prey' Programmierprojekt Prof. Wagner & private Notes https://itp.uni-frankfurt.de/~mwagner/teaching/C_WS19/projects/Predator_Prey_proj.pdf

[4] Runge Kutta Methods, Wikipedia https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods , https://de.wikipedia.org/wiki/Lotka-Volterra-Gleichungen