[SOLVED] CS python assignment_4

assignment_4

MSE3114 Computational Methods for Physicists and Materials Engineers

Assignment 4

Q1. (Gaussian elimination) Write a python code for solving a system of linear equations by Gaussian elimination. Written in matrix form, a system of linear equations is expressed as $mathbf{A}mathbf{x} = mathbf{b}$.

Define a function forward_elim(A, b) which takes in $mathbf{A}$ and $mathbf{b}$, does forward elimination, and returns the new $mathbf{A}$ and $mathbf{b}$ after foward elimination.
Define a function backward_subs(A, b) which takes in the $mathbf{A}$ and $mathbf{b}$ obtained by foward elimination, does backward substitution, and returns the solution $mathbf{x}$.
Define a function gaussian_elim(A, b) which takes in the initial $mathbf{A}$ and $mathbf{b}$ and returns the solution $mathbf{x}$. In this function, we need to call the two functions defined in Q1.1 and Q1.2.
Apply the function gaussian_elim(A, b) defined in Q1.3 to solve the following equations:
$$
left(begin{array}{cccc}
2 & 1 & -2 \
1 & -1 & -1 \
1 & 1 & 3
end{array}right)
left(begin{array}{c}
x_1 \ x_2 \ x_3
end{array}right)
=
left(begin{array}{c}
3 \ 0 \ 12
end{array}right).
$$

(25 marks)

Q2. (LU decomposition) Write a python code for solving a system of linear equations by LU decomposition. Written in matrix form, a system of linear equations is expressed as $mathbf{A}mathbf{x} = mathbf{b}$. The pivoted LU decomposition on $mathbf{A}$ gives $mathbf{A} = mathbf{P}mathbf{L}mathbf{U}$. Then, the equations become $mathbf{P}mathbf{L}mathbf{U}mathbf{x} = mathbf{b}$. We can firstly solve $mathbf{L}mathbf{z} = mathbf{P}^Tmathbf{b}$ for $mathbf{z}$ by the forward substitution, and finally solve $mathbf{U}mathbf{x} = mathbf{z}$ for $mathbf{x}$ by the backward substitution.

Define a function plu_decomposition(A) which takes in $mathbf{A}$, does pivoted LU decomposition by scipy.linalg.lu(), and returns the permutation matrix $mathbf{P}$, the lower triangular matrix $mathbf{L}$ and the upper triangular matrix $mathbf{U}$.
Define a function forward_subs(L, Pb) which takes in $mathbf{L}$ and $mathbf{P}^Tmathbf{b}$, does forward substitution, and returns the result after forward substitution $mathbf{z}$.
Define a function solve_by_lu_decomp(A, b) which takes in $mathbf{A}$ and $mathbf{b}$, does LU decomposition by calling plu_decomposition(A) defined in Q2.1, print out the result of LU decomposition (i.e., $mathbf{P}$, $mathbf{L}$ and $mathbf{U}$), does forward substitution by calling forward_subs() defined in Q2.2 on $mathbf{L}$ and $mathbf{P}^Tmathbf{b}$ and returns $mathbf{z}$, does backward_substitution by calling backward_subs() defined in Q1.2 on $mathbf{U}$ and $mathbf{z}$ and returns the solution $mathbf{x}$.
Apply the function solve_by_lu_decomp(A, b) defined in Q2.3 to solve the following equations:
$$
left(begin{array}{cccc}
3 & 2 & -3 \
1 & -4 & -1 \
7 & 1 & 31
end{array}right)
left(begin{array}{c}
x_1 \ x_2 \ x_3
end{array}right)
=
left(begin{array}{c}
1 \ 9 \ -12
end{array}right).
$$

Solve the same equations in Q2.4 by scipy.linalg.solve() directly.

(25 marks)

Q3. (QR decomposition) Write a python code for solving a system of linear equations by QR decomposition. Written in matrix form, a system of linear equations is expressed as $mathbf{A}mathbf{x} = mathbf{b}$. The QR decomposition on $mathbf{A}$ gives $mathbf{A} = mathbf{Q}mathbf{R}$. Then, the equations become $mathbf{Q}mathbf{R}mathbf{x} = mathbf{b}$. We can solve $mathbf{R}mathbf{x} = mathbf{Q}^Tmathbf{b}$ for $mathbf{x}$ by the backward substitution.

Define a function qr_decomposition(A) which takes in $mathbf{A}$, does QR decomposition by scipy.linalg.qr(), and returns the orthogonal matrix $mathbf{Q}$ and the upper triangular matrix $mathbf{R}$.
Define a function solve_by_qr_decomp(A, b) which takes in $mathbf{A}$ and $mathbf{b}$, does QR decomposition by calling qr_decomposition(A) defined in Q3.1, print out the result of QR decomposition (i.e., $mathbf{Q}$ and $mathbf{R}$), does backward_substitution by calling backward_subs() defined in Q1.2 on $mathbf{R}$ and $mathbf{Q}^Tmathbf{b}$ and returns the solution $mathbf{x}$.
Apply the function solve_by_qr_decomp(A, b) defined in Q3.2 to solve the following equations:
$$
left(begin{array}{cccc}
13 & 2 & -3 \
1 & -44 & -1 \
7 & 1 & 1
end{array}right)
left(begin{array}{c}
x_1 \ x_2 \ x_3
end{array}right)
=
left(begin{array}{c}
1 \ 0.9 \ -1
end{array}right).
$$

Solve the same equations in Q3.3 by scipy.linalg.solve() directly.

(25 marks)

Q4. From the previous questions, we have four methods to solve the system of linear equations: gaussian_elim() defined in Q1, solve_by_lu_decomp() defined in Q2, solve_by_qr_decomp() defined in Q3 and scipy.linalg.solve() in SciPy. Now, we will compare the efficiencies of the four methods.

Copy the functions gaussian_elim(), solve_by_lu_decomp(), solve_by_qr_decomp() defined in the previous questions and all the other functions on which they depend.

Define a function create_rand_mat(N) for creating a random invertible matrix $mathbf{A}$ of size $Ntimes N$ and a random matrix $mathbf{b}$ of size $Ntimes 1$. The method is as follows. First, create a random lower triangular matrix:

L_rand = np.tril(np.random.rand(N, N))

Second, create a random upper triangular matrix:

U_rand = np.triu(np.random.rand(N, N))

Third, create the random matrix $mathbf{A}$:

A = np.matmul(L_rand, U_rand) + 2*np.eye(N)

The $mathbf{A}$ created in this way must be invertible. np.eye(N) creates a $Ntimes N$ identity matrix. Addition of 2*np.eye(N) guarantees that $mathbf{A}$ is well-condictioned or the inverse matrix $mathbf{A}^{-1}$ can be solved with high accuracy. Finally, create the random matrix $mathbf{b}$:

b = np.random.rand(N, 1)

Loop over the matrix size $N = 100, 110, 120, cdots, 190, 200$. In each cycle with $N$, create random $mathbf{A}$ and $mathbf{b}$ by create_rand_mat(N) defined in Q4.2. Count the times spent on runing gaussian_elim(A, b), solve_by_lu_decomp(A, b), solve_by_qr_decomp(A, b) and scipy.linalg.solve(A, b). 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^3$ and the four times mentioned above to five lists.

Plot two figures. In the first figure, plot the times spent on four methods vs. $N^3$; so, you will show four curves. In the second figure, ignore the data of the time spent on gaussian_elim() and plot the times spent on the other three methods vs. $N^3$; so, you will show three curves. Show the legend on each figure to label the meaning of each curve. According to the figures, rank the efficiencies of the four methods from high to low.

(25 marks)