Modeling and Simulation in Python

Implementation of Lotka-Volterra using Euler and run_ode_solver

Copyright 2018 Allen Downey

License: Creative Commons Attribution 4.0 International


In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline

# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import functions from the modsim.py module
from modsim import *

Euler with implicit dt=1


In [2]:
def run_simulation(system, update_func):
    """Runs a simulation of the system.
        
    system: System object
    update_func: function that updates state
    
    returns: TimeFrame
    """
    unpack(system)
    
    frame = TimeFrame(columns=init.index)
    frame.row[t0] = init
    
    for t in linrange(t0, t_end):
        frame.row[t+1] = update_func(frame.row[t], t, system)
    
    return frame

In [3]:
def update_func(state, t, system):
    """Update the Lotka-Volterra model.
    
    state: State(x, y)
    t: time
    system: System object
    
    returns: State(x, y)
    """
    unpack(system)
    x, y = state

    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    
    x += dxdt
    y += dydt
    
    return State(x=x, y=y)

In [4]:
init = State(x=1, y=1)


Out[4]:
values
x 1
y 1

In [5]:
system = System(alpha=0.05,
                beta=0.1,
                gamma=0.1,
                delta=0.1,
                t0=0,
                t_end=200)


Out[5]:
values
alpha 0.05
beta 0.10
gamma 0.10
delta 0.10
t0 0.00
t_end 200.00

In [6]:
update_func(init, 0, system)


Out[6]:
values
x 0.95
y 1.00

In [7]:
results = run_simulation(system, update_func)
results.head()


Out[7]:
x y
0.0 1 1
1.0 0.95 1
2.0 0.9025 0.995
3.0 0.857826 0.985299
4.0 0.816196 0.97129

In [8]:
results.plot()


Out[8]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa73f3c7c88>

Using the ODE solver


In [9]:
def slope_func(state, t, system):
    """Compute slopes for the Lotka-Volterra model.
    
    state: State(x, y)
    t: time
    system: System object
    
    returns: pair of derivatives
    """
    unpack(system)
    x, y = state

    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    
    return dxdt, dydt

In [10]:
system.set(init=init, t_end=200)

In [11]:
results, details = run_ode_solver(system, slope_func)
details


Out[11]:
values
sol None
t_events []
nfev 134
njev 0
nlu 0
status 0
message The solver successfully reached the end of the...
success True

In [12]:
results.plot()


Out[12]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa72c15b9b0>

In [13]:
system.set(init=init, t_end=200)
results, details = run_ode_solver(system, slope_func, max_step=2)
details


Out[13]:
values
sol None
t_events []
nfev 608
njev 0
nlu 0
status 0
message The solver successfully reached the end of the...
success True

In [14]:
results.plot()


Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa72c06cf60>

Euler, running for a longer duration


In [15]:
system.set(t_end=2000)
results = run_simulation(system, update_func)
results.head()


Out[15]:
x y
0.0 1 1
1.0 0.95 1
2.0 0.9025 0.995
3.0 0.857826 0.985299
4.0 0.816196 0.97129

In [16]:
results.plot()


Out[16]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa72c052240>

Running the ODE solver for a longer duration


In [17]:
results, details = run_ode_solver(system, slope_func)
details


Out[17]:
values
sol None
t_events []
nfev 1172
njev 0
nlu 0
status 0
message The solver successfully reached the end of the...
success True

In [18]:
results.plot()


Out[18]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa72bfd2630>

In [ ]: