FIZ353 - Numerical Analysis | 30/10/2020
Emre S. Tasci emre.tasci@hacettepe.edu.tr
import numpy as np
np.random.seed(353) # we are fixing the random seed,
# so you'll get the same "random" numbers
# as of this lecture note's
M = np.array([[1,2,3],[4,5,6],[7,8,9]])
M
Even though they are much easier to define (and afterwards converted to np.arrays), the usage of np.matrix class are disencouraged because:
It is no longer recommended to use this class, even for linear algebra. Instead use regular arrays. The class may be removed in the future. Source
MM = np.matrix("1 2 3; 4 5 6; 7,8,9")
MM
# either like this:
M = np.array(MM)
# or like this:
M = MM.A
M
O = np.zeros((3,4))
O
L = np.ones((3,4)) * 3.2
L
I = np.eye(3,4)
I
D = np.diag((5,3,7,2))
D
# Upper
DU = np.diag((5,3,7,2),1)
DU
# Lower
DL = np.diag((5,3,7,2),-1)
DL
# Uniform distribution [0,1)
R1 = np.random.random((4,5))
R1
# Uniform distribution, integers [a,b)
R2 = np.random.randint(6,10,(4,5))
R2
M = np.array([[1,2,3],[4,5,6],[7,8,9]])
M
M[2,2]
M[2,:]
M[:,[0,2]]
M[:,1] = -1
M
When copying an array object, be careful that they are not copied over like scalars:
M = np.array([[1,2,3],[4,5,6],[7,8,9]])
M
N = M
N
M[2,2] = -9
M
N
Instead, use the copy() method:
M = np.array([[1,2,3],[4,5,6],[7,8,9]])
M
N = M.copy()
N
M[2,2] = -9
M
N
M
M>5
# The elements that are greater than 5:
M[M > 5]
# The elements that are even:
M[M % 2 == 0]
M>=3
M<=8
np.logical_and(M>=3, M<=8)
# The elements that are between 3 and 8 (inclusive):
M[np.logical_and(M>=3, M<=8)]
# The elements that are smaller than 3 or greater than 7 (inclusive)
M[np.logical_or(M<=3, M>=7)]
# Intersection of two arrays:
N = np.random.randint(-5,10,(3,3))
N
np.intersect1d(M,N)
M
5.21 * M
v = np.array([[9],[8],[7]])
v
np.dot(M,v)
v.T # 'T' for "transpose"
np.dot(v.T,M)
Be very careful!
np.dot(M,v) is not equal to M*v and,
np.dot(v.T,M) is not equal to v.T*M
Multiplying
Compare the following results with the above ones!
M
v
M*v
v.T*M
(These just do element-wise multiplications!)
M = np.array([[2,1,3],[4,6,3],[7,8,9]])
M
M_det = np.linalg.det(M)
print("Determinant of M is: {:.3f}".format(M_det))
M_inv = np.linalg.inv(M)
M_inv
np.dot(M,M_inv)
np.dot(M_inv,M)
Consider a (3x1) vector (let's think it as a position vector. It can be multiplied from left by a (3x3) matrix, and the result will be a (3x1) vector. We can interpret this operation as a transformation operation on the vector. As, by definition, a vector has a magnitude and a direction, a transformation can change these two properties. First, let's investigate these effects seperately by operating on a 2-dimensional vector, and then we'll combine them.
This is very straightforward as we can change a vector's magnitude by simply multiplying it with a scalar.
$\vec{v} = 3\hat{\imath}+4\hat{\jmath}$
v = np.array([3,4]) # a 2 dimensional vector
v
# It's magnitude:
v_mag = np.linalg.norm(v)
print("The magnitude of vector v: {:}".format(v_mag))
# We could as well have calculated it via Pythagorean Theorem:
v_mag = (v[0]**2 + v[1]**2)**0.5
print("The magnitude of vector v: {:}".format(v_mag))
# or, more practically then the above
# (for which, we have to enter each dimension exclusively)
# we could (and should) go for direct evaluation:
v_mag = (np.sum(v*v))**0.5
print("The magnitude of vector v: {:}".format(v_mag))
(examine and understand how we have done - this is better than the second solution but using predefined functions such as norm is almost always more optimized)
Now that we know our vector's initial magnitude, let's multiply it by a scalar, s = 2 and check the resultant vector $\vec{r}$'s magnitude:
s = 2
r = 2*v
print(r)
print("The magnitude of vector r: {:}".format(np.linalg.norm(r)))
How can we write this size change operation as a multiplication of a matrix by the vector? It will sound dumb, but we can insert the identity matrix in the middle, as it doesn't change a vector it is being multiplied to, so nothing should change:
$$ s\cdot\vec{v} = \underbrace{s\cdot(\mathbb{1})}_{S}\cdot\vec{v}$$From here, we can apply the scalar to the identity vector, which will be equal to a diagonal vector with the diagonal values being equal to $s$:
I = np.eye(2,2) # as our vector is 2 dimensional,
# corresponding identity vector is (2x2)
print(I)
S = s*I
S
The multiplication has now taken the form of a matrix multiplied by a vector:
$$\underbrace{\begin{bmatrix}2 & 0 \\ 0 & 2\end{bmatrix}}_{S}\cdot\underbrace{\begin{bmatrix}3 \\ 4\end{bmatrix}}_{\vec{v}}$$r = np.dot(S,v)
print(r)
print("The magnitude of vector r: {:}".format(np.linalg.norm(r)))
So, in summary, if we change the size of a vector, we just multiply it with a diagonal matrix, whose diagonal elements' values are equal to the ratio of the change of size.
The rotation operator R for a 2D vector is defined as:
$$ R_{\theta} = \begin{pmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{pmatrix}$$where $\theta$ is the rotation angle, counter-clockwise.
So, if we would rotate a $\vec{v}$ vector of, let's say, 5 units long along the positive x-direction by an angle of $90^o$, we would get a vector with the same magnitude, pointing alont the positive y-direction:
$$R_{90^o}.\vec{v} = R_{90^o}.(5\hat{\imath}) = 5\hat{\jmath}$$Let's put this in action:
theta = np.deg2rad(90)
R_theta = np.array([[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
R_theta
v = np.array([[5],[0]])
v
v_rotated_90 = np.dot(R_theta,v)
v_rotated_90
# Let's rotate this resultant vector once again by 90 degrees:
v_rotated_90_90 = np.dot(R_theta,v_rotated_90)
v_rotated_90_90
We can write the last operation explicitly as:
$$R_{90^o}.\vec{v}' =R_{90^o}.\left(R_{90^o}\vec{v}\right)$$changing the order of multiplication:
$$R_{90^o}.\vec{v}' =\left(R_{90^o}.R_{90^o}\right)\vec{v} = R_{90^o}^2.\vec{v}$$R_90_squared = np.dot(R_theta,R_theta)
R_90_squared
v_rotated_90_90 = np.dot(R_90_squared,v)
v_rotated_90_90
We can rotate a vector, and then change its magnitude, or, we can change its magnitude and then rotate it -- from physical point of view, the order shouldn't matter. And as it turns out perfectly, even though $A\cdot B \ne B\cdot A$ almost always, for these two rotation & scaling operation matrices, both ordering yields the same result (you can easily prove this for two general (2x2) matrices by hand). So, we can write a general transformation matrix $T$ as a multiplication (combination) of the two:
theta = np.deg2rad(30)
R_30 = np.array([[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
R_30
S_3 = 3*np.eye(2,2)
S_3
T = np.dot(R_30,S_3)
T
T = np.dot(S_3,R_30)
T
What we have ourselves here is a transformation matrix that rotates a vector by $30^o$ and triples its size. Let's roll! 8)
# For our first case, let's transform the unit vector
# along the positive x-direction:
v = np.array([[1],[0]])
v
vp = np.dot(T,v)
vp
# Its size:
np.linalg.norm(vp)
# We can calculate its angle wrt x-axis by calculating
# the arctan (vp_y/vp_x):
angle_rad = np.arctan(vp[1]/vp[0])
print("The angle in radians: {:}".format(angle_rad))
print("The angle in degrees: {:}".format(np.rad2deg(angle_rad)))
(For a much more "intelligent" arctan process, I strongly recommend the arctan2 function, which saves you from further consideration of the quadrants (i.e., Q1 vs. Q3 and Q2 vs. Q4)
angle_rad = np.arctan2(vp[1],vp[0])
print("The angle in radians: {:}".format(angle_rad))
print("The angle in degrees: {:}".format(np.rad2deg(angle_rad)))
So, the size is indeed tripled and our vector has been rotated by 30 degrees! 8)
## Now let's roll with a random vector:
v = (np.random.rand(1,2)*10).T
v
# initial size:
size_before = np.linalg.norm(v)
print("The vector's size is: {:}".format(size_before))
angle_deg_before = np.rad2deg(np.arctan2(v[1],v[0]))
print("The vector's angle (wrt x-axis) is: {:.3f} degress"
.format(angle_deg_before[0]))
# Transform it!
vp = np.dot(T,v)
vp
# Its size after the transformation:
size_after = np.linalg.norm(vp)
print("The vector's size after transformation is: {:}".format(size_after))
# The ratios of after/before:
print("The size has changed by a factor of: {:}".format(size_after/size_before))
# Its angle after the transformation:
angle_deg_after = np.rad2deg(np.arctan2(vp[1],vp[0]))
print("The vector's angle (wrt x-axis) is: {:.3f} degress"
.format(angle_deg_after[0]))
# The difference between the final and initial angles:
print("The vector has rotated by {:.3f} degrees"
.format(angle_deg_after[0] - angle_deg_before[0]))
Suppose that we are given a vector $\vec{v}'$ and we are told that this vector has been rotated by an angle of $\theta$ and its magnitude is scaled by a factor of $n$, and we are asked to calculate the initial vector $\vec{v}$. So here is what we have:
$$\vec{v}' = \begin{bmatrix}n & 0 \\ 0 & n\end{bmatrix}\begin{bmatrix}\cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{bmatrix}\vec{v}=\underbrace{\begin{bmatrix}n\cos\theta & -n\sin\theta \\ n\sin\theta & n\cos\theta\end{bmatrix}}_{T}\vec{v}$$To find $\vec{v}$, we would first rotate it (back) clockwise by 30 degrees (i.e., $-30^0$) and then would scale it $\vec{v}'$ by a factor of $\tfrac{1}{3}$. In mathematical representation this would be:
$$\vec{v} = \begin{bmatrix}\tfrac{1}{n} & 0 \\ 0 & \tfrac{1}{n}\end{bmatrix}\begin{bmatrix}\cos{(-\theta)} & -\sin{(-\theta)} \\ \sin{(-\theta)} & \cos{(-\theta)}\end{bmatrix}\vec{v}'=\underbrace{\begin{bmatrix}\tfrac{1}{n}\cos\theta & \tfrac{1}{n}\sin\theta \\ -\tfrac{1}{n}\sin\theta & \tfrac{1}{n}\cos\theta\end{bmatrix}}_{T'}\vec{v}'$$Let's put some numbers in it, for example $n=3$ & $\theta=30^o$, $\vec{v}' = 1.7942\hat{\imath} + 14.8923\hat{\jmath} $:
n = 3
angle = np.deg2rad(30)
T = np.array([[n*np.cos(angle),-n*np.sin(angle)],[n*np.sin(angle),n*np.cos(angle)]])
T
Now let's calculate the reversed (inverted?) version from the formula:
n = 1/3
angle = np.deg2rad(-30)
Tp = np.array([[n*np.cos(angle),-n*np.sin(angle)],[n*np.sin(angle),n*np.cos(angle)]])
Tp
Thus, by multiplying $T'$ with $\vec{v}'$, we get:
vp = np.array([[1.7942],[14.8923]])
vp
v = np.dot(Tp,vp)
v
which means that, our original vector was: $\vec{v} = 3\hat{\imath}+4{\jmath}$
Ready for a surprise?
Here is the inverse of the transformation matrix, directly calculated by the linalg.inv() method:
np.linalg.inv(T)
We saw it somewhere before, right? ;)
Thus:
$$\begin{align*}T\cdot \vec{v} &= \vec{v}'\\ \underbrace{T^{-1}\cdot T}_{\mathbb{1}}\cdot \vec{v} &=T^{-1}\cdot \vec{v}'\\ \vec{v} &=T^{-1}\cdot \vec{v}'\end{align*}$$np.dot(np.linalg.inv(T),vp)
(*As long as we're keeping the 'new coordinates' perpendicular)
We could have been directly given the transformation matrix T:
T = np.array([[ 2.59807621, -1.5],[ 1.5 , 2.59807621]])
T
and we could have been asked to extract $n$ & $\theta$. What we could have done would be to calculate the determinant of the transformation matrix, which would be:
$$n^2.\cos^2{\theta} + n^2.\sin^2{\theta} \equiv n^2$$Then, we would divide the matrix by $n$ to get the rotation matrix and thus, the angle:
n = np.sqrt(np.linalg.det(T))
print("The scaling factor is: {:.3f}".format(n))
R = T/n
R
angle = np.rad2deg(np.arctan2(R[1,0],R[0,0]))
print("The rotation angle is: {:.3f} degrees".format(angle))
Example $\vec{v}=\begin{pmatrix}v_x\\v_y\end{pmatrix}$ vector is rotated and scaled by a transformation matrix defined as:
$$T=\begin{bmatrix}2.59807621&-1.5\\1.5 &2.59807621\end{bmatrix}$$and if the resulting vector $\vec{v}'$ is given as $\vec{v}' = 1.7942\hat{\imath} + 14.8923\hat{\jmath} $, calculate the initial vector $\vec{v}$.
So we can write it as: $$T\cdot \vec{v}=\vec{v}'$$ $$\begin{bmatrix}2.59807621&-1.5\\1.5 &2.59807621\end{bmatrix}\begin{pmatrix}v_x\\v_y\end{pmatrix}=\begin{pmatrix}1.7942\\14.8923\end{pmatrix}$$
But we have already solved this problem and found the answer to be:
$$\vec{v}=\begin{pmatrix}v_x\\v_y\end{pmatrix}=\begin{pmatrix}3\\4\end{pmatrix} = 3\hat{\imath} + 4\hat{\jmath}$$How very nice! 8)
Now consider this set of two equations with two unknowns: $$\begin{gather*}4x - y = 10\\3x + 2y=13\end{gather*}$$
Can we write it as a multiplication of a matrix by a vector like the following?
$$\begin{bmatrix}4&-1\\3&2\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix} = \begin{bmatrix}10\\13\end{bmatrix}$$(please verify that, we indeed can...)
So, even though it is a set of 2 equations with two unknowns, we can interpret it as if it is a vector transformation! (albeit, a transformation, that doesn't keep the angle between the new axes perpendicular)
Sooooo:
vp = np.array([[10],[13]])
T = np.array([[4,-1],[3,2]])
np.dot(np.linalg.inv(T),vp)
We have just solved our set of equations!
(As this solving of linear operations are very common in scientific calculations, there is a direct function called linalg.solve() defined for it, i.e.,
np.linalg.solve(T,vp)
A = np.array([[2,3],[4,5]])
A
[eigvalues, eigvectors] = np.linalg.eig(A)
print("Eigenvalues: ",eigvalues,"\n-----\nEigenvectors:\n",eigvectors)
E1 = eigvectors[:,0]
E2 = eigvectors[:,1]
l1 = eigvalues[0]
l2 = eigvalues[1]
vp1 = np.dot(A,E1)
vp1
np.cross(vp1,E1)
np.dot(vp1,E1)
-0.79681209* 0.21905736 + 0.60422718 * -0.16611246
ctheta = np.dot(vp1,E1)/(np.linalg.norm(vp1)*np.linalg.norm(E1) )
np.rad2deg(np.arccos(ctheta))
vp1 / l1
vp1 / E1
np.linalg.norm(E1)
np.diag(eigvalues)
np.linalg.multi_dot([eigvectors,np.diag(eigvalues),np.linalg.inv(eigvectors)])