CS计算机代考程序代写 flex python algorithm Computational Methods

Computational Methods
for Physicists and Materials Engineers
2
Basic programming II NumPy, matplotlib, SymPy

From Lecture 0: Algorithm
What is algorithm?
“An algorithm is a finite sequence of well‐defined, computer‐ implementable instructions.” (Wikipidea)
Computers only know two states: “0” and “1”. Naturally, computers can only deal with discrete numbers, rather than continuous functions. So, a continuous function must be represented by an array of discrete numbers. E.g. force field on a line f(s) is continuous; but in a computer, it must be stored in a discrete form:
f(s) f(s)  f(s) 11121N
f(s) f(s)  f(s) 21222N
f3(s1) f3(s2)  f3(sN)
In math, the array of numbers is vector (1D), matrix (2D), or higher‐ order tensor (nD) [think of the force field on a bulk material]
The Python package that can handle the arrays is NumPy

Installation
Import
https://numpy.org/
> pip install numpy
# at the head of a code
import numpy as np

Creating arrays
Dimension is the number of “[…]” 2D arrays are matrices, e.g. d, e, f
NumPy
numpy.array()
a = np.array(42)
b = np.array([42])
c = np.array([1, 2, 3])
d = np.array([[1, 2, 3]])
e = np.array([[1], [2], [3]])
f = np.array([[1, 2, 3], [4, 5, 6]])
g = np.array([[[1, 2, 3], [4, 5, 6]], [[1, 2, 3], [4, 5, 6]]])
for x in [a, b, c, d, e, f, g]:
print(f”This array is \n{x}”)
print(f”Its type is {type(x)}; its dimension is {x.ndim}.”)
1 1 2 3
d1 2 3 , e2 , f 
13  456 3 23
 31

Check the type of data in a NumPy array
arr = np.array([[1, 2, 3], [4, 5, 6]])
print(arr.dtype)
Define a NumPy array with specified data type
arr = np.array([[1, 2, 3], [4, 5, 6]], dtype=’f’) # i: integer, f: float, b: boolean, S: string print(arr.dtype)
Change the data type of an existing array
arr = np.array([[1, 2, 3], [4, 5, 6]], dtype=’i’) print(arr.dtype)
arr = arr.astype(‘f’)
print(arr.dtype)
NumPy
.dtype .astype()

Creating arrays
a = np.empty((2, 3))
b = np.zeros((3, 4))
c = np.ones((4, 5))
# the argument is a tuple; don’t forget “(…)”
for x in [a, b, c]:
print(f”Matrix is \n{x}; \n\
shape is {x.shape}; size is {x.size}.”)
d = np.array([1, 2, 3], [4, 5, 6])
e = np.zeros(d.shape); print(e)
NumPy
numpy.empty(), numpy.zeros(), numpy.ones()

Creating 1D arrays
import numpy as np
import matplotlib.pyplot as plt
NumPy
numpy.arange(), numpy.linspace()
# np.arange(start, stop, step); numbers in range [start, stop)
f = np.arange(3, 10, 2); print(f); print(f”dimension is {f.ndim}”) g = np.arange(3, 11, 2); print(g)
h = np.arange(3, 12, 2); print(h)
# arguments can be float;
# but it’s better to use np.linspace for this case
# np.linspace(start, stop, number of points, if include end point)
num_points = 40
x = np.linspace(0, 2*np.pi, num_points, endpoint=True) y = np.sin(x)
plt.plot(x, y, ‘o’)
plt.show()

Creating random arrays
NumPy
numpy.random.rand(), numpy.random.randint(), numpy.random.normal()
# generate an array of random floats from uniform distribution over [0, 1)
print(np.random.rand(3, 2))
# generate an array of random integers from uniform distribution over [low, high)
print(np.random.randint(3, 9, size=(3, 2)))
# generate an array of random floats from normal distribution with specified mean, standard deviation and number of points mean = 1; stddev = 2; num_points = 10000000
y = np.random.normal(mean, stddev, num_points)
plt.hist(y, bins=np.arange(-6, 8, 0.1), density=True) plt.show()

Indexing arrays
NumPy
a = np.array(42)
b = np.array([42])
c = np.array([1, 2, 3])
d = np.array([[1, 2, 3]])
e = np.array([[1], [2], [3]])
f = np.array([[1, 2, 3], [4, 5, 6]])
g = np.array([[[1, 2, 3], [4, 5, 6]], [[1, 2, 3], [4, 5, 6]]])
print(“a: “); print(a)
print(“b: “); print(b[0])
print(“c: “); print(c[0:-1])
print(“d: “); print(d[0]); print(d[0, 1]); print(d[:, 1]) print(“e: “); print(e[1]); print(e[1, 0]); print(e[1, :]) print(“f: “); print(f[1, 2]); print(f[0:2, 0:-1]) print(“g: “); print(g[0, 1, 2])
[]

Searching arrays
arr = np.array([[10001, -12352, 234, 32], [23, 1, 104, -4],
print(arr)
# set the limits
[-23, -1003, 24, -98]])
lower = -100
upper = 100
# set all elements which > upper to upper
arr[np.where(arr > upper)] = upper
# set all elements which < lower to lower arr[np.where(arr < lower)] = lower print(arr) Sorting arrays arr = np.random.randint(1, 10, 10) print(np.sort(arr)) NumPy numpy.where() # explicit element-wise operations for i, _ in enumerate(a): sum1[i] = a[i] + b[i] prod1[i] = a[i] * b[i] power1[i] = a[i]**2 compare1[i] = a[i] == b[i] NumPy Element‐wise operations a = np.arange(1, 11) b = np.arange(2, 12) print(a); print(b) sum1 = np.copy(a); prod1 = np.copy(a); power1 = np.copy(a) compare1 = np.empty(a.shape, dtype=bool) # for NumPy arrays, "+", "-", "*", "/", "**", and comparisons are all element-wise operations sum2 = a + b prod2 = a * b power2 = a**2 compare2 = a == b print(sum1); print(sum2) print(prod1); print(prod2) print(power1); print(power2) print(compare1); print(compare2) Element‐wise operations More complicated NumPy operations arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) # explicit element-wise operations f1 = np.empty(arr.shape) for i, x in np.ndenumerate(arr): f1[i] = (math.sin(math.sqrt(abs(x))) - np.pi)**3 # NumPy operations are element-wise f2 = (np.sin(np.sqrt(np.abs(arr))) - np.pi)**3 print(f1) print(f2) NumPy Element‐wise operations Advantage of NumPy element‐wise operations: much faster import numpy as np import matplotlib.pyplot as plt import math, time def numpy_elemwise_comparison(arr_size): arr = np.linspace(1, 100, arr_size) time0 = time.time_ns() # explicit element-wise operations f1 = np.empty(arr.shape) # NumPy operations are element-wise f2 = (np.sin(np.sqrt(np.abs(arr))) - np.pi)**3 time2 = time.time_ns() return (time2 - time1) / (time1 - time0) arr_size_list = np.arange(10000, 1000001, 10000) time_lapse_list = [] for arr_size in arr_size_list: NumPy for i, x in np.ndenumerate(arr): f1[i] = (math.sin(math.sqrt(abs(x))) - np.pi)**3 time1 = time.time_ns() time_lapse_list.append(numpy_elemwise_comparison(arr_size)) NumPy Element‐wise operations NumPy element‐wise operation is about 16 times faster than explicit element‐wise operation (why?) Physical quantities coordinate of a particle r = (x1, x2, x3) velocity of a particle v = (v1, v2, v3) force f = (f1, f2, f3) electric field E = (E1, E2, E3) O e1 e2 mass flux/electric current J = (J1, J2, J3) Gradients of scalar fields concentration is a scalar field c(r) r(x1, x2), g(x1, x2), b(x1, x2) Its gradient is c x1  c(r)c x  Mass flows from high‐ to low‐c region Fick's1stlaw: J(r)Dc(r) x1 Recap: linear algebra e3 Vectors in 3D space 2 c x  x 32 Geometry features Material properties Recap: linear algebra normalofasurfacen=(n ,n ,n )andn2 n2 n2 1 123123 tangentofacurvel=(l ,l ,l )andl2 l2 l2 1 123123 Vectors in 3D space pyroelectricity: temperature change induces the change in the state of electric polarization Electric polarization P is a vector; temperature change ΔT is a scalar P  ΠΔT P Π ‒ ‒ ‒ ‒ p = qd ++++ 11 P P Π ΔT ‒‒‒‒ ++++ 22 PΠ 33 Vectors in nD space A function can be represented by a vector ψ(x) ψ(x) Δx ψ(x1) When Δx → 0 and N →, the vector is exactly same as the function (Hilbert space) ψ(x ) ψ(x) 2   Δxψ (ket) ψ(x ) N Why Δx? Recap: linear algebra x xxxxxxxxxxx 1 2 3 4 5 6 7 8 9 10 Vectors in nD space A function can be represented by a vector In quantum mechanics, ψ(x) is a wave function ψ(x1) Recap: linear algebra ψ(x ) 2 ket:ψΔx,bra:ψψ(x1)ψ(x2)ψ(xN) Δx ψ(x ) N Inner product of a wave function  ψψψ(x1) ψ(x2)  ψ(xN)  Δx ψ(x)ψ(x)Δx 2 ii Δx0   ψ (x)ψ(x)dx ψ(x1) ψ(x ) N ψ(x ) N  i1 Physical quantities stress Recap: linear algebra Matrices in 3D space σ σ σFA FA FA 11 12 131 1 1 2 1 3 σσ σ σ F A F A F A 21 22 232 1 2 2 2 3 σσσFAFAFA 31 32 333 1 3 2 3 3 εεε112113 strain εε ε ε  ΔL ΔL  ΔL L ΔL ΔL   ΔL L ΔL ΔL  ΔL ΔL  21 31 11 12 13  1 1 2 L L 2 L L  112 123 2122232LL 222LL 21 32 113123   ε31 ε32 ε33ΔLΔLΔLΔL ΔLL  2LL2LL33 31 32 Gradients of vector fields u (x ,x ) 112 u2(x1 ,x2) x2  x 1 Recap: linear algebra Matrices in 3D space displacement field u(r) is a vector field; its gradient is u1 x1 u1 x2 u1 x3  1T uu2 x1 u2 x2 u2 x3 , ε2uu  ux ux ux 313233 u(x1 , x2 ) Geometry features Recap: linear algebra Matrices in 3D space At a point of a surface, establish the coordinate system with e3 parallel to the normal n and e1 and e2 parallel to the tangential directions. The curvature of this point is n n 11 κnx1 x2  n n 22 x1 x2  e3 e1 e2 Matrices in 3D space A crystal may be polarized by an electric field Material properties electric polarizability P ‒‒‒‒ ++++ ‒‒‒‒ ++++ ‒‒‒‒ + + + +p P ‒‒‒‒ ++++p PαE or P α  1   11 α α E 12 13  1  Recap: linear algebra Pα  2   21 22 23  2  α α E PαααE  3   31 32 33  3  E Material properties electrical conductivity Recap: linear algebra Matrices in 3D space electric current density can be induced by electric field J1  σ11 σ12 σ13 E1  JσE or Jσ σ σ E  2   21 22 J  σ σ 23  2  σ E  33  3   3   31 32 σ is electrical conductivity tensor This is Ohm’s law. Conduction in a wire is of 1D, J = σE. J is current per area: J = I/A. Voltage is the potential difference: V = EL. So, J  σE  I  σ V  V  I L  IR A L σA A L x1 x2 2H0, 2H0, 2H 0 Recap: linear algebra Matrices in nD space Hessian matrix: 2nd‐order partial derivatives of scalar fields Landscape H(x1, x2); what is the local curvature? local H minimum x2 x2 x x 1212 x1 x2 2H0, 2H0, 2H 0 Recap: linear algebra Matrices in nD space Hessian matrix: 2nd‐order partial derivatives of scalar fields Landscape H(x1, x2); what is the local curvature? local H maximum x2 x2 x x 1212 x1 x2 2H0, 2H0, 2H 0 Recap: linear algebra Matrices in nD space Hessian matrix: 2nd‐order partial derivatives of scalar fields Landscape H(x1, x2); what is the local curvature? H x2 x2 x x 1212 saddle point Hessianmatrix:  1  2H Recap: linear algebra Matrices in nD space Hessian matrix: 2nd‐order partial derivatives of scalar fields Landscape H(x1, x2); what is the local curvature?  2H 2H  xx 2H x2 1 2  2H  x2 12 2 x2  At a critical point (x1, x2) , rotate the coordinate system (e1, e2) such that, under the new system, 2H x2 0   x x 2H  1  2H  x2 2 0 x2 Then, from the signs of components, we can know if the point is local minimum, local maximum or saddle point Hessian matrix: 2H x   x x x2 x x  Recap: linear algebra Matrices in nD space Hessian matrix: 2nd‐order partial derivatives of scalar fields A scalar function can depend on multiple variables H(x1, x2, ∙∙∙, xN) E.g. potential energy of a material is a function of the coordinates of all atoms E(r1, r2, ∙∙∙, rN)  2H 2H 2H   x2 xx  xx  1121N  2H 2H 2H   212 2 2N   2H 2H  2H  x x x x x2  1N 2N N Recap: linear algebra Matrices in nD space Jacobian matrix: 1st‐order partial derivatives of vector fields A vector function can depend on multiple variables f f(x,x,,x) 1112 N 2212N f f  f (x ,x ,,x )f(x)  f  f (x ,x ,,x ) MM12 N Gradient of f is Jacobian matrix: f f f  f  f ff 111 x1 x2 xN  222 Jacobianmatrix: xx1 x2 xN   ff f MMM xx x 12N Matrices in nD space An operator turns a function to another function An operator can be represented by a matrix First derivative d dx From Lecture 0, dψ dx  ψ(xi1)ψ(xi1) i 2Δx ψ(x) dψ dx ψ(x2)ψ(x1) 1 Δx ψ(xN)ψ(xN1) dxN Δx dψ x x x x x x x x x x x Recap: linear algebra 1 2 3 4 5 6 7 8 9 10 First derivative  dψ   11  ψ(x1)  ψ(x) dx1Δx Δx  dψ   2  ψ(x3)   1 0 1  dx 2   2Δx 2Δx    ψ(xi1) ψ(x)  dψ 101 2Δx 2Δx i ψ(xi1) dx i      1 0 1  dψ 2Δx ψ(xN2) ψ(x )  dxN1 2Δx  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra Matrices in nD space First derivative dx1Δx Δx  dψ   2  ψ(x3)   1 0 1  dx 2   2Δx 2Δx    ψ(xi1) ψ(x)  dψ 101 2Δx 2Δx i ψ(xi1) dx i      1 0 1  dψ 2Δx ψ(xN2) ψ(x )  dxN1 2Δx  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra Matrices in nD space  dψ  dψ ψ(x2)ψ(x1) 11 dx1 Δx  ψ(x1)  ψ(x) First derivative dx1Δx Δx  dψ   2  ψ(x3)   1 0 1  dx 2   2Δx 2Δx    ψ(xi1) ψ(x)  dψ 101 2Δx 2Δx i ψ(xi1) dx i      1 0 1  dψ 2Δx ψ(xN2) ψ(x )  dxN1 2Δx  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra Matrices in nD space  dψ  dψ ψ(x3)ψ(x1)  1 1 dx2 2Δx ψ(x1) ψ(x) First derivative dx1Δx Δx  dψ   2  ψ(x3)   1 0 1  dx 2   2Δx 2Δx    ψ(xi1) ψ(x)  dψ 101 2Δx 2Δx i ψ(xi1) dx i      1 0 1  dψ 2Δx ψ(xN2) ψ(x )  dxN1 2Δx  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra Matrices in nD space  dψ  dψ  ψ(xi1)ψ(xi1) 11 dxi 2Δx  ψ(x1)  ψ(x) First derivative Matrices in nD space dψ ψ(xN)ψ(xN2)  dψ   ψ(x1)  ψ(x)   1 1 dxN1 2Δx dx1Δx Δx 2  ψ(x3)   dψ    1 0 1  dx 2   2Δx 2Δx     ψ(xi1) ψ(x) dψ 101 dx 2Δx 2Δx i ψ(xi1) i   dψ    1 0 1   dxN1 2Δx 2Δx ψ(xN2) ψ(x )  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra First derivative dx1Δx Δx  dψ   2  ψ(x3)   1 0 1  dx 2   2Δx 2Δx    ψ(xi1) ψ(x)  dψ 101 2Δx 2Δx i ψ(xi1) dx i      1 0 1  dψ 2Δx ψ(xN2) ψ(x )  dxN1 2Δx  dψ 1 1 N1 dxN   Δx Δx  ψ(xN)  Recap: linear algebra Matrices in nD space  dψ  dψ ψ(xN)ψ(xN1) 11dxN Δx  ψ(x1)  ψ(x) First derivative Matrices in nD space d dψ(x) dx dx ψ(x) ψ(x1) dψ  11 ψ(x)   2  dx1 Δx Δx  dψ      1 dx2 2Δx 0 1 2Δx   ψ(x3)   Recap: linear algebra    ψ(xi1) 1 0 1 ψ(x) dψ 2Δx2Δx i   ψ(xi1) dx i  11  dψ   ψ(xN2) 2Δx    2Δx 0 dxN1  1  ψ(xN1)   dψ   1 Δx Δx N  dxN    ψ(x )  dx dx ψ(x) dψ dx ˆ ψ Recap: linear algebra Matrices in nD space dψ(x) d First derivative A Matrices in nD space An operator can be represented by a matrix Second derivative 2 d From Lecture 0, dx2 d2ψ  ψ(xi1)2ψ(xi )ψ(xi1) What is this matrix? (Homework) Schrödinger equation for a particle in 1D space: Inmatrixform H ψEψ  ij j i Recap: linear algebra 22 dx i (Δx) ˆ Hψ(x)   2 d2 V(x)ψ(x)  Eψ(x) 2  2mdx  Eigenvalue problem: there exist infinite E values s.t. the equation holds. These E’s are eigenvalues; they are possible energies Creation a 1 e3 aa2   a a  M 1 a = np.array([-1, 1, 1]) b = np.array([2, -2, 0]) ‐2 O ‐1 2 1 e2 Magnitude or l2 norm b2 e1 M aa a2 2 i i1 a_norm = np.linalg.norm(a) NumPy Vectors represented by 1D arrays There are other norms; but they do not correspond to the magnitude of a vector in 3D Vectors represented by 1D arrays Multiplication with a scalar  λa  e3 a 1 λaλa2   a λa  M |a| ‐2 O print(a / np.linalg.norm(a)) # create unit vector of a print(0.5 * b) 0.5b e2 NumPy b2 e1 Inner product (dot product) M abab e3 a Cross product (vector product) Only in 3D space O e2 ii i1 print(np.inner(a, b)) e e eabab 1 2 323 32 abdeta a a ab ab  1 2 331 13 bbb abab 1 2 312 21 print(np.cross(a, b)) NumPy Vectors represented by 1D arrays bec 1 Creation AAA  11 12 1N  A = np.array([[-1, 1, 1], [2,-2, 0]]) B = np.array([[0.5, 1, 0], [ 0,0.5, 1], Transpose [ 1, 0, 0]]) #Bisasquarematrix AAA print(np.transpose(A)) NumPy Matrices represented by 2D arrays AAAA 21 22 2N  AAA  M1 M2 MN MN #Aisa2x3matrix  11 21 M1  T 12 22 M2 AAAA  AAA  1N M2 NM NM Matrices represented by 2D arrays Trace (only for square matrices) print(np.trace(B)) NumPy trA  A Norm We may simply generalize the defini􏱤on of magnitude of a vector → Frobenius norm M MN aaa2A A2 A_norm = np.linalg.norm(A, 'fro') B_norm = np.linalg.norm(B, 'fro') M ii i1  2 i1 F i1 k1 ik i The norm means some kind of “magnitude” of a matrix. There are other definitions of norms Matrices represented by 2D arrays Multiplication with a scalar NumPy λA λA  λA   11 12 1N λAλA λA  λA  21 22 2N print(A / np.linalg.norm(A, 'fro')) print(np.trace(B) * B)  λA λA  λA   M1 M2 MNMN Matrix product NumPy Matrices represented by 2D arrays AAABBB  11 12 1N   11 12 1K  21 22 2N 21 22 2K AA A A ,BB B B  AAA BBB  M1 M2 MN MN  N1 N2 NK NK NNN AB ABAB   i1 i1 i1  1i i1  1i i2 1i iK NNN ABAB ABAB   i1 i1 i1    2i i1 2i i2 2i iK N N N  ABABAB     i1 i1 i1 MK Mi i1 Mi i2 Mi iK Matrix product Matrices represented by 2D arrays A = np.array([[-1, 1, 1], [ 2, -2, 0]]) B = np.array([[0.5, 1, 0], [0,0.5, 1], [ 1, 0, 0]]) M, N_A = A.shape; N_B, K = B.shape # first check if A can be right multiplied by B if N_A != N_B: # raise error and stop the code NumPy raise Exception("Dimensions of two matrices do not match!") print(np.matmul(A, B)) Inner product of 2 vectors is equivalent to the matrix product of a row vector (1M matrix) and a column vector (M1 matrix) a = np.array([[-1, 1, 1]]) b = np.array([[ 2], # a row vector [ 0]]) c = np.array([[2, -2, 0]]) # a column vector # another row vector [-2], print(np.matmul(a, b)) print(np.matmul(a, np.transpose(c))) Matrices represented by 2D arrays Element‐wise multiplication Two matrices must have the same shape A = np.array([[-1, 1, 1], [ 2, -2, 0]]) B = np.array([[0.5, 1, 0], [ 0,0.5, 1]]) # first check if A and B have the same shape if A.shape != B.shape: # raise error and stop the code NumPy raise Exception("Dimensions of two matrices do not match!") print(A * B) Matrices represented by 2D arrays Recall the inner product of 2 vectors (of same shape) Frobenius inner product of 2 matrices (of same shape) A = np.array([[-1, 1, 1], [ 2, -2, 0]]) B = np.array([[0.5, 1, 0], [ 0,0.5, 1]]) ik ik  ki ik # first check if A and B have the same shape if A.shape != B.shape: # raise error and stop the code NumPy M abab raise Exception("Dimensions of two matrices do not match!") print((A * B).sum()) print(np.trace(np.matmul(np.transpose(A), B))) ii i1 MN NM N A:BAB  ATB ATBtrATB  i1 k1 k1i1  k1 kk Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix AAA  11 12 13  detAdetA A A   21 22 23  AAA  31 32 33  NumPy 11 12 11 22 12 21 AA Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix AAA  11 12 13  detAdetA A A   21 22 23  AAA  31 32 33  AAAAAAAAA 11 22 33 12 23 31 13 21 32 NumPy 11 12 11 22 12 21 AA Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix AAA  11 12 13  detAdetA A A   21 22 23  AAA  31 32 33  AAAAAAAAA 11 22 33 12 23 31 13 21 32 A A A A A A A A A 13 22 31 23 32 11 33 12 21 NumPy 11 12 11 22 12 21 AA Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix AAA 11 12 13 AA 11 22 23   detAdet A A A (1) A det NumPy 11 12 11 22 12 21 AA 21 22 23 11 A A AAA  31 32 33   32 33  Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix 21 22 23 AAA 11 A A  31 32 33   32 33  12 21 (1) A det 12 A A AA  31 33  NumPy 11 12 11 22 12 21 AA AAA 11 12 13 AA 11 22 23   detAdet A A A (1) A det 23 Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix   detAdet A A A (1) A det NumPy 11 12 11 22 12 21 AA AAA 11 12 13 AA 11 22 23 21 22 23 11 A A  32 33  AA AA AAA  31 32 33  12 21 23 13 21 22 (1) A det  (1) A det 12 AA 13 A A  31 33   31 32  Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix detAdetA A A A A A  21 22  Determinant of a 33 matrix AAA  11 12 13  NumPy 11 12 11 22 12 21 AA detAdetA A A A (A A A A ) 21 22 23 11 22 33 23 32 AAA  31 32 33  A (A A A A )A (A A A A ) 12 21 33 23 31 13 21 32 22 31 Determinant of a NN matrix? Matrices represented by 2D arrays Why define determinant in this way? Geometric meaning Area of a parallelogram in 2D a  b  1 1 e2 e1 e2 e3 0  b Areaab, aa, bb 2 2 0 0 abdeta a 0  0  12 b b 0abab a 1 2 12 21 Area ab a b a b e1 Extension 12 21 e3 to 3D deta bdeta babba 111212 ab 2 2 Determinant is area of the parallelogram in 2D NumPy Matrices represented by 2D arrays Why define determinant in this way? Geometric meaning Volume of a parallelepiped in 3D Volume  Areaheight  ab c ab   cab  ab  c  a b a b  e3 c 123 32 cabc ab ab  231 13 c a b a b 31 2 c (a b a b ) 2 1 Proj.ofc onto normal of ab plane b 12332  c (a b  a b ) 21331 c (a b a b ) O e1 e2 a 31221 a b c 111 deta b cdeta b c 222 abc 332 Determinant is volume of the parallelepiped in 3D NumPy Matrices represented by 2D arrays Determinant (only for square matrices) Determinant of a 22 matrix A = np.array([[-1, 1], [ 2, -3]]) print(np.linalg.det(A)) Determinant of a 33 matrix B = np.array([[-1, 1, 0], [2,-3, 2], [ 0, 0, 1]]) print(np.linalg.det(B)) Determinant of a 55 matrix C = np.array([[-1, 1, 0, 0, 0], print(np.linalg.det(C)) [2,-3, 2, 0, 0], [ 0, 0, 1, 0, 0], [ 0, 0, 0, 1, 0], [0, 0, 0, 0, 2]]) NumPy Matrices represented by 2D arrays Inverse of a square matrix The inverse of a square matrix A is a square matrix B such that AB = I We denote B as A‒1 A1 = np.array([[-1, 1], [ 2, -3]]) A1inv = np.linalg.inv(A) print(A1inv) print(np.matmul(A1, A1inv)) Some matrices are not invertible try: A2 = np.array([[-1, 1], [ 2, -2]]) try to run this block try: what if error occurs except: A2inv = np.linalg.inv(A2) print(A2inv) print("This matrix is not invertible!") NumPy except: Import https://matplotlib.org/ Installation > pip install matplotlib
# at the head of a code
import matplotlib.pyplot as plt

import numpy as np
# prepare four trigonometric functions
matplotlib
Template
Although there are many scripts for plotting figures by matplotlib, here
we provide one of the templates
step = 0.2
x = np.arange(0, 2*np.pi + step, step)
sinx = np.sin(x)
cosx = np.cos(x)
tanx = np.tan(x)
cotx = 1/np.tan(x)
func = (sinx, cosx, tanx, cotx)
title = (r'(a) $\sin(x)$’, r'(b) $\cos(x)$’, r'(c) $\tan(x)$’, r'(d) $\cot(x)$’)

import itertools
fig, axs = plt.subplots(2, 2)
matplotlib
for n, index in enumerate(itertools.product(range(2), repeat=2)): ax = axs[index]
ax.plot(x, func[n], ‘o-‘, linewidth=2, color=’#0000FF’)
ax.set_xlim([0, 2*np.pi])
ax.set_ylim([-1.2, 1.2])
ax.set_xticks([0, 0.5*np.pi, np.pi, 1.5*np.pi, 2*np.pi]) ax.set_xticklabels([r’$0$’, r’$\frac{\pi}{2}$’, r’$\pi$’,
r’$\frac{3\pi}{2}$’, r’$2\pi$’]) ax.set_yticks([-1, 0, 1]) ax.set_yticklabels([r’$-1$’, r’$0$’, r’$1$’]) ax.tick_params(axis=’x’, labelsize=10) ax.tick_params(axis=’y’, labelsize=10) ax.set_xlabel(r’$x$ (rad)’, fontsize=12) ax.set_ylabel(r’$y$’, fontsize=12) ax.set_title(title[n], fontsize=12) ax.grid(True)
fig.set_size_inches(6,5)
fig.suptitle(rf’Trigonometric functions ($\Delta x =$ {step})’) plt.tight_layout()
plt.show()

fig, axs = plt.subplots(2, 2)
matplotlib
for n, index in enumerate(itertools.product(range(2), repeat=2)): ax = axs[index]
axs[0,0] axs[0,1]
axs[1,0] axs[1,1]

ax.plot(x, func[n], ‘o-‘, linewidth=2, color=’#0000FF’)
matplotlib
circles and solid line hexadecimal color code
https://htmlcolorcodes.com/

ax.set_xlim([0, 2*np.pi])
ax.set_ylim([-1.2, 1.2])
y  [‐1.2, 1.2]
x  [0, 2π]
matplotlib

xticks
r’$\frac{\pi}{2}$’
raw string:
LaTeX syntax
math mode:
fraction a/b:
Greek letter:
r’…’
$…$
\frac{a}{b}
\pi
LaTeX tutorial https://gking.harvard.edu/files/lshort2.pdf
matplotlib
ax.set_xticks([0, 0.5*np.pi, np.pi, 1.5*np.pi, 2*np.pi]) ax.set_xticklabels([r’$0$’, r’$\frac{\pi}{2}$’, r’$\pi$’, r’$\frac{3\pi}{2}$’, r’$2\pi$’])
ax.tick_params(axis=’x’, labelsize=10)
yticks

ax.set_xlabel(r’$x$ (rad)’, fontsize=12) ax.set_ylabel(r’$y$’, fontsize=12) ax.set_title(title[n], fontsize=12)
ylabel
title
xlabel
matplotlib

matplotlib
fig.set_size_inches(6,5)
fig.suptitle(rf’Trigonometric functions ($\Delta x =$ {step})’)
suptitle

import itertools
fig, axs = plt.subplots(2, 2)
matplotlib
for n, index in enumerate(itertools.product(range(2), repeat=2)): ax = axs[index]
ax.plot(x, func[n], ‘o-‘, linewidth=2, color=’#0000FF’)
ax.set_xlim([0, 2*np.pi])
ax.set_ylim([-1.2, 1.2])
ax.set_xticks([0, 0.5*np.pi, np.pi, 1.5*np.pi, 2*np.pi]) ax.set_xticklabels([r’$0$’, r’$\frac{\pi}{2}$’, r’$\pi$’,
r’$\frac{3\pi}{2}$’, r’$2\pi$’]) ax.set_yticks([-1, 0, 1]) ax.set_yticklabels([r’$-1$’, r’$0$’, r’$1$’]) ax.tick_params(axis=’x’, labelsize=10) ax.tick_params(axis=’y’, labelsize=10) ax.set_xlabel(r’$x$ (rad)’, fontsize=12) ax.set_ylabel(r’$y$’, fontsize=12) ax.set_title(title[n], fontsize=12) ax.grid(True)
fig.set_size_inches(6,5)
fig.suptitle(rf’Trigonometric functions ($\Delta x =$ {step})’) plt.tight_layout()
plt.show()
Type of plotting

2D plot Plot a function f(x) as a curve
E.g. f(x) = tanh(x)
x = np.arange(-10, 10+0.01, 0.01) _, axs = plt.subplots() axs.plot(x, np.tanh(x), linewidth=4, color=’#A81551′) plt.tight_layout()
plt.show()
Plot a polar equation r(θ)
E.g. r(θ) = sin(6θ) + 2
theta = np.arange(0, 2*np.pi+0.01, 0.01) _, axs = plt.subplots(subplot_kw ={‘projection’: ‘polar’}) axs.plot(theta, np.sin(6*theta) + 2, linewidth=4, color=’#A81551′) plt.tight_layout()
plt.show()
matplotlib

Plot data from a file as a curve. E.g. “HSI.csv”
fig, axs = plt.subplots()
axs.plot(stocks[“index”], stocks[“open”], linewidth=1, color=’#FFC000′, label=’open’) axs.plot(stocks[“index”], stocks[“high”], linewidth=1, color=’#A81551′, label=’high’) axs.plot(stocks[“index”], stocks[“low”], linewidth=1, color=’#0000FF’, label=’low’) axs.plot(stocks[“index”], stocks[“close”], linewidth=1, color=’#00B050′, label=’close’)
axs.legend()
plt.tight_layout()
plt.show()
matplotlib
2D plot

Plot histogram
E.g. Normal distribution of random numbers
matplotlib
mean = 1; stddev = 1; num_points = 10000
y = np.random.normal(mean, stddev, num_points)
fig, axs = plt.subplots()
axs.hist(y, bins=np.arange(-6, 8+0.2, 0.2), density=True, color=’#A81551′)
# bins are [-6,-5.8), [-5.8,-5.6), …, [7.6,7.8), [7.8,8] # density=True: normailized
plt.tight_layout()
plt.show()
2D plot

E.g. H(x,y) (x0.5)2 y2 1  (x0.5)2 y2 1
step = 0.01
x = np.arange(-1, 1+step, step)
y = np.arange(-1, 1+step, step)
X, Y = np.meshgrid(x, y)
H = (np.sqrt((X+0.5)**2+Y**2) – 1)**2\ + (np.sqrt((X-0.5)**2+Y**2) – 1)**2 fig, axs = plt.subplots() axs.contourf(X, Y, H,
levels=500, cmap=’rainbow’) axs.contour(X, Y, H, levels=np.arange(0,2,0.05)) plt.tight_layout()
plt.show()
matplotlib
2D plot
Plot contour  2  2

Import the library
from mpl_toolkits.mplot3d import Axes3D
E.g.
x(z2 1)sin(t) y(z 1)cos(t)
We can pick t = z

z5 z4 z3 z2 z1 z0
matplotlib
3D plot
A curve in 3D
A curve in 3D can be represented by a system of parametric equations
 x  x ( t )  y  y ( t ) z  z(t)
z z9
2z
zt 2π 
7 z6
x  x(z) y  y(z)
z  z 
z8

x(z2 1)sin(t) 2
y(z 1)cos(t) zt 2π

num_points = 100
t = np.linspace(-4 * np.pi, 4 * np.pi, num_points)
z = np.linspace(-2, 2, num_points)
r = z**2 + 1
x = r * np.sin(t)
y = r * np.cos(t)
fig, axs = plt.subplots(subplot_kw={‘projection’:’3d’}) axs.plot(x, y, z, ‘o-‘, linewidth=4, color=’#A81551’) plt.tight_layout()
plt.show()
matplotlib
3D plot
A curve in 3D

E.g.
y  sin(u)sin(v)
 x  c o s ( u ) s i n ( v ) z  cos(v)
We can pick u = x and v = y  x  x
y
 y  y
z  z(x,y)
x
matplotlib
3D plot
A surface in 3D
A surface in 3D can be represented by a system of parametric equations.
Different from the case of 3D curves, there are two parameters, u and v x  x(u,v)
 y  y ( u , v )
z  z(u,v) z

E.g.
H(x,y) (x0.5)2 y2 12  (x0.5)2 y2 12 
step = 0.002
x = np.arange(-1, 1+step, step)
y = np.arange(-1, 1+step, step)
X, Y = np.meshgrid(x, y)
H = (np.sqrt((X+0.5)**2+Y**2) – 1)**2\ + (np.sqrt((X-0.5)**2+Y**2) – 1)**2 fig, axs = plt.subplots(subplot_kw={‘projection’: ‘3d’})
axs.plot_surface(X, Y, H,
shade=True, linewidth=0) plt.tight_layout()
plt.show()
matplotlib
3D plot
A surface in 3D

matplotlib
Ellipse
Basic quadric surfaces Ellipsoid
 x 2  y 2  z 2
   1 a b c
Note the relationship: sin2 θ  cos2 θ  1
 x  a c o s ( u ) s i n ( v ) y  bsin(u)sin(v) z  ccos(v)

matplotlib
Hyperbola
Basic quadric surfaces Hyperboloid of one sheet
 x 2  y 2  z 2
   1 a b c
Note the relationships:
sin2 θ  cos2 θ  1
cosh2 θ  sinh2 θ  1 hyperbolic functions
 x  a c o s ( u ) c o s h ( v ) y  bsin(u)cosh(v) z  csinh(v)

matplotlib
Hyperbola
Basic quadric surfaces Hyperboloid of two sheets
 x 2  y 2  z 2
   1
a b c
Note the relationships: sin2 θ  cos2 θ  1
cosh2 θ  sinh2 θ  1
 x  a c o s ( u ) s i n h ( v ) y  bsin(u)sinh(v) z  ccosh(v)
two branches

Basic quadric surfaces Elliptic paraboloid
 x 2  y 2 z  
Ellipse x‐y
a b y
Parabola y‐z
Parabola z z‐x
z=1
x
matplotlib
z x=0yy=0x

Basic quadric surfaces Hyperbolic paraboloid
 x 2  y 2 z  
Hyperbola x‐y
y
Parabola y‐z
z
Parabola z z‐x
a b
z=2
x
z = ‒2
x=0
yy=0x
matplotlib

# convert an RGB image to a grayscale image
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
matplotlib
Open, read and plot an image
def rgb2gray(rgb_img):
r, g, b = rgb_img[:, :, 0], rgb_img[:, :, 1], rgb_img[:, :, 2] return 0.2990 * r + 0.5870 * g + 0.1140 * b
img = mpimg.imread(“./monet.jpg”)
print(type(img)) # numpy array print(img.shape) # 3D array: pixels, pixels, 3 gray_img = rgb2gray(img)
plt.imshow(gray_img, cmap=plt.get_cmap(‘gray’)) plt.tick_params(axis=’both’, which=’both’, bottom=False, top=False, right=False, left=False, labelbottom=False, labelleft=False)
plt.axis(‘off’)
plt.savefig(“monet_gray.png”, dpi=1200, bbox_inches=’tight’, pad_inches=0.0)
plt.show()

matplotlib
Open, read and plot an image
monet.jpg (RGB) monet_gray.png (grayscale)
After we transform the image to a numpy array, we can process the image by maths, e.g. Fourier transform, convolution, photoshop, machine learning, etc.

https://www.sympy.org/en/index.html
Installation
> pip install sympy
Import
# at the head of a code
from sympy import *

Evaluate the indefinite integral eax2 dx with a > 0
Answer:
Taylor expansion of sin(ax) about x = 0
Answer:
Solve the differential equation
(2n1)!
d f(x) 0, f(x 0)a, f(x L)b
Answer:
dx2
f (x)  (b  a)x / L  a
sin(ax)  n0
SymPy
Computation of mathematical objects symbolically (analytical solution) We usually see this kind of problems in math exams
E.g.
Solve the equation 4a2x2 + 4ax + 1 = 0, x is unknown and a is parameter Answer: if a = 0, no solution; else, x = ‒1/2a
π erf axconstant 4a
n 2n1 (1) (ax)
2

Wolfram Mathematica Powerful but not free
https://www.wolfram.com/mathematica
WolframAlpha
Easy access and flexible syntax
Free use limited by computation time
https://www.wolframalpha.com/
SymPy (Python library) Powerful and free
brief introduction in this lecture
More details:
https://docs.sympy.org/latest/index.html
SymPy
Computation of mathematical objects symbolically (analytical solution) We can let computers to solve the exam problems for us

pprint((1/2)**2 + 8/10)
# “pprint” is printing in pretty form pprint(Rational(1,2)**2 + Rational(8,10)) pprint(sqrt(1 + Rational(3,2)**2)) pprint(sqrt(1.1 + Rational(3,2)**2))
M1 = Matrix([[1, 2, 3],
[3, 2, 1]])
M2 = Matrix([0, 1, 1])
# by default, a 1D list is considered as a column vector M3 = Matrix([[Rational(1,3), 0, 1],
[2, -sin(4), 3],
[4, 3, 2]]) pprint(M3.col(0)); pprint(M3.row(1))
# access columns/rows
# multiplied by a scalar
# matrix product
# transpose
# determinant
# inverse
pprint(Rational(1,2) * M1)
pprint(M1 * M2)
pprint(M1.T)
pprint(M3.det())
pprint(M3**-1)
SymPy
Try simple algebra of numbers
Try simple algebra of matrices

Define constants
coeff1 = 2
coeff2 = 4.1
Define symbols
x, y, z = symbols(‘x y z’, real=True)
a, b, c = symbols(‘a b c’, real=True, positive=True)
Define expressions
expr = coeff1 * x**2 + coeff2 * sqrt(abs(y)) + sin(z) # note: “sqrt” is “sympy.sqrt”, “sin” is “sympy.sin” print(type(expr)) # type is “sympy.core.add.Add” pprint(expr)
SymPy
Definition of symbols and expressions
Substitution and evaluation
pprint(expr.subs(x,2)) # substitute x by 2 pprint(expr.subs([(x,2), (y,3), (z,4)]))# multiple substitution
# evaluate an expression
pprint(expr.subs([(x,2), (y,3), (z,4)]).evalf()) pprint(expr.evalf(subs={x:2, y:3, z:4}))

func_name = lambdify([x, y, z], expr, “numpy”)
”’
def func_name([x, y, z]):
SymPy
Convert SymPy expressions to Python functions
coeff1 = 2
coeff2 = 4.1
x, y, z = symbols(‘x y z’, real=True)
expr = coeff1 * x**2 + coeff2 * sqrt(abs(y)) + sin(z)
return expr
the input [x, y, z] and the output expr are all NumPy arrays ”’
x_arr = np.arange(-1, 1, 0.01)
y_arr = np.arange(-1, 1, 0.01)
z_arr = np.arange(-1, 1, 0.01) print(func_name(x_arr, y_arr, z_arr))

pprint(diff(c + a*exp(b*x*y*z), x, y, y, z, z, z, z))
 ax2
e dx,
   
Integrals
  x 2  y 2   x
SymPy
Derivatives
d327   aebx , c  aebxyz
dx3 xy2z4 pprint(diff(a*exp(b*x**2), x, x, x))
exp  dxdy, xdx ab

pprint(integrate(exp(-a*x**2), x)) # indefinite integral pprint(integrate(exp(-(x/a)**2 – (y/b)**2), (x, -oo, oo), (y, -oo, oo))) # definite integral
# what if the integral cannot be computed?
pprint(integrate(x**x, x))

x0 pprint(limit(expr_sinc, x, 0))
x
expr_sinc = sin(x)/x
pprint(expr_sinc.subs(x, 0)) # nan: not a number
SymPy
Limits
lim sin(x)
Series expansion
f (x)  sin(ax), f (x,y)  x4  y3 ex2 y2 12
Expand f1(x) about x = 0 to the 10th order and f2(x, y) about x = 0, y = 0 to the 10th order
pprint((sin(a*x)).series(x, 0, 10))
expr_expxy = (x**4 + y**3) * exp(x**2 + y**2) pprint(expand(expr_expxy.series(x, 0, 10).removeO().series(y, 0, 10).removeO()).coeff(x, 8).coeff(y, 9)) # expand the expression and pick the coefficient of the x8y9 term

Summation of geometric series
n
zk  km
zm zn1 1z
k, m, n = symbols(‘k m n’, integer=True, nonnegative=True) z = symbols(‘z’, real=True)
expr = z**k
pprint(simplify(summation(expr, (k, m, n))))
SymPy
Summation of sequence
(4k)! 1 1z

k0 24k 2(2k)!(2k1)!
zk 
z
, |z|1
k = symbols(‘k’, integer=True, nonnegative=True)
z = symbols(‘z’, real=True)
expr = factorial(4*k) * z**k / (2**(4*k) * sqrt(2) * factorial(2*k) * factorial(2*k + 1)) pprint(simplify(summation(expr, (k, 0, oo))))

x = symbols(‘x’, real=True) pprint(solveset(Eq(sin(x), 1), x))
SymPy
Solve equations sin(x)  1
Solve ordinary differential equations
2
d f(x)0, f(x0)a, f(xL)b
dx2
x, y = symbols(‘x y’, real=True)
a, b, L = symbols(‘a b L’, real=True)
f, g = symbols(‘f g’, cls=Function)
# f and g are undefined functions
pprint(f(x).diff(x, x)) # derivative without evaluation pprint(g(x, y).diff(x, x, y))
eq = Eq(f(x).diff(x, x), 0)
pprint(dsolve(eq, f(x), ics={f(0):a, f(L):b}))
# ics: take in initial/boundary conditions in form of dictionary

SymPy
Solve ordinary differential equations: An example
From Lecture 0: damped oscillation
r1 r2
0 l0 Ordinary differential equation
d2 (r  l )   k (r  l )  β d(r  l ) 20 20
dt2 m20mdt Spring force Damping
Initial conditions 2
r(0)r , dr v
2 2,0 dt 2,0
t0

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
# r2(t) t = symbols(‘t’, real=True, nonnegative=True) # time t
r = symbols(‘r’, cls=Function)
SymPy
Damped oscillation
b, m, k, l0, r0 = symbols(‘b m k l0 r0’, real=True, positive=True) v0 = symbols(‘v0’, real=True)
eq = Eq(r(t).diff(t,t) + (b/m)*r(t).diff(t) + (k/m)*(r(t)-l0), 0) sol = dsolve(eq, r(t), ics={r(0):r0, r(t).diff(t).subs(t,0):v0}, simplify=True).rhs # dsolve returns “equation”; its right-hand side “rhs” is expression
t_arr = np.arange(0, 80, 0.1)
sol = sol.subs([(m,1), (k,0.04), (l0,1), (r0,0.5), (v0,0)]) fig, axs = plt.subplots()
b_list = [0.08, 0.4, 1]
for b_val in b_list:
sol_b = sol.subs(b, b_val)
rt_func = lambdify(t, sol_b, “numpy”)
axs.plot(t_arr, rt_func(t_arr).real) plt.show()

Underdamped
Critically damped
SymPy
Damped oscillation
Overdamped



Final remarks
For most cases of practical interest in science and engineering, symbolic/analytical solutions do not exist
E.g. motion of large amounts of atoms

dt2 11 N dt
2
mdr F(r,,r)βdr
11
2
mdr F(r,,r)βdr dt2 21 N dt
22
 NN
2
mdr F(r,,r)βdr dt2 N1 N dt
SymPy doesn’t work. The functions {Fi} is too complicated. Even with simple {Fi}, the number of ODEs is too large
This problem can be solved numerically (molecular dynamics)