This notebook aims to collect a very brief summary of simple ways to solve common engineering problems in Python
For a far more comprehensive discussion of chemical engineering problems solved in Python, see the Kitchin Group pycse file.
Problem: find P where T=130.
T | P |
---|---|
100 | 200 |
200 | 500 |
300 | 1000 |
import numpy
T = [100, 200, 300]
P = [200, 500, 1000]
numpy.interp(150, T, P)
Problem: Find the viscosity from this table for $T=230$ and $P=800$
T\P | 200 | 500 | 1000 |
---|---|---|---|
100 | 5 | 4 | 2 |
200 | 3 | 3 | 2 |
300 | 2 | 2 | 1 |
import scipy.interpolate
T = [100, 200, 300]
P = [200, 500, 1000]
visc = [[5, 4, 2],
[3, 3, 2],
[2, 2, 1]]
interpolator = scipy.interpolate.interp2d(T, P, visc)
interpolator(230, 800)
array([ 2.1])
Plotting is normally handled by the matplotlib library. The pyplot interface aims for rough Matlab equivalence
import matplotlib.pyplot as plt
The notebook allows for inline graphics
%matplotlib inline
plt.plot(T, P)
plt.xlabel('T')
plt.ylabel('P')
<matplotlib.text.Text at 0x113e10128>
plt.contour(T, P, visc)
<matplotlib.contour.QuadContourSet at 0x11419d6d8>
The sympy library supplies the ability to manipulate equations.
import sympy
sympy.init_printing() # This will make the expressions be printed in a more pretty way
x = sympy.Symbol('x')
equation = sympy.expand((x + 1)**4)
equation
y = sympy.Symbol('y')
equation.subs(x, y + 1)
a, b, c = sympy.symbols('a, b, c')
We state the equations in the form $f(a, b, c) = 0$.
sympy.solve([2*a + 3*b + 4*c - 5,
2*b - c - 1,
a + b + c - 1], [a, b, c])
A = numpy.matrix([[2, 3, 4],
[0, 2, -1],
[1, 1, 1]])
c = numpy.matrix([5, 1, 1]).T # .T for transpose
Direct calculation via inverse:
A.I*c
matrix([[-1.], [ 1.], [ 1.]])
For larger matrices, it can be faster to use a dedicated solver:
scipy.linalg.solve(A, c)
array([[-1.], [ 1.], [ 1.]])
t, s = sympy.symbols('t, s')
sympy.laplace_transform(sympy.sin(t), t, s)
(1/(s**2 + 1), 0, True)
sympy.inverse_laplace_transform(1/(s + 1), s, t)
sympy.solve([a**2 + b,
a/b + 1], [a, b])
import scipy.optimize
def eqs(x):
a, b = x
return [a**2 + b,
a/b + 1]
scipy.optimize.fsolve(eqs, [2, -1])
array([ 1., -1.])
def sys(x, t):
a, b = x
return [-2*a, -a - 2*b]
times = numpy.linspace(0, 5)
scipy.integrate.odeint?
solution = scipy.integrate.odeint(sys, [1, 1], times)
plt.plot(times, solution)
[<matplotlib.lines.Line2D at 0x113e2b6d8>, <matplotlib.lines.Line2D at 0x113e2b860>]
Let's say we're trying to integrate the surface under the curve $f(x) = \sin(x)/x$ between 0.5 and 1
def f(x):
return numpy.sin(x)/x
x = numpy.linspace(0.01, 10)
plt.plot(x, f(x))
[<matplotlib.lines.Line2D at 0x114550b70>]
integral, error = scipy.integrate.quad(f, 0.5, 1)
print(integral)
0.45297565232411635
xdata = numpy.linspace(0.5, 1)
ydata = f(xdata)
scipy.integrate.trapz(ydata, xdata) # note it's _not_ x, y
scipy.integrate.simps(ydata, xdata)
x = numpy.linspace(0, 2*numpy.pi)
y = numpy.sin(x)
derivative_estimate = numpy.gradient(y, x[1])
plt.plot(x, y, x, derivative_estimate)
[<matplotlib.lines.Line2D at 0x114a4ee80>, <matplotlib.lines.Line2D at 0x114a8c8d0>]
That's fine for smooth data, but what if you have some noise added?
y_noise = y + numpy.random.randn(len(y))*0.1
plt.plot(x, y_noise)
[<matplotlib.lines.Line2D at 0x114b7a400>]
derivative_estimate = numpy.gradient(y_noise, x[1])
plt.plot(x, y_noise, x, derivative_estimate)
[<matplotlib.lines.Line2D at 0x114cadf98>, <matplotlib.lines.Line2D at 0x114cb39e8>]
That's clearly not great.
import scipy.signal
smooth_derivative_estimate = scipy.signal.savgol_filter(y, 5, 2, deriv=1, delta=x[1])
plt.plot(x, y_noise, x, smooth_derivative_estimate)
[<matplotlib.lines.Line2D at 0x114de9a90>, <matplotlib.lines.Line2D at 0x114de9c50>]
That's more like it!
Let's minimise $\sin(x)/x$
plt.plot(x, f(x))
/Users/alchemyst/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in true_divide from ipykernel import kernelapp as app
[<matplotlib.lines.Line2D at 0x114153860>]
By inspection, we can see the minimum lies between 4 and 5.
scipy.optimize.minimize(f, 4.5)
fun: -0.21723362821071318 hess_inv: array([[ 4.61057157]]) jac: array([ -4.67523932e-07]) message: 'Optimization terminated successfully.' nfev: 12 nit: 3 njev: 4 status: 0 success: True x: array([ 4.49340729])
Simple polynomial regression
x = numpy.linspace(1, 4, 20)
y = f(x)
plt.plot(x, y, '.')
[<matplotlib.lines.Line2D at 0x1165d26d8>]
quadratic = numpy.polyfit(x, y, 2)
smoothx = numpy.linspace(1, 4, 1000)
plt.plot(x, y, '.')
plt.plot(smoothx, numpy.polyval(quadratic, smoothx))
[<matplotlib.lines.Line2D at 0x11670cf98>]