# Tutorial 5 for WWCS2016

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)


### Derive the Laplacian matrix

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


mypos = nx.graphviz_layout(graph,&amp;#039;neato&amp;#039;,root=0)
fig = plt.figure(0)
ax0.axis(&amp;#039;off&amp;#039;)
nx.draw_networkx(graph,ax=ax0,pos=mypos,
node_size=700, linewidths=1.,
node_color=&amp;#039;white&amp;#039;,
with_labels=True,
edge_color=&amp;#039;green&amp;#039;,
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(&amp;#039;dopri5&amp;#039;)
solver.set_initial_value(u0,t0)


### Run the integrator

while solver.successful() and solver.t &amp;lt; 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()

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 = """&lt;video controls&gt;
&lt;source src="data:video/x-webm;base64,{0}" type="video/webm"&gt;
Your browser does not support the video tag.
&lt;/video&gt;"""

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'])
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)