assignment_5
MSE3114 Computational Methods for Physicists and Materials Engineers
Assignment 5
Q1. (Vector norm) Vector norm is a measure of the “size” of a vector. We have learned three norms: $\ell_1$-norm, $\ell_2$-norm and $\ell_\infty$-norm. In this question we will understand these norms by some examples.
Figure (a) shows the google map of Manhattan of New York City. Consider how to measure the distance between two locations on this map. To simplify the problem, we only consider the simplified map shown in Fig. (b) and measure the distance between “Start” point and “End” point. Suppose that the address of “Start” is $i$-th St $j$-th Ave, the address of “End” is $p$-th St $q$-th Ave, and the size of block is $L_1$ by $L_2$. Then, the coordinate of “Start” is $\mathbf{x}_\text{start} = (iL_1, jL_2)$ and the coordinate of “End” is $\mathbf{x}_\text{end} = (pL_1, qL_2)$. The straight distance between “Start” and “End” is the $\ell_2$-norm: $\lVert \mathbf{x}_\text{end} – \mathbf{x}_\text{start}\rVert_2$; i.e., the blue line in Fig. (b). Write a python code to define a function straight_distance(i, j, p, q, L1, L2), which takes in the values of $i$, $j$, $p$, $q$, $L_1$ and $L_2$ and returns the straight distance.
However, a taxi cab cannot run along the straight line; it must run along the streets and avenues. For the taxi cab, the distance between “Start” and “End” is the $\ell_1$-norm: $\lVert \mathbf{x}_\text{end} – \mathbf{x}_\text{start}\rVert_1$; i.e., the red line in Fig. (a). Define a function taxi_cab_distance(i, j, p, q, L1, L2), which takes in the values of $i$, $j$, $p$, $q$, $L_1$ and $L_2$ and returns the distance for a taxi cab.
If the address of “Start” point is 10th St 4th Ave, the address of “End” point is 3th St 7th Ave, and the size of block is $L_1 = 300$ m and $L_2 = 100$ m. Apply the function straight_distance() defined in Q1.1 to calculate the straight distance. Apply the function taxi_cab_distance() defined in Q1.2 to calculate the distance for a taxi cab.
Overhead crane is used to move an object from one place to another; see Fig. (a). The object can be moved along $\mathbf{e}_1$- and $\mathbf{e}_2$-directions simultaneously with the same velocity $v$. Consider how much time an overhead crane takes to move an object for “Start” point to “End” point in Fig. (b). Suppose that the coordinate of “Start” is $\mathbf{x}_\text{start} = (x_1, x_2)$ and the coordinate of “End” is $\mathbf{x}_\text{end} = (y_1, y_2)$. The time for the overhead crane takes to move an object from $\mathbf{x}_\text{start}$ to $\mathbf{x}_\text{end}$ is $t = \lVert \mathbf{x}_\text{end} – \mathbf{x}_\text{start}\rVert_\infty / v$, i.e., the $\ell_\infty$-norm divided by the velocity $v$. Define a function overhead_crane_time(x1, x2, y1, y2, v), which takes in the coordinates $(x_1, x_2)$ and $(y_1, y_2)$ and the velocity $v$ and returns the time for an overhead crane to move an object from $(x_1, x_2)$ to $(y_1, y_2)$. Apply this function to calculate the time to move an object from $(2.1, 0.3)$ to $(11.3, 31.0)$ (unit: m) with the velocity $v = 0.1$ m/s.
(25 marks)
Q2. (Contraction) If an operator $\mathbf{f}$ is a contraction, there must exist $q \in [0, 1)$ such that
$$
\lVert \mathbf{f}(\mathbf{x}) – \mathbf{f}(\mathbf{y}) \rVert
\le q \lVert \mathbf{x} – \mathbf{y} \rVert
$$
for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}$.
You can choose any norm convenient for your calculation, e.g. $\ell_1$-norm, $\ell_2$-norm, $\ell_\infty$-norm, etc. Check if the following operators are contractions. (You don’t need python code for solving this question.)
Consider the following linear equations:
$$
\mathbf{f}(\mathbf{x})
\equiv \left(\begin{array}{c}
0.1 x_1 + 0.2 x_2 – 10 \\
0.3 x_1 – 0.6 x_2 + 1
\end{array}\right)
= \left(\begin{array}{c}
x_1 \\ x_2
\end{array}\right)
\equiv \mathbf{x}.
$$
Is $\mathbf{f}$ a contraction for $\mathbf{x} \in \mathbb{R}$? Prove your answer.
Hint: (i) you can use $\ell_1$-norm or $\ell_\infty$-norm. (ii) Recall the inequality: $|x + y| \le |x| + |y|$ and $|x – y| \le |x| + |-y| = |x| + |y|$.
Kepler’s equation is
$$
M + e \sin E = E,
$$
where $M$ is the mean anomaly, $E$ is the eccentric anomaly and $e$ is the eccentricity.
Consider the equation $f(E) = E$, where $f(E) \equiv M + e\sin E$.
Is $f$ a contraction for $E \in \mathbb{R}$? Prove your answer.
Hint: (i) Note the relationship:
$$
\sin x – \sin y = 2 \sin\left(\frac{x – y}{2}\right) \cos\left(\frac{x + y}{2}\right).
$$
(ii) Note the inequalities: $|\sin x| \le |x|$ and $|\cos x| \le 1$.
(iii) Finally, you may find that if $f$ is a contraction depends on the value of $e$.
(25 marks)
Q3. (Jacobi and Gauss-Seidel methods) Write a python code for solving a system of linear equations by Jacobi method and Gauss-Seidel method. Written in matrix form, a system of linear equations is expressed as $\mathbf{A}\mathbf{x} = \mathbf{b}$.
Define a function jacobi(A, b, x, eps, max_n), where A is the matrix $\mathbf{A}$, b is the vector $\mathbf{b}$, x is the initial guess of the solution, eps is the criterion for $\ell_2$-norm of residual (i.e., $\lVert \mathbf{x}_\text{current_step} – \mathbf{x}_\text{previous_step}\rVert_2$), and max_n is the maximum iteration steps allowed. This function solves the equations by Jacobi iteration. It returns the solution when $\lVert \mathbf{x}_\text{current_step} – \mathbf{x}_\text{previous_step}\rVert_2 \lt$ eps, or when the iteration step is larger than max_n.
Define a function gauss_seidel(A, b, x, eps, max_n), where A is the matrix $\mathbf{A}$, b is the vector $\mathbf{b}$, x is the initial guess of the solution, eps is the criterion for $\ell_2$-norm of residual (i.e., $\lVert \mathbf{x}_\text{current_step} – \mathbf{x}_\text{previous_step}\rVert_2$), and max_n is the maximum iteration steps allowed. This function solves the equations by Gauss-Seidel iteration. It returns the solution when $\lVert \mathbf{x}_\text{current_step} – \mathbf{x}_\text{previous_step}\rVert_2 \lt$ eps, or when the iteration step is larger than max_n.
Solve the following equations:
$$
\left(\begin{array}{cccc}
20 & 1.3 & -2 \\
1 & -13 & -7 \\
0.1 & 0.7 & 3
\end{array}\right)
\left(\begin{array}{c}
x_1 \\ x_2 \\ x_3
\end{array}\right)
=
\left(\begin{array}{c}
3.3 \\ 10 \\ 52
\end{array}\right)
$$
by three methods: (i) calling the function jacobi(A, b, x, eps, max_n) defined in Q3.1 with the initial guess x = $\mathbf{0}$, eps $=10^{-10}$ and max_n $=1000000$, (ii) calling the function gauss_seidel(A, b, x, eps, max_n) defined in Q3.2 with the initial guess x = $\mathbf{0}$, eps $=10^{-10}$ and max_n $=1000000$, and (iii) applying the direct method in SciPy scipy.linalg.solve().
(25 marks)
Q4. (GMRES and CG) Write a python code to do the following operations.
Define a function create_sparse_rand_mat(N, zero_frac) for creating (i) A: a sparse random symmetric matrix $\mathbf{A}$ of size $N\times N$ with the fraction of non-zero elements zero_frac stored in format of compressed sparse row matrix (CSR), (ii) A_arr: the same matrix $\mathbf{A}$ stored in format of numpy array, and (iii) b: a random matrix $\mathbf{b}$ of size $N\times 1$ by
A = scipy.sparse.random(N, N, zero_frac) # sparse random matrix in format of CSR
A = A.toarray() # convert CSR format to numpy array format
A_arr = A + np.transpose(A) + 4*np.eye(N) # addition of “4*np.eye(N)” guaranttes that A is well-conditioned
A = scipy.sparse.csr_matrix(A_arr) # convert numpy array format to CSR format
b = np.random.rand(N, 1)
Create the matrices $\mathbf{A}$ in both CSR format and array format and $\mathbf{b}$ by calling the function create_sparse_rand_mat(N, zero_frac) defined in Q4.1 with $N = 100$ and zero_frac $=0.005$. Plot the matrix $\mathbf{A}$ thus created by .imshow() in matplotlib.
Loop over the matrix size $N = 100, 150, 200, \cdots, 750, 800$. In each cycle with $N$, create random $\mathbf{A}$ and $\mathbf{b}$ by create_sparse_rand_mat(N, zero_frac) defined in Q4.1 with zero_frac $=0.005$. Count the times spent on runing scipy.sparse.linalg.gmres(A, b, x0, tol, maxiter), scipy.sparse.linalg.cg(A, b, x0, tol, maxiter), and scipy.linalg.solve(A, b). The parameters are x0 $=\mathbf{0}$, tol $=10^{-10}$ and maxiter $= 1000000$. You can learn how to count the time by time.time_ns() (in unit of nanosecond) from Page 13 of the note of lecture 2. For each cycle, append the values of $N$ and the three times mentioned above to four lists.
Plot the computation times spent on three methods vs. $N$. Show the legend to label the meaning of each curve.
(25 marks)