Mechpy Tutorials

a mechanical engineering toolbox

source code - https://github.com/nagordon/mechpy
documentation - https://nagordon.github.io/mechpy/web/

Neal Gordon
2017-02-20

Math

Linear Algebra
Signal Processing

plotting

This is more or less my notes on using python for a variety of engineering tasks, including plotting, solving equations, linear algebra etc. Some of this work is my own, but many of it is from the scipy documentation.

Python Initilaization with module imports

In [ ]:
In [1]:
# setup 
import numpy as np
import sympy as sp
import scipy
from scipy import linalg
from pprint import pprint
sp.init_printing(use_latex='mathjax')

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12, 8)  # (width, height)
plt.rcParams['font.size'] = 14
plt.rcParams['legend.fontsize'] = 16
from matplotlib import patches

#get_ipython().magic('matplotlib')  # seperate window
get_ipython().magic('matplotlib inline') # inline plotting

Math with Python

index

WHen using ipython, variables and modules can be stored in the console. Use %reset to remove any modules and variables from the kernal

The modules used for math are numpy , scipy and matplotlib

All the these modules can be accessed with the module pylab, but is not updated as regular as the other packages.

the linear algebra system in python via numpy and scipy is a bit awkward at first, but once the nuances are understood it has most of the functionality of MATLAB with zero the cost

The matrix class of scipy is handy but limits the type of data stored. As with MATLAB all "arrays" or "matrices" are just matrices. This is not true with python. Thus, to avoid confusion, it is discouraged to use the matrix class

TO get started lets looks at numpy arrays, and then later we will use python powerful symbolic math package sympy for a different way to perform calculations

Linear Algebra with Python

index

Python's numpy package allows python, a generic computing language to perform powerful mathematical calculations. Although python's math syntax is not as obvious as MATLAB's, the functionality is comparable. This document is designed to be an intro to that syntax

Some references

http://nbviewer.ipython.org/github/carljv/cython_testing/blob/master/cython_linalg.ipynb

We can either use scipy, which includes numpy, http://docs.scipy.org/doc/

or use numpy directly http://docs.scipy.org/doc/numpy/

Since there are many ways to solve linear algebra problems, (eg Octave/Matlab, julia, scipy, numpy) I tend to prefer the most matlabesc approaches due to the ubiquity of Matlab and the simplicity of the syntax, which frankly, python suffers with.

The major difference between arrays and matrices in python is that arrays are n-dimensions, where matrices are only up to 2-dimensions
m

Numpy arrays

In [2]:
# using the range function to create a numpy array
x = np.array(range(10))
print('x = ',x, 'as a',x.dtype)

x =  [0 1 2 3 4 5 6 7 8 9] as a int32
In [3]:
# using numpy function arange
x = np.arange(10)
print('x = ',x, 'as a',x.dtype)

x =  [0 1 2 3 4 5 6 7 8 9] as a int32
In [4]:
x = np.linspace(0,9,10)
print('x = ',x, 'as a',x.dtype)

x =  [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9.] as a float64
In [5]:
type(x)
Out[5]:
numpy.ndarray
In [6]:
# multiple array operations
x*3
Out[6]:
array([  0.,   3.,   6.,   9.,  12.,  15.,  18.,  21.,  24.,  27.])
In [7]:
# exponent array operations
x**2
Out[7]:
array([  0.,   1.,   4.,   9.,  16.,  25.,  36.,  49.,  64.,  81.])
In [8]:
# trig array operations
np.sin(x)
Out[8]:
array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ,
       -0.95892427, -0.2794155 ,  0.6569866 ,  0.98935825,  0.41211849])
In [9]:
x.sum()
Out[9]:
$$45.0$$
In [10]:
print([k for k in dir(x) if k[:2] != '__' ])

['T', 'all', 'any', 'argmax', 'argmin', 'argpartition', 'argsort', 'astype', 'base', 'byteswap', 'choose', 'clip', 'compress', 'conj', 'conjugate', 'copy', 'ctypes', 'cumprod', 'cumsum', 'data', 'diagonal', 'dot', 'dtype', 'dump', 'dumps', 'fill', 'flags', 'flat', 'flatten', 'getfield', 'imag', 'item', 'itemset', 'itemsize', 'max', 'mean', 'min', 'nbytes', 'ndim', 'newbyteorder', 'nonzero', 'partition', 'prod', 'ptp', 'put', 'ravel', 'real', 'repeat', 'reshape', 'resize', 'round', 'searchsorted', 'setfield', 'setflags', 'shape', 'size', 'sort', 'squeeze', 'std', 'strides', 'sum', 'swapaxes', 'take', 'tobytes', 'tofile', 'tolist', 'tostring', 'trace', 'transpose', 'var', 'view']

Logical Indexing

In [11]:
# logical indexing
3 < x
Out[11]:
array([False, False, False, False,  True,  True,  True,  True,  True,  True], dtype=bool)
In [12]:
# negate logical indexing
~(3 < x)
Out[12]:
array([ True,  True,  True,  True, False, False, False, False, False, False], dtype=bool)
In [13]:
# logical indexing and 
(3 < x) & (x < 5)
Out[13]:
array([False, False, False, False,  True, False, False, False, False, False], dtype=bool)
In [14]:
# logical index or
(3 < x) | (x < 5)
Out[14]:
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)
In [15]:
# logical indexing and 
(3 < x) & (x < 5)
Out[15]:
array([False, False, False, False,  True, False, False, False, False, False], dtype=bool)
In [16]:
z = (3 < x) & (x < 5)
z
Out[16]:
array([False, False, False, False,  True, False, False, False, False, False], dtype=bool)
In [17]:
# if 0,1 is preferred, changethe datatype
z.astype(np.int)
Out[17]:
array([0, 0, 0, 0, 1, 0, 0, 0, 0, 0])

numpy arrays as a matrix

In [18]:
# createa 1x3 array
c = np.array([1,2,3])
print('c =',c,' of shape',c.shape)

c = [1 2 3]  of shape (3,)
In [19]:
#does not  convert to 3x1 since it is only a 1d array
# this is not what I would expected, so I will call it a nuance of arrays in python
c.transpose()
Out[19]:
array([1, 2, 3])
In [20]:
# we must reshape it
c = c.reshape((1,3))
print('c =',c,' of shape',c.shape)

c = [[1 2 3]]  of shape (1, 3)
In [21]:
c.transpose()
Out[21]:
array([[1],
       [2],
       [3]])
In [22]:
# or when we defined it, add brackets in it to create a true 1x3 array
c = np.array([[1,2,3]])
print('c =',c,' of shape',c.shape)

c = [[1 2 3]]  of shape (1, 3)
In [23]:
# now we can transpose it.
c.transpose()
Out[23]:
array([[1],
       [2],
       [3]])
In [24]:
print('c =',c,' of shape',c.shape)

c = [[1 2 3]]  of shape (1, 3)
In [25]:
c.shape
Out[25]:
$$\left ( 1, \quad 3\right )$$
In [26]:
# lets creat some matrices (but recall use the array)
a = np.array([[1, 2], [3, 4]])
b = np.array([[1, 4], [5, 7]])
In [27]:
a
Out[27]:
array([[1, 2],
       [3, 4]])
In [28]:
b
Out[28]:
array([[1, 4],
       [5, 7]])
In [29]:
a * b
Out[29]:
array([[ 1,  8],
       [15, 28]])
In [30]:
# to perform a matrix multiplecation on an array, use 
a @ b
Out[30]:
array([[11, 18],
       [23, 40]])
In [31]:
# or 
np.dot(a,b)
Out[31]:
array([[11, 18],
       [23, 40]])
In [32]:
# to solve ax=b, use 
linalg.solve(a,b)
Out[32]:
array([[ 3. , -1. ],
       [-1. ,  2.5]])
In [33]:
# or slower and less accurate
linalg.inv(a) @ b
Out[33]:
array([[ 3. , -1. ],
       [-1. ,  2.5]])
In [34]:
# or 
linalg.inv(a).dot(b)
Out[34]:
array([[ 3. , -1. ],
       [-1. ,  2.5]])

$ 3 x_0 + x_1 = 9 $

$ x_0 + 2x_1 = 8 $

$ Ax=b $

$ \left[\begin{matrix}3 & 1\\1 & 2\end{matrix}\right] \left[\begin{matrix}x_0\\x_1\end{matrix}\right] = \left[\begin{matrix}9\\8\end{matrix}\right] $

$ \left[\begin{matrix}x_0\\x_1\end{matrix}\right] = \left[\begin{matrix}2\\3\end{matrix}\right] $

In [35]:
# 3 * x0 + x1 = 9 and x0 + 2 * x1 = 8
a = np.array([[3,1], [1,2]])
b = np.array([[9],[8]])
x = np.linalg.solve(a, b)
x
Out[35]:
array([[ 2.],
       [ 3.]])
In [36]:
# lets create a 2x5 zeros matrix
c = np.zeros((2,5))
c
Out[36]:
array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])
In [37]:
# With numpy, we can create arrays of numbers as well as functions that accept arrays
f = lambda x,y: np.sqrt(x**2+y**2)
In [38]:
f(3,4)
Out[38]:
$$5.0$$
In [39]:
a = 3.1
b = np.arange(15,35, 0.75)
f(a,b)
Out[39]:
array([ 15.31698404,  16.05218054,  16.78868667,  17.52633732,
        18.26499384,  19.00453893,  19.74487275,  20.48590979,
        21.2275764 ,  21.96980883,  22.7125516 ,  23.45575622,
        24.19938016,  24.9433859 ,  25.68774027,  26.43241381,
        27.1773803 ,  27.92261628,  28.66810074,  29.41381478,
        30.15974138,  30.90586514,  31.65217212,  32.39864966,
        33.14528624,  33.89207134,  34.63899537])
In [40]:
# We can also perform numerical differentiation and integration
x = np.arange(-np.pi,np.pi, 0.01)
y = np.sin(x)*np.exp(-x)
plt.plot(x,y);
plt.xlabel('x')
plt.ylabel('$f(x)$')
plt.title('numerical integration example')
plt.show()
In [41]:
# Using the trapezoidal rule, we can numerically integrate the function
np.trapz(y,x)
Out[41]:
$$-11.5485471009$$

Symbolic Math with Sympy

in ipython use help(integrate) or integrate? to check syntax or learn about the command. To clear the previous modules use %reset to clear variables and imports. It is also very import to not import over other modules or difficult errors to troubleshoot may appear. This is especially important with numpy and sympy

In [42]:
%reset

Once deleted, variables cannot be recovered. Proceed (y/[n])? y

Now we can initialize the console to have very nice looking mathematical expressions. If we dont do this, the text will be simple ascii characters

In [43]:
import sympy as sp
import numpy as np
In [44]:
sp.init_printing()
In [45]:
x,y = sp.symbols('x y')

I really enjoy using sympy because I can type my equations in, evaluate them, and it shows them very nice, which is easier to read then lines of text. Simple expressions are easy to make

In [46]:
y = (sp.pi + x)**2
y
Out[46]:
$$\left(x + \pi\right)^{2}$$
In [47]:
y.subs(x,1.3)
Out[47]:
$$\left(1.3 + \pi\right)^{2}$$
In [48]:
# or to get a number
sp.N(y.subs(x,1.3))
Out[48]:
$$19.7277453004228$$
In [49]:
# calculus is also easy
sp.diff(y)
Out[49]:
$$2 x + 2 \pi$$
In [50]:
f = sp.Function('f')(x)
In [51]:
f = sp.sin(x)*sp.exp(-x)
f
Out[51]:
$$e^{- x} \sin{\left (x \right )}$$
In [52]:
F = sp.integrate(f)
F
Out[52]:
$$- \frac{e^{- x}}{2} \sin{\left (x \right )} - \frac{e^{- x}}{2} \cos{\left (x \right )}$$
In [53]:
F1 = sp.integrate(f,(x,-sp.pi,sp.pi ))
# as shown numerically, we can symbolically evaluate the function 
#  and get the same result
F1
Out[53]:
$$- \frac{e^{\pi}}{2} + \frac{1}{2 e^{\pi}}$$
In [54]:
float(F1)
Out[54]:
$$-11.548739357257748$$
In [55]:
sp.plot(f,(x,-sp.pi,sp.pi));
In [56]:
sp.plot(F,(x,-sp.pi,sp.pi));
In [57]:
# Other various math tasks can be much easier when using sympy
In [58]:
n = sp.symbols('n')
In [59]:
sum_fcn = sp.Sum(1/n**2, (n,1,10))
sum_fcn
Out[59]:
$$\sum_{n=1}^{10} \frac{1}{n^{2}}$$
In [60]:
# to evaluate the summation simply add the evalf method
sum_fcn.evalf()
Out[60]:
$$1.54976773116654$$
In [61]:
sp.latex(sum_fcn)
Out[61]:
'\\sum_{n=1}^{10} \\frac{1}{n^{2}}'

Since this is a windows machine, the latex command printed two backslashes, if we change that to single backslashes we can get the latex version of the summation $$ \sum_{n=1}^{10}\frac{1}{n^{2}} $$

In [62]:
from IPython.display import Latex
Latex('$sum\_fcn=' + sp.latex(sum_fcn) + '$' )
Out[62]:
$sum\_fcn=\sum_{n=1}^{10} \frac{1}{n^{2}}$
In [63]:
# linear algebra
In [64]:
s11, n22, e33 = sp.symbols("sigma11 nu22 epsilon33")
A = sp.Matrix([[s11],[n22],[e33]])
A
Out[64]:
$$\left[\begin{matrix}\sigma_{11}\\\nu_{22}\\\epsilon_{33}\end{matrix}\right]$$
In [65]:
m11, m12, m21, m22 = sp.symbols("m11, m12, m21, m22")
b1, b2 = sp.symbols("b1, b2")
In [66]:
A = sp.Matrix([[m11, m12],[m21, m22]])
A
Out[66]:
$$\left[\begin{matrix}m_{11} & m_{12}\\m_{21} & m_{22}\end{matrix}\right]$$
In [67]:
b = sp.Matrix([[b1], [b2]])
b
Out[67]:
$$\left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]$$
In [68]:
A**2
Out[68]:
$$\left[\begin{matrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\end{matrix}\right]$$
In [69]:
A.det()
Out[69]:
$$m_{11} m_{22} - m_{12} m_{21}$$
In [70]:
A.inv()
Out[70]:
$$\left[\begin{matrix}\frac{1}{m_{11}} + \frac{m_{12} m_{21}}{m_{11}^{2} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)} & - \frac{m_{12}}{m_{11} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)}\\- \frac{m_{21}}{m_{11} \left(m_{22} - \frac{m_{12} m_{21}}{m_{11}}\right)} & \frac{1}{m_{22} - \frac{m_{12} m_{21}}{m_{11}}}\end{matrix}\right]$$

mixed symbols and numerals linear algebra

In [71]:
x = sp.Symbol('x')   # x = var('x')
M = sp.Matrix([[2,x],[x,3]])
M
Out[71]:
$$\left[\begin{matrix}2 & x\\x & 3\end{matrix}\right]$$
In [72]:
M.eigenvals()
Out[72]:
$$\left \{ - \frac{1}{2} \sqrt{4 x^{2} + 1} + \frac{5}{2} : 1, \quad \frac{1}{2} \sqrt{4 x^{2} + 1} + \frac{5}{2} : 1\right \}$$
In [73]:
M.eigenvects()
Out[73]:
$$\left [ \left ( - \frac{1}{2} \sqrt{4 x^{2} + 1} + \frac{5}{2}, \quad 1, \quad \left [ \left[\begin{matrix}- \frac{x}{\frac{1}{2} \sqrt{4 x^{2} + 1} - \frac{1}{2}}\\1\end{matrix}\right]\right ]\right ), \quad \left ( \frac{1}{2} \sqrt{4 x^{2} + 1} + \frac{5}{2}, \quad 1, \quad \left [ \left[\begin{matrix}- \frac{x}{- \frac{1}{2} \sqrt{4 x^{2} + 1} - \frac{1}{2}}\\1\end{matrix}\right]\right ]\right )\right ]$$
In [74]:
M.eigenvects()[1][0]
Out[74]:
$$\frac{1}{2} \sqrt{4 x^{2} + 1} + \frac{5}{2}$$
In [75]:
Mval = M.eigenvects()[1][0]
Mval.evalf(subs={x:3.14})
Out[75]:
$$5.67955971794838$$
In [76]:
print(sp.latex(M))

\left[\begin{matrix}2 & x\\x & 3\end{matrix}\right]
In [77]:
from IPython.display import Latex
Latex('$M=' + sp.latex(M) + '$' )
Out[77]:
$M=\left[\begin{matrix}2 & x\\x & 3\end{matrix}\right]$
In [ ]:

Pythons list is a generic data storage object. it can be easily extended to a numpy array, which is specialized for numerical and scientific computation

In [78]:
np.zeros((5,3))
Out[78]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
In [79]:
np.array([[1,2],[3,4]])
Out[79]:
array([[1, 2],
       [3, 4]])
In [80]:
np.array([[1,2],[3,4]])
Out[80]:
array([[1, 2],
       [3, 4]])
In [81]:
# Matrix multiplication can be achieved using the dot method
i = [[1,0,0],[0,1,0],[0,0,1]]  # identiy matrix
a = [[4,3,1],[5,7,2],[2,2,2]]
np.dot(i,a)
Out[81]:
array([[4, 3, 1],
       [5, 7, 2],
       [2, 2, 2]])
In [82]:
#Or, matrix multiplication can be done if a matrix is explicitly defined
np.array(i) @ np.array(a)
Out[82]:
array([[4, 3, 1],
       [5, 7, 2],
       [2, 2, 2]])
In [83]:
# Notice, when arrays are mutliplied, we get the dot product 
m = np.array(i) * np.array(a)
m
Out[83]:
array([[4, 0, 0],
       [0, 7, 0],
       [0, 0, 2]])
In [84]:
m.T  # transpose
Out[84]:
array([[4, 0, 0],
       [0, 7, 0],
       [0, 0, 2]])
In [85]:
m**2
Out[85]:
array([[16,  0,  0],
       [ 0, 49,  0],
       [ 0,  0,  4]])
In [86]:
np.array(a)**2
Out[86]:
array([[16,  9,  1],
       [25, 49,  4],
       [ 4,  4,  4]])
In [87]:
m
Out[87]:
array([[4, 0, 0],
       [0, 7, 0],
       [0, 0, 2]])
In [88]:
m[:,2]
Out[88]:
array([0, 0, 2])
In [89]:
m[2,:]
Out[89]:
array([0, 0, 2])
In [90]:
m[:2,:2]
Out[90]:
array([[4, 0],
       [0, 7]])
In [91]:
m[1:,1:]
Out[91]:
array([[7, 0],
       [0, 2]])

Symbolic mathematics with sympy

In [92]:
# import sympy library and initialize latex printing
import sympy as sp
#sp.init_printing()
#sp.init_printing(use_latex='matplotlib')
sp.init_printing(use_latex='mathjax')
In [93]:
# add a symbolic character
x = sp.Symbol('x')
In [94]:
sp.sqrt(x**2)
Out[94]:
$$\sqrt{x^{2}}$$
In [95]:
r = sp.Rational(11, 13)
r
Out[95]:
$$\frac{11}{13}$$
In [96]:
float(r)
Out[96]:
$$0.8461538461538461$$
In [97]:
f = sp.Function('f')
f
Out[97]:
f
In [98]:
f(x)
Out[98]:
$$f{\left (x \right )}$$
In [99]:
h = sp.Lambda(x,x**2)
h
Out[99]:
$$\left( x \mapsto x^{2} \right)$$
In [100]:
w = 2*(x**2-x)-x*(x+1)
w
Out[100]:
$$2 x^{2} - x \left(x + 1\right) - 2 x$$
In [101]:
w.args
Out[101]:
$$\left ( - 2 x, \quad 2 x^{2}, \quad - x \left(x + 1\right)\right )$$
In [102]:
sp.simplify(w)
Out[102]:
$$x \left(x - 3\right)$$
In [103]:
sp.factor(x**2-1)
Out[103]:
$$\left(x - 1\right) \left(x + 1\right)$$
In [104]:
#partial fractions
y = 1/(x**2+3*x+2)
y
Out[104]:
$$\frac{1}{x^{2} + 3 x + 2}$$
In [105]:
sp.apart(y,x)
Out[105]:
$$- \frac{1}{x + 2} + \frac{1}{x + 1}$$
In [106]:
f = sp.Function('f')(x)
sp.diff(f,x)
Out[106]:
$$\frac{d}{d x} f{\left (x \right )}$$
In [107]:
y = sp.Symbol('y')
g = sp.Function('g')(x,y)
g.diff(x,y)
Out[107]:
$$\frac{\partial^{2}}{\partial x\partial y} g{\left (x,y \right )}$$
In [108]:
a,b,c,d = sp.symbols("a b c d")
M = sp.Matrix([[a,b],[c,d]])
M
Out[108]:
$$\left[\begin{matrix}a & b\\c & d\end{matrix}\right]$$
In [109]:
M*M
Out[109]:
$$\left[\begin{matrix}a^{2} + b c & a b + b d\\a c + c d & b c + d^{2}\end{matrix}\right]$$
In [136]:
# if ipython is to be used as a calculator initialize with 
from sympy import init_session
init_session() 

IPython console for SymPy 1.0 (Python 3.5.2-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/
In [111]:
from sympy import oo, Function, dsolve, Eq, Derivative, sin,cos,symbols
from sympy.abc import x
import sympy as sp
import numpy as np
import matplotlib.pyplot as mp
get_ipython().magic('matplotlib inline')
# this will print output as unicode
In [112]:
# assign a sympy variable
x = sp.var('x')
x
Out[112]:
$$x$$
In [113]:
#assign a function 
f =  sp.sin(6*x)*sp.exp(-x)
f
Out[113]:
$$e^{- x} \sin{\left (6 x \right )}$$
In [114]:
f.subs(x,3)
Out[114]:
$$\frac{1}{e^{3}} \sin{\left (18 \right )}$$
In [115]:
float(f.subs(x,3))
Out[115]:
$$-0.037389453398415345$$
In [116]:
sp.plot(f,(x,-5,5))
Out[116]:
<sympy.plotting.plot.Plot at 0xada2b00>
In [117]:
# a onetime pretty print
sp.pprint(f)

 -x         
ℯ  ⋅sin(6⋅x)
In [118]:
#or we can print the latex rendering
sp.latex(f)
Out[118]:
'e^{- x} \\sin{\\left (6 x \\right )}'
In [119]:
# first derivative
df = f.diff()
df
Out[119]:
$$- e^{- x} \sin{\left (6 x \right )} + 6 e^{- x} \cos{\left (6 x \right )}$$
In [120]:
# differentaite f'' wrt x
sp.diff(f,x,1)
Out[120]:
$$- e^{- x} \sin{\left (6 x \right )} + 6 e^{- x} \cos{\left (6 x \right )}$$
In [121]:
# substitute x with pi
f.subs(x,np.pi)
Out[121]:
$$-3.17530720082064 \cdot 10^{-17}$$
In [122]:
#%% Numeric Computation from the documentation
from sympy.abc import x
In [123]:
# lambdify using the math module, 10^2 faster than subs
expr = sp.sin(x)/x
f = sp.lambdify(x,expr)
f(3.14)
Out[123]:
$$0.0005072143046136395$$
In [124]:
# lambdify using numpy
f = sp.lambdify(x,expr, "numpy")
f(np.linspace(1,3.14,20))
Out[124]:
array([  8.41470985e-01,   8.06076119e-01,   7.67912588e-01,
         7.27262596e-01,   6.84424864e-01,   6.39711977e-01,
         5.93447624e-01,   5.45963742e-01,   4.97597617e-01,
         4.48688937e-01,   3.99576866e-01,   3.50597122e-01,
         3.02079129e-01,   2.54343238e-01,   2.07698064e-01,
         1.62437944e-01,   1.18840569e-01,   7.71647744e-02,
         3.76485431e-02,   5.07214305e-04])

Signal Processing

In [125]:
## Signal Processing
#Page 174 Introduction for python for Science - David Pine

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline') # inline plotting
width = 2.0
freq = 0.5
t = np.linspace(-10, 10, 101) # linearly space time array
g = np.exp(-np.abs(t)/width)*np.sin(2.0 * np.pi * freq * t)
dt = t[1]-t[0] # increment between times in time array
G = fftpack.fft(g) # FFT of g
f = fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]
f = fftpack.fftshift(f) # shift frequencies from min to max
G = fftpack.fftshift(G) # shift G order to coorespond to f
fig = plt.figure(1, figsize=(8,6), frameon=False)
ax1 = fig.add_subplot(211)
ax1.plot(t, g)
ax1.set_xlabel('t')
ax1.set_ylabel('g(t)')
ax2 = fig.add_subplot(212)
ax2.plot(f, np.real(G), color='dodgerblue', label='real part')
ax2.plot(f, np.imag(G), color='coral', label='imaginary part')
ax2.legend()
ax2.set_xlabel('f')
ax2.set_ylabel('G(f)')
plt.show()
In [126]:
from pylab import *
from scipy import fft
get_ipython().magic('matplotlib inline') # inline plotting
N = 2**9
F = 25
t = arange(N)/float(N)
x = cos(2*pi*t*F) + rand(len(t))*3
subplot(2,1,1)
plot(t,x)
ylabel('x []')
xlabel('t [seconds]')
title('A cosine wave')
grid()

subplot(2,1,2)
f = t*N
xf = fft(x)
plot(f,abs(xf))
title('Fourier transform of a cosine wave')
xlabel('xf []')
ylabel('xf []')
xlim([0,N])
grid()
show()

# note the spike at 25 hz and 512-25

Differential Equations

In [127]:
# see here for example scripts
# C:\Users\Neal\Anaconda3\Lib\site-packages\sympy\mpmath\tests
from sympy import Function, dsolve, Eq, Derivative, sin,cos,symbols
from sympy.abc import x
import numpy as np
import sympy as sp
import matplotlib.pyplot as mp
import matplotlib.pyplot as plt
In [128]:
f = Function('f')
deq = dsolve(Derivative(f(x), x,x) + 9*f(x), f(x))
deq
Out[128]:
$$f{\left (x \right )} = C_{1} \sin{\left (3 x \right )} + C_{2} \cos{\left (3 x \right )}$$
In [129]:
diffeq1_sym = deq.args[1]
diffeq1_sym
Out[129]:
$$C_{1} \sin{\left (3 x \right )} + C_{2} \cos{\left (3 x \right )}$$
In [130]:
diffeq1 = diffeq1_sym.subs({'C1':1, 'C2':0.5})
diffeq1
Out[130]:
$$\sin{\left (3 x \right )} + 0.5 \cos{\left (3 x \right )}$$
In [131]:
diffeq1_f = sp.lambdify(x,diffeq1, "numpy")
diffeq1_f
Out[131]:
<function numpy.<lambda>>
In [132]:
diffeq1arr = diffeq1_f(np.linspace(1,3.14,20))
diffeq1arr
Out[132]:
array([-0.35387624, -0.68544104, -0.93948885, -1.08728921, -1.11212728,
       -1.0111941 , -0.79590429, -0.49060511, -0.12982305,  0.24564078,
        0.59332492,  0.87390954,  1.05566313,  1.11803104,  1.05396004,
        0.87069598,  0.5889643 ,  0.24062625, -0.1349244 , -0.49521635])
In [133]:
# numpy plot
plt.plot(diffeq1_f(np.linspace(-10,10,2000)));
plt.title('plot of the numpy array');
In [134]:
# sympy plot
sp.plot(diffeq1, title='plot of the sympy function');
In [ ]:

Differential Equations Quiver plot

In [135]:
# quiver plot
f=lambda x: [x[0]**2 - 2*x[0] - x[1] + 0.5, x[0]**2 + 4*x[1]**2 - 4]
x,y=np.mgrid[-0.5:2.5:24j,-0.5:2.5:24j]
U,V=f([x,y])
plt.quiver(x,y,U,V,color='r', \
         linewidths=(0.2,), edgecolors=('k'), \
         headaxislength=5)
plt.show()