#least square regression
import numpy as np
X = np.column_stack([x, x**0])
XTX = np.dot(X.T, X)
XTy = np.dot(X.T, y)
p = np.dot(np.linalg.inv(XTX), XTy)
slope, intercept = p # note the order in X
print ’The slope is {0} \nand intercept is {1}’.format(slope, intercept)
import matplotlib.pyplot as plt
plt.plot(x, y, ’bo’)
plt.plot(x, np.dot(X, p), ’r--’)
L
alpha = 1 - 0.95
p, pint, se = regress(X, y, alpha)
slope_interval, intercept_interval = pint
#derivation
import numpy as np
np.set_printoptions(precision=3)
from pycse import deriv, regress
import matplotlib.pyplot as plt
t = np.array([0, 50, 100, 150, 200, 250, 300], dtype=np.float)
Ca = np.array([0.05, 0.038, 0.0306, 0.0256, 0.0222, 0.0195, 0.0174])
dCadt = deriv(t, Ca)
x = np.log(Ca)
y = np.log(-dCadt)
X = np.column_stack([x**0, x])
p, pint, se = regress(X, y, 0.05)
intercept_range = pint[0]
alpha_range = pint[1]
#Polyfit
import numpy as np
from scipy.optimize import curve_fit
x = np.array([0.5, 0.387, 0.24, 0.136, 0.04, 0.011])
y = np.array([1.255, 1.25, 1.189, 1.124, 0.783, 0.402])
def func(x, a, b):
return a * x / (b + x)
initial_guess = [1.2, 0.03]
pars, pcov = curve_fit(func, x, y, p0=initial_guess)
a,b = pars
17 print ’a = {0} and b={1}’.format(a,b)
18 import matplotlib.pyplot as plt
19 plt.plot(x,y,’bo ’)
20 xfit = np.linspace(min(x), max(x))
21 yfit = func(xfit, *pars)
22 plt.plot(xfit,yfit,’b-’)
#nonlinfit
import numpy as np
2 from pycse import nlinfit
3 import matplotlib.pyplot as plt
4
5 t = np.array([0, 50, 100, 150, 200, 250, 300])
6 Ca = np.array([0.05, 0.038, 0.0306, 0.0256, 0.0222, 0.0195, 0.0174])
7
8 Ca0 = 0.05
9
10 def model(t, k, alpha):
11 return 1.0 / (-k * t +
12 k * t * alpha +
13 np.exp(-alpha * np.log(Ca0)) * Ca0)**(1.0/(alpha-1.0))
14
15 guess = [0.1, 2.0] # from the linear regression
16
17 p, pint, se = nlinfit(model, t, Ca, guess, 0.05)
18
19 k_range = pint[0]
20 alpha_range = pint[1]
21
22 print ’alpha = {0} \n at the 95% confidence level’.format(alpha_range)
23 print ’k = {0} \n at the 95% confidence level’.format(k_range)
24
25 # always visually inspect the fit
26 tfit = np.linspace(min(t), max(t))
27 plt.plot(t, Ca,’ko ’)
28 plt.plot(tfit, model(tfit, *p))
#monte carlo
import numpy as np
2 from scipy.optimize import fsolve
3
4 N = 10000 # number of MC samples
5 V = 66000 # L
6
7 Fa0 = np.random.normal(5, 0.05, (1, N))
8 v0 = np.random.normal(10.0, 0.1, (1, N))
9 k = np.random.normal(3.0, 0.2, (1, N))
10
11 # create the array to store the results in
12 # It is usually more efficient to create arrays then fill them in.
13 SOL = np.empty(k.shape)
14
15 for i in range(N):
16 def func(Ca):
17 ra = -k[0,i] * Ca**2
18 return Fa0[0,i] - v0[0,i] * Ca + V * ra
19 guess = 0.1 * Fa0[0,i] / v0[0,i] # guessing 90% conversion
SOL[0,i] = fsolve(func, guess)[0]
#uncertainty
import uncertainties as u
2 from scipy.optimize import fsolve
3
4 V = 66000 # reactor volume L^3
5 Fa0 = u.ufloat(5.0, 0.05) # mol / h
6 v0 = u.ufloat(10., 0.1) # L / h
7 k = u.ufloat(3.0, 0.2) # rate constant L/mol/h
8
10 def func(Ca, v0, k, Fa0, V):
11 "Mole balance for a CSTR. Solve this equation for func(Ca)=0"
12 Fa = v0 * Ca # exit molar flow of A
13 ra = -k * Ca**2 # rate of reaction of A L/mol/h
14 return Fa0 - Fa + V * ra
15
16
17 def Ca_solve(v0, k, Fa0, V):
21 guess = 0.1 * Fa0 / v0
22 sol = fsolve(func, guess, args=(v0, k, Fa0, V))[0]
23 return sol
24
25
26 Ca_exit = u.wrap(Ca_solve)(v0, k, Fa0, V)
27 print ’The exit concentration is {0}’.format(Ca_exit)
#intermidates
1 import numpy as np
2 from scipy.integrate import odeint
3
4 k1 = 1.5e-3 #/u.s;
5 k2 = 2.3e6 #*u.L/u.mol/u.s;
6 k3 = 5.7e4 #/u.s;
7 k4 = 9.5e8 #*u.L/u.mol/u.s;
8 k5 = 2.0e9 #*u.L/u.mol/u.s;
9
10 def dCdt(C, t):
11 Ca, Cb, Cc, Cd, Ce = C
12
13 r1 = k1 * Ca
14 r2 = k2 * Ca * Cb
15 r3 = k3 * Cc
16 r4 = k4 * Ca * Cd
17 r5 = k5 * Cc**2
18
19 ra = -r1 - r2 - r4
20 rb = 2*r1 - r2
21 rc = r2 - r3 + r4 - 2*r5
22 rd = r3 - r4
23 re = r3 # this is the ethylene
24
25 dCadt = ra
26 dCbdt = rb
27 dCcdt = rc
28 dCddt = rd
29 dCedt = re
30 return [dCadt, dCbdt, dCcdt, dCddt, dCedt]
31
32 C0 = [0.1, 0, 0, 0, 0]
33
34 tspan = np.linspace(0, 12)
35 sol = odeint(dCdt, C0, tspan, rtol=1e-12)
36
37 import matplotlib.pyplot as plt
38 plt.plot(tspan, sol)
39 plt.legend([’A’, ’B’, ’C’ ,’D’, ’E’], loc=’best’)
40 plt.xlabel(’Time (s)’)
41 plt.ylabel(’C (mol/L)’)
42 plt.savefig(’images/ethane-cracking-1.png’)
43
44 # now plot the intermediates
45 plt.figure()
46 plt.semilogy(tspan, sol[:,1:4])
47 plt.legend([’B’, ’C’ ,’D’], loc=’best’)