## This sheet simulates two reversible first-order chemical reactions in series:
## A1 <-> A2 <-> A3
## The initial value problem is solved by Runge-Kutta and the dynamics are graphed
var('A10 A20 A30 k1 k2 k_1 k_2')
var('A1 A2 A3 t')
A10 = 0.01
A20 = 0.00
A30 = 0.00
k1 = 1
k_1 = 1
k2 = 1
k_2 = 1
## Calculation Parameters
end_points = 5
stepsize = 0.01
steps = end_points/stepsize
ics = [0,A10,A20,A30]
## Equations
r1 = (diff(A1,t) == k_1*A2 - k1*A1)
r2 = (diff(A2,t) == k1*A1 - k_1*A2 - k2*A2 + k_2*A3)
r3 = (diff(A3,t) == k2*A2 - k_2*A3)
des = [r1.rhs(), r2.rhs(), r3.rhs()]
sol = desolve_system_rk4(des,[A1,A2,A3],ics,ivar=t,end_points=end_points,step=stepsize)
|
|
## Process the output into a form that can be graphed
sols_1=[]
sols_2=[]
sols_3=[]
for i in range(steps):
sols_1.append([sol[i][0],sol[i][1]])
sols_2.append([sol[i][0],sol[i][2]])
sols_3.append([sol[i][0],sol[i][3]])
################################################
#### Unnecessarily Fancy Plotting Stuff ####
################################################
## Create a plot object
a = plot([])
## Set the plot parameters
title = "Two-step, series, reversible reaction" ## Graph Title
xmin = 0 ## X minimum
xmax = end_points ## X maximum
ymin = 0 ## Y minimum
ymax = 0.01 ## Y maximum
## Add a title to the plot
a += text(title,(xmax/1.8,ymax),color='black',fontsize=15)
## Add the desired lines to the plot
a += list_plot(sols_1,color='orange',legend_label='[A1]')
a += list_plot(sols_2,color='purple',legend_label='[A2]')
a += list_plot(sols_3,color='green',legend_label='[A3]')
## For more information on plots in general, evaluate 'plot?'
## For a list of legend options, evaluate 'a.set_legend_options?'
## For a list of Sage predefined colors, evaluate 'sorted(colors)'
a.set_axes_range(xmin,xmax,ymin,ymax)
a.axes_labels(['time','concentration'])
a.axes_label_color('grey')
a.set_legend_options(ncol=3,borderaxespad=5,back_color='whitesmoke',fancybox=true)
show(a)
|
|
