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.

# nbi:hide_in

import numpy as np

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

from IPython.display import clear_output, display

import matplotlib.pyplot as plt
%matplotlib inline


def rk4(fx, fy, x0, y0, x1, n):
    
    vx = [0] * (n + 1)
    vy = [0] * (n + 1)
    vi = [0] * (n + 1)
    
    vx[0] = x = x0
    vy[0] = y = y0
    for i in range(1, n + 1):
    
        kx1 = fx(x, y)
        ky1 = fy(x, y)

        kx2 = fx(x + 0.5 * h * kx1, y + 0.5 * h * ky1)
        ky2 = fy(x + 0.5 * h * kx1, y + 0.5 * h * ky1)

        kx3 = fx(x + 0.5 * h * kx2, y + 0.5 * h * ky2)
        ky3 = fy(x + 0.5 * h * kx2, y + 0.5 * h * ky2)

        kx4 = fx(x + h * kx3, y + h * ky3)
        ky4 = fy(x + h * kx3, y + h * ky3)

        vx[i] = x = x + h * (kx1 + 2.*kx2 + 2.*kx3 + kx4) / 6.
        vy[i] = y = y + h * (ky1 + 2.*ky2 + 2.*ky3 + ky4) / 6.
        vi[i] = h*i    
    return vx, vy, vi

def fx(x, y):
    return (a-b*y)*x

def fy(x, y):
    return (-m+g*x)*y

def update(a_par=0.5,b_par=0.5, m_par=0.5, g_par=0.5):
    global a
    global b
    global m
    global g
    global e
    
    
    a = a_par
    b = b_par
    m = m_par
    g = g_par
    e = (a*m)**.5
    
    
    
    vx, vy, vi = rk4(fx, fy, 2, 1, 10, 500)
    
    fig = plt.figure(figsize =(10,6))##
    ax = fig.add_subplot(1, 2, 1)

    ax.set_ylabel('Population (blue=Prey,orange=Predator)')
    ax.set_xlabel('Time')

    ax2 = fig.add_subplot(1, 2, 2)
    ax2.set_ylabel('Population Predator')
    ax2.set_xlabel('Population Prey')
    
    

    line1, = ax.plot(vi, vx)
    line2, = ax.plot(vi,vy)

    dot1, = ax2.plot([m/g], [a/b], 'bo')

    line3, = ax2.plot( vx, vy, c='g')
    line4, = ax2.plot (2, 1, 'go')
    
    ax2.text(1., 5., "Eigenvalues: $\\pm i \\cdot$" + "{:.2f}".format(e), fontsize=12)
    
    ax.set_xlim([0.,50.])
    ax.set_ylim([0.,15.])
    
    ax2.set_xlim([0.,6.])
    ax2.set_ylim([0.,6.])
    
    
    
    #clear_output(True)
    
    #print ("Eigenvalues:") 
    #print("+/- i*",e,)

    

a = 0.5 #ReproduktionPrey
b = 0.1 #Mortalityprey
g = 0.5 #ReproduktionPredator
m = 0.5 #MortalityPredator
h = 0.1

widgets.interact(update,a_par=(0.,1.),
                 b_par=(0.,1.),
                 m_par=(0.,1.),
                 g_par=(0.,1.))
<function __main__.update(a_par=0.5, b_par=0.5, m_par=0.5, g_par=0.5)>

Theory

the Runge-Kutta Methods:

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

$ \frac {dy} {dt} = f_y(x,y) = -m \cdot y +n \cdot 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}$

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

Constant of motion

constant%20of%20motion.PNG

You can see the known phase diagram with different initial Values, and every initial value leads to periodic solutions

The Periodicity is independent from the initial Values.

mathematical approach:

$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