Tutorial 5 for WWCS2016

By nkoub / December 17, 2015 / Bistable systems, Networks, software, workshops / No Comments

IPython Notebooks for the Tutorial 5 of the Winter Workshop on Complex Systems 2016

This python notebook solves a bistable system on a complex network and visualizes the solution

Import networkx library

import networkx as nx

For the numerical integration of the system we use the ode routine from scipy.integrate

import numpy as np
from scipy.integrate import ode

we plot and animate the results with matplotlib

import matplotlib.pyplot as plt
from matplotlib import animation

Define the right-hand-side of the bistable system’s ODE f(x)=x(h-x)(x-a)

def systemRHS(t0, y0, params):
    laplacian = params[0]
    a = params[2]
    h = params[3]
    diff = coupling * laplacian * y0 
    f = y0 * (h - y0) * (y0 - a) + diff
    return f

System’s parameters

h = 0.25           # the unstable state
a = 1.0            # the stable state
coupling = 0.05    # the diffusive coupling

Create the network. Here we show a regular tree with branching ratio 2 and 5 shells

graph = nx.balanced_tree(2,3) 
graph.add_edges_from([(5,51),(5,21),(5,53),(5,54),(5,55)])

Derive the Laplacian matrix

lapl = (-1.0) * nx.laplacian_matrix(graph) # nx: L = D - A


mypos = nx.graphviz_layout(graph,'neato',root=0)
fig = plt.figure(0)
ax0 = fig.add_subplot(111)
ax0.axis('off')
nx.draw_networkx(graph,ax=ax0,pos=mypos,
                 node_size=700, linewidths=1.,
                 node_color='white',
                 with_labels=True,
                 edge_color='green',
                 width=2,alpha=1)
plt.show()

Initial conditions

init_u = [0]                              # active state
u0 = np.zeros(graph.number_of_nodes())    # passive state
for i in init_u:
    u0[i] = a

Integrator parameters

t0 = 0.0
tmax = 200
dt = 0.1
params = (lapl,coupling,a,h) # tuple with 4 elements

Integrate the system

Declare the integrator

sol = [u0]
solver = ode(systemRHS).set_f_params(params)
solver.set_integrator('dopri5')  
solver.set_initial_value(u0,t0)

Run the integrator

while solver.successful() and solver.t < tmax:
    solver.integrate(solver.t+dt)
    sol.append(solver.y)
sol = np.array(sol)

Plot the results

Animation of the solution by sampling every ::n snapshots

sol_anim = sol[::10]
sol_anim.shape

Define the animation function

fig_anim = plt.figure()
ax_anim = fig_anim.add_subplot(111)

def animate(t):
    ax_anim.cla()
    ax_anim.axis('off')
    ax_anim.set_title('t={}'.format(t))
    pl = nx.draw_networkx(graph,ax=ax_anim,pos=mypos,
                          vmin=0.0,vmax=1.0,
                          node_size=200, linewidths=1.,
                          node_color=sol_anim[t],
                          cmap=plt.cm.binary,
                          with_labels=False,
                          edge_color='green',
                          width=2,alpha=1)

Call the animate function

matplotlib inline
from tempfile import NamedTemporaryFile


VIDEO_TAG = """<video controls>
<source src="data:video/x-webm;base64,{0}" type="video/webm">
Your browser does not support the video tag.
</video>"""

def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with NamedTemporaryFile(suffix='.webm') as f:
            anim.save(f.name, fps=25, extra_args=['-vcodec', 'libvpx'])
            video = open(f.name, "rb").read()
        anim._encoded_video = video.encode("base64")
    
    return VIDEO_TAG.format(anim._encoded_video)


from IPython.display import HTML

def display_animation(anim):
    plt.close(anim._fig)
    return HTML(anim_to_html(anim))


display_animation(anim)