Introduction:

Predator-prey dynamics have been a fascinating subject of study in ecological systems, and understanding their behavior is crucial for managing biodiversity. In this blog post, we delve into numerical methods to simulate and analyze a predator-prey model known as the Lotka-Volterra equations. We’ll explore the Euler method, Modified Euler method, Midpoint method, Heun’s method, and the Runge-Kutta method, showcasing their applications and differences in solving differential equations.

The Lotka-Volterra Model:

The Lotka-Volterra model describes the population dynamics between prey and predator species. The system of differential equations for this model is given by:

\[\begin{align*} \frac{dx}{dt} &= \alpha x - \beta xy \\ \frac{dy}{dt} &= \delta xy - \gamma y \\ \end{align*}\]

Here, \(x\) and \(y\) represent the densities of prey and predator populations, respectively, and \(\alpha, \beta, \gamma, \delta\) are parameters that influence the interaction between the two species.

Numerical Methods Implementation: We implement four different numerical methods - Euler, Modified Euler, Midpoint, Heun’s, and Runge-Kutta - to solve the Lotka-Volterra equations numerically. Each method iteratively computes the populations at discrete time steps, providing insights into how the prey and predator populations evolve over time.

Results and Visualizations:

To visually compare the methods, we use a numerical example with initial conditions and parameters. The populations’ evolution is plotted over time using the Runge-Kutta method as the baseline. We observe the trajectories of prey and predator populations obtained by each numerical method, highlighting their unique characteristics.

Analysis of Numerical Methods:

We evaluate the accuracy and stability of each method by comparing their results at specific time points. A comparison table is presented, summarizing the results obtained by each method. This table provides a concise overview of the strengths and limitations of each numerical approach in solving the Lotka-Volterra equations.

Conclusion:

In conclusion, this blog post provides a comprehensive exploration of numerical methods applied to predator-prey dynamics. Numerical simulations serve as powerful tools in understanding complex ecological systems, offering a bridge between theoretical models and real-world observations.

Numerical Methods:

1. Euler Method: The Euler method is a simple numerical technique for solving ordinary differential equations. Here’s the euler_method function:

def euler_method(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha
    
    ts = [t]
    ws = [w]

    for i in range(1, N+1):
        k = h * f(t, w)
        w = w + k
        t = a + i * h
        ts.append(t)
        ws.append(w)

    return ts, ws


2. Modified Euler Method: The modified Euler method, also known as Heun’s method, improves upon the Euler method by incorporating a weighted average of two slopes. Below is the modified_euler_method function:

def modified_euler_method(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha

    ts = [t]
    ws = [w]

    for i in range(1, N+1):
        k1 = h * f(t, w)
        k2 = h * f(t + h, w + k1)
        w = w + 0.5 * (k1 + k2)
        t = a + i * h
        ts.append(t)
        ws.append(w)

    return ts, ws


3. Midpoint Method: The midpoint method refines the Euler method by evaluating the slope at the midpoint of each step. Here’s the midpoint_method function:

def midpoint_method(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha

    ts = [t]
    ws = [w]

    for i in range(1, N+1):
        k = h * f(t + h/2, w + (h/2) * f(t, w))
        w = w + k
        t = a + i * h
        ts.append(t)
        ws.append(w)

    return ts, ws


4. Heun’s Method: Heun’s method is a predictor-corrector technique that combines slopes at the beginning and end of the interval. The heuns_method function is implemented as follows:

def heuns_method(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha

    ts = [t]
    ws = [w]

    for i in range(1, N+1):
        k1 = h * f(t, w)
        k2 = h * f(t + (h/3), w + (1/3) * k1)
        k3 = h * f(t + (2*h/3), w + (2/3) * k2)
        w = w + (1/4) * (k1 + 3*k3)
        t = a + i * h
        ts.append(t)
        ws.append(w)

    return ts, ws


5. Runge-Kutta Method: The Runge-Kutta method is a robust and accurate numerical technique. Here’s the runge_kutta function:

def runge_kutta(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha

    ts = [t]
    ws = [w]

    for i in range(1, N+1):
        K1 = h * f(t, w)
        K2 = h * f(t + h/2, w + K1/2)
        K3 = h * f(t + h/2, w + K2/2)
        K4 = h * f(t + h, w + K3)
        w = w + (K1 + 2*K2 + 2*K3 + K4) / 6
        t = a + i * h
        ts.append(t)
        ws.append(w)

    return ts, ws


Comparison and Visualization: To compare results, we create a dataframe and plot prey and predator densities over time using the Runge-Kutta method as a reference:

import pandas as pd
import matplotlib.pyplot as plt

ts,ws = runge_kutta(lv, a, b, N, alpha)

x1 = [te[0] for te in ws]
x2 = [te[1] for te in ws]
plt.plot(ts,x1,label='prey')
plt.plot(ts,x2,label='predator');
plt.ylabel('Density')
plt.xlabel('Time');
plt.legend();


df = pd.DataFrame(columns=["Time", "Euler", "Modified Euler", "Midpoint", "Heuns", "Runge-Kutta"])
ts,ws = euler_method(lv, a, b, N, alpha)
df.Time = ts
df['Euler'] = [(round(wsi[0][0],4),round(wsi[1][0],4)) for wsi in ws]

ts,ws = modified_euler_method(lv, a, b, N, alpha)
df['Modified Euler'] = [(round(wsi[0][0],4),round(wsi[1][0],4)) for wsi in ws]

ts,ws = midpoint_method(lv, a, b, N, alpha)
df['Midpoint'] = [(round(wsi[0][0],4),round(wsi[1][0],4)) for wsi in ws]

ts,ws = heuns_method(lv, a, b, N, alpha)
df['Heuns'] = [(round(wsi[0][0],4),round(wsi[1][0],4)) for wsi in ws]

ts,ws = runge_kutta(lv, a, b, N, alpha)
df['Runge-Kutta'] = [(round(wsi[0][0],4),round(wsi[1][0],4)) for wsi in ws]

df.head(20)