# Modeling and Simulation in Python

Implementation of Lotka-Volterra using Euler and run_ode_solver

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

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)

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)

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 [ ]: