USING BOKEH TO RECONSTRUCT THE BRUSSELATOR OSCILLATIONS NUMERICALLY (ODEINT): A WORK IN PROGRESS

SEE HERE FOR MORE INFORMATION: http://www.idea.wsu.edu/OscilChem/#Brusselator%20Model



In [1]:

    
%matplotlib inline
import numpy as np
from scipy.integrate import odeint
from pylab import *
from bokeh.plotting import *
from bokeh.models import ColumnDataSource
#from bokeh.io import output_notebook

# DEFINE A FUNCTION THAT OUTPUTS CHEMICAL SPECIES CONCENTRATION TIME SERIES, 
#GIVEN AN ARRAY OF TIME POINTS (x) AND SOME REACTION RATES (A,B)

def g(x,A,B):
    #A = 1
    #B = 3
    #DEFINE OUR SYSTEM OF DIFFERENTIAL EQUATIONS dZ/dt = f(Z, t)
    def f(Z,t):
        X = Z[0]
        Y = Z[1]
        #the model equations
        f0 = A + np.power(X,2)*Y - B*X - X
        f1 = B*X - np.power(X,2)*Y
        return [f0 , f1]
    #initial conditions
    X0 = 1.1
    Y0 = 3
    Z0 = [ X0 , Y0 ]
    #t  = np.linspace(0, 50, 1000)   # time grid
    #Solve the DEs
    return odeint( f , Z0 , x)
#GENERATE THE FIRST SOLUTION
x = np.linspace(0, 50, 1000)
soln = g(x , 1 , 3)
y = soln[:,0]
z = soln[:,1]



In [2]:

    
#LOAD THE NOTEBOOK
output_notebook()
#GET ALL DATA TOGETHER TO DEFINE THE 1ST PLOTS
source = ColumnDataSource( data = dict(x=x,y=y,z=z))
#output_file("lines1.html")
p1 = figure( width=350, height=350, title='Chemical species Y')
p1.line('x','y', line_width = 4 , source = source )
p1.xaxis.axis_label = 'time'
p1.yaxis.axis_label = 'concentration'
p2 = figure(width=350, height=350, title='Chemical species Z')
p2.line('x','z', line_width = 4 , source = source , color = 'red')
p2.xaxis.axis_label = 'time'
p2.yaxis.axis_label = 'concentration'
p3 = figure( width=350, height=350, title='Phase plane')
p3.line('y','z', line_width = 2 , source = source , color = 'green')
p3.yaxis.axis_label = '[Z]'
p3.xaxis.axis_label = '[Y]'
p = gridplot([[p1 , p2 ], [p3]])

#THIS IS HOW WE UPDATE THE PLOTS
def update(A=1 , B =3):
    source.data['y'] = g(x , A, B)[:,0]
    source.data['z'] = g(x , A, B)[:,1]
    source.push_notebook()
#LETS LOOK AT THE CONCENTRATION PROFILES AND PHASE PLANE NOW 
show(p)









    




    
        
        
        
    
        
        BokehJS successfully loaded.



In [3]:

    
#HERE WE GENERATE SOME SLIDERS THAT ALLOW US TO CHANGE THE REACTION RATES A,B ON THE FLY
from IPython.html.widgets import interact
interact( update , A = (1,5,0.1), B = (2,10,0.1))
#output_file("lines1.html")



In [ ]: