# Common engineering numeric problems
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](http://kitchingroup.cheme.cmu.edu/pycse/pycse.html).

# Table lookup/interpolation

## 1d: numpy.interp

Problem: find P where T=130.

<table>
<tr><th>T</th><th>P</th></tr>
<tr><td>100</td><td>200</td></tr>
<tr><td>200</td><td>500</td></tr>
<tr><td>300</td><td>1000</td></tr>
</table>


```python
import numpy

T = [100, 200, 300]
P = [200, 500, 1000]

numpy.interp(150, T, P)
```

## 2d: scipy.interpolate.interp2d

Problem: Find the viscosity from this table for $T=230$ and $P=800$

<table>
<tr><th>T\P</th><td>200</td><td>500</td><td>1000</td></tr>
<tr><td>100</td><td>5</td><td>4</td><td>2</td></tr>
<tr><td>200</td><td>3</td><td>3</td><td>2</td></tr>
<tr><td>300</td><td>2</td><td>2</td><td>1</td></tr>
</table>


```python
import scipy.interpolate

T = [100, 200, 300]
P = [200, 500, 1000]

visc = [[5, 4, 2],
        [3, 3, 2],
        [2, 2, 1]]
```


```python
interpolator = scipy.interpolate.interp2d(T, P, visc)
interpolator(230, 800)
```




    array([ 2.1])



# Plotting

Plotting is normally handled by the matplotlib library. The pyplot interface aims for rough Matlab equivalence


```python
import matplotlib.pyplot as plt
```

The notebook allows for inline graphics


```python
%matplotlib inline
```


```python
plt.plot(T, P)
plt.xlabel('T')
plt.ylabel('P')
```


```python
plt.contour(T, P, visc)
```

# Algebra (manipulation of equations)

The sympy library supplies the ability to manipulate equations.


```python
import sympy
sympy.init_printing()  # This will make the expressions be printed in a more pretty way
```


```python
x = sympy.Symbol('x')
```


```python
equation = sympy.expand((x + 1)**4)
equation
```


```python
y = sympy.Symbol('y')
```


```python
equation.subs(x, y + 1)
```

# Systems of linear equations
Systems of linear equations can be solved directly without iteration. In both cases we reason via matrix algebra.

$$\begin{align}
2a + 3b + 4c &= 5 \\
2b + c &= 1 \\
a + b + c &= 1 \\
\end{align}$$

## Analytic solution
### sympy.solve


```python
a, b, c = sympy.symbols('a, b, c')
```

We state the equations in the form $f(a, b, c) = 0$.


```python
sympy.solve([2*a + 3*b + 4*c - 5,
             2*b - c - 1,
             a + b + c - 1], [a, b, c])
```

### scipy.linalg.solve


```python
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:


```python
A.I*c
```




    matrix([[-1.],
            [ 1.],
            [ 1.]])



For larger matrices, it can be faster to use a dedicated solver:


```python
scipy.linalg.solve(A, c)
```




    array([[-1.],
           [ 1.],
           [ 1.]])



# Laplace transforms
## sympy.laplace_transform


```python
t, s = sympy.symbols('t, s')
```


```python
sympy.laplace_transform(sympy.sin(t), t, s)
```




    (1/(s**2 + 1), 0, True)



## sympy.inverse_laplace_transform


```python
sympy.inverse_laplace_transform(1/(s + 1), s, t)
```

# Systems of nonlinear equations
## Analytic solution (sometimes):
### sympy.solve


```python
sympy.solve([a**2 + b,
             a/b + 1], [a, b])
```

## Numeric solution (often)
### scipy.optimize.fsolve


```python
import scipy.optimize
```


```python
def eqs(x):
    a, b = x
    return [a**2 + b,
            a/b + 1]
```


```python
scipy.optimize.fsolve(eqs, [2, -1])
```




    array([ 1., -1.])



# Systems of ODEs (IVP)
## scipy.integrate.odeint


```python
def sys(x, t):
    a, b = x
    return [-2*a, -a - 2*b]
```


```python
times = numpy.linspace(0, 5)
```


```python
scipy.integrate.odeint?
```


```python
solution = scipy.integrate.odeint(sys, [1, 1], times)
```


```python
plt.plot(times, solution)
```

# Numerical integration (aka quadrature)
## Given a function
### scipy.integrate.quad

Let's say we're trying to integrate the surface under the curve $f(x) = \sin(x)/x$ between 0.5 and 1


```python
def f(x):
    return numpy.sin(x)/x
```


```python
x = numpy.linspace(0.01, 10)
```


```python
plt.plot(x, f(x))
```


```python
integral, error = scipy.integrate.quad(f, 0.5, 1)
print(integral)
```

    0.45297565232411635


## Given smooth data
### Trapezoidal rule: scipy.signal.trapz


```python
xdata = numpy.linspace(0.5, 1)
ydata = f(xdata)
```


```python
scipy.integrate.trapz(ydata, xdata)  # note it's _not_ x, y
```

### Simpson's rule: scipy.signal.simps


```python
scipy.integrate.simps(ydata, xdata)
```

# Numerical derivatives
## Given smooth data
### central difference estimate: numpy.gradient


```python
x = numpy.linspace(0, 2*numpy.pi)
```


```python
y = numpy.sin(x)
```


```python
derivative_estimate = numpy.gradient(y, x[1])
```


```python
plt.plot(x, y, x, derivative_estimate)
```

## Given noisy data
### Savitzky-Golay scipy.signal.savgol_filter

That's fine for smooth data, but what if you have some noise added?


```python
y_noise = y + numpy.random.randn(len(y))*0.1
```


```python
plt.plot(x, y_noise)
```


```python
derivative_estimate = numpy.gradient(y_noise, x[1])
```


```python
plt.plot(x, y_noise, x, derivative_estimate)
```

That's clearly not great.


```python
import scipy.signal
```


```python
smooth_derivative_estimate = scipy.signal.savgol_filter(y, 5, 2, deriv=1, delta=x[1])
```


```python
plt.plot(x, y_noise, x, smooth_derivative_estimate)
```

That's more like it!

# Optimisation
## Nonlinear
### scipy.optimize.minimize

Let's minimise $\sin(x)/x$ 


```python
plt.plot(x, f(x))
```

By inspection, we can see the minimum lies between 4 and 5.


```python
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])



# Statistical analysis
## Regression

Simple polynomial regression


```python
x = numpy.linspace(1, 4, 20)
```


```python
y = f(x)
```


```python
plt.plot(x, y, '.')
```


```python
quadratic = numpy.polyfit(x, y, 2)
```


```python
smoothx = numpy.linspace(1, 4, 1000)
```


```python
plt.plot(x, y, '.')
plt.plot(smoothx, numpy.polyval(quadratic, smoothx))
```

181b48fed3d86058bda4349d5f9e4ce0c68fbc23
ipynb
Jupyter Notebook
Useful references/Common problems solved in Python.ipynb
afcarl/Useful-python
5d1947052fb25b2388704926e4692511cc162031
MIT
Useful references/Common problems solved in Python.ipynb
afcarl/Useful-python
5d1947052fb25b2388704926e4692511cc162031
MIT
Useful references/Common problems solved in Python.ipynb
afcarl/Useful-python
5d1947052fb25b2388704926e4692511cc162031
MIT
2018-09-05T21:48:57.000Z
2018-09-05T21:48:57.000Z
Qwen/Qwen-72B
1. YES
2. YES