<jupyter_start><jupyter_text>Lecture 3 Part 1 - Linear Equations - In-Class Exercises Exercise 1 - Gaussian elimination and back substitutionBelow you find part of an example program to solve the following set of equations: $$ \left(\begin{array}{cccc} 2 & 1 & 4 & 1 \\ 3 & 4 & -1 & -1 \\ 1 & -4 & 1 & 5 \\ 2 & -2 & 1 & 3 \end{array} \right) \left(\begin{array}{c} w \\ x \\ y \\ z\end{array} \right) = \left(\begin{array}{c} -4 \\ 3 \\ 9 \\ 7\end{array} \right)$$1. Complete the Gaussian elimination part of the program.2. Add a print statement that prints the matrix at every step to check that the program is eliminating the lower triangle correctly. You can add a read statement that waits for you to press enter to proceed.3. Check your final solution for the vector $\mathbf{x}=(w,x,y,z)$ after backsubstitution. Talking points1. What do you observe?2. Is your solution correct?Example program to complete:<jupyter_code>from numpy import array,empty

A = array([[ 2,  1,  4,  1 ],
           [ 3,  4, -1, -1 ],
           [ 1, -4,  1,  5 ],
           [ 2, -2,  1,  3 ]], float)
v = array([ -4, 3, 9, 7 ],float)
N = len(v)

def print_Av(A,v):
    for i in range(N):
        print(A[i], " ", v[i])
    print()

# Gaussian elimination
print_Av(A,v)

for m in range(N):
    # Divide by the diagonal element
    print("Dividing row", m, "by", A[m,m])
    v[m] /= A[m,m]
    A[m] /= A[m,m]
    print_Av(A,v)

    # Now subtract from the lower rows
    print("Subtracting row", m, "from the lower rows")
    for n in range(m+1, N):
        v[n] -= A[n,m] * v[m]
        A[n] -= A[n,m] * A[m]
    print_Av(A,v)

    
# Backsubstitution
x = empty(N,float)
for m in range(N-1,-1,-1):
    x[m] = v[m]
    for i in range(m+1,N):
        x[m] -= A[m,i]*x[i]<jupyter_output>[2. 1. 4. 1.]   -4.0
[ 3.  4. -1. -1.]   3.0
[ 1. -4.  1.  5.]   9.0
[ 2. -2.  1.  3.]   7.0

Dividing row 0 by 2.0
[1.  0.5 2.  0.5]   -2.0
[ 3.  4. -1. -1.]   3.0
[ 1. -4.  1.  5.]   9.0
[ 2. -2.  1.  3.]   7.0

Subtracting row 0 from the lower rows
[1.  0.5 2.  0.5]   -2.0
[ 0.   2.5 -7.  -2.5]   9.0
[ 0.  -4.5 -1.   4.5]   11.0
[ 0. -3. -3.  2.]   11.0

Dividing row 1 by 2.5
[1.  0.5 2.  0.5]   -2.0
[ 0.   1.  -2.8 -1. ]   3.6
[ 0.  -4.5 -1.   4.5]   11.0
[ 0. -3. -3.  2.]   11.0

Subtracting row 1 from the lower rows
[1.  0.5 2.  0.5]   -2.0
[ 0.   1.  -2.8 -1. ]   3.6
[  0.    0.  -13.6   0. ]   27.2
[  0.    0.  -11.4  -1. ]   21.8

Dividing row 2 by -13.6
[1.  0.5 2.  0.5]   -2.0
[ 0.   1.  -2.8 -1. ]   3.6
[-0. -0.  1. -0.]   -2.0
[  0.    0.  -11.4  -1. ]   21.8

Subtracting row 2 from the lower rows
[1.  0.5 2.  0.5]   -2.0
[ 0.   1.  -2.8 -1. ]   3.6
[-0. -0.  1. -0.]   -2.0
[ 0.  0.  0. -1.]   -0.9999999999999964

Dividing row 3 by -1.0
[1.  0.5 2.  0.5]   -2.0
[ 0.   1.  [...]<jupyter_text>Exercise 2 - LU decompositionSolve the following set of equations: $$ \left(\begin{array}{cccc} 2 & 1 & 4 & 1 \\ 3 & 4 & -1 & -1 \\ 1 & -4 & 1 & 5 \\ 2 & -2 & 1 & 3 \end{array} \right) \left(\begin{array}{c} w \\ x \\ y \\ z\end{array} \right) = \left(\begin{array}{c} -4 \\ 3 \\ 9 \\ 7\end{array} \right)$$ with LU decomposition. The LU decomposition for the matrix is $$	\mathbf{L}=\left(\begin{array}{cccc}2 & 0 & 0 & 0 \\ 3 & 2.5 & 0 & 0 \\ 1 & -4.5 & -13.6 & 0 \\ 2 & -3 & -11.4 & -1  \end{array} \right) \quad	\mathbf{U}=\left(\begin{array}{cccc}1 & 0.5 & 2 & 0.5 \\ 0 & 1 & -2.8 & -1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1  \end{array} \right)$$1. Verify that $\mathbf{LU}$ gives the matrix $\mathbf{A}$. You can use the $\texttt{numpy}$ routine $\texttt{matmul}$.2. Perform the double back substitution $\mathbf{Ly}=\mathbf{v}$ and $\mathbf{Ux}=\mathbf{y}$ with the $\texttt{numpy.linalg}$ routine $\texttt{solve}$.3. Verify your result.4. Apply the LU decomposition to the new vectors $\mathbf{v}_1=(1,0,0,0)$, $\mathbf{v}_2=(0,1,0,0)$, $\mathbf{v}_3=(0,0,1,0)$, and $\mathbf{v}_4=(0,0,0,1)$.5. Check your result with $\texttt{solve}$. Talking points1. What do you observe?2. What happens when you apply the LU decomposition to the vectors $\mathbf{v}_1$ to $\mathbf{v}_4$?  Exercise 3 - Linear spring modelSolve the linear masses connected by springs model $$\left( \begin{array}{cccccc} (\alpha-k) & -k \\ -k & \alpha & -k \\ & -k & \alpha & -k \\ & & \ddots & \ddots & \ddots & \\ & & & -k & \alpha & -k \\ & & & & -k & (\alpha-k) \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ \vdots \\ x_{N-1} \\ x_N \end{array} \right) = \left( \begin{array}{c}  C \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{array} \right) $$ where $\alpha=2k-m\omega^2$ for 26 masses with $C=1$, $m=1$, $k=2$ and $\omega=2$. Below is a skeleton program that calls the subroutine $\texttt{banded}$ for the solution of a banded matrix.1. Copy the file $\texttt{banded.py}$ from MyCourses to your directory.2. Complete the program by initializing the problem for the spring model.3. Add plot statements to plot your results in a graph. Talking points1. What do you observe?2. What can you say about the amplitudes of the vibrating masses?<jupyter_code>from numpy import empty,zeros
from banded import banded # the description of the banded subroutine is given below or at the top of the file banded.py
import matplotlib.pyplot as plt

# Constants

# SET THE PARAMETERS OF THE SPRING MODEL FOR 26 MASSES, C=1, m=1, k=6 and omega=2

# Set up the initial values of the arrays

# BUILD THE MATRIX A AND THE VECTOR V

# Solve the equations - call to subroutine banded that solves a banded matrix equation
x = banded(A,v,1,1)

# Make a plot using both dots and lines

# PLOT YOUR RESULTS

######################################################################
#
# Function to solve a banded system of linear equations using
# Gaussian elimination and backsubstitution
#
# x = banded(A,v,up,down)
#
# This function returns the vector solution x of the equation A.x = v,
# where v is an array representing a vector of N elements, either real
# or complex, and A is an N by N banded matrix with "up" nonzero
# elements above the diagonal and "down" nonzero elements below the
# diagonal.  The matrix is specified as a two-dimensional array of
# (1+up+down) by N elements with the diagonals of the original matrix
# along its rows, thus:
#
#   (  -   -  A02 A13 A24 ...
#   (  -  A01 A12 A23 A34 ...
#   ( A00 A11 A22 A33 A44 ...
#   ( A10 A21 A32 A43 A54 ...
#   ( A20 A31 A42 A53 A64 ...
#
# Elements represented by dashes are ignored -- it doesn't matter what
# these elements contain.  The size of the system is taken from the
# size of the vector v.  If the matrix A is larger than NxN then the
# extras are ignored.  If it is smaller, the program will produce an
# error.<jupyter_output>/home/hlappal/projects/compphys/2021/Lecture 3 - Linear and Nonlinear Equations/banded.py:51: RuntimeWarning: divide by zero encountered in double_scalars
  v[m] /= div
/home/hlappal/projects/compphys/2021/Lecture 3 - Linear and Nonlinear Equations/banded.py:54: RuntimeWarning: invalid value encountered in double_scalars
  v[m+k] -= A[up+k,m]*v[m]
/home/hlappal/projects/compphys/2021/Lecture 3 - Linear and Nonlinear Equations/banded.py:60: RuntimeWarning: divide by zero encountered in double_scalars
  A[i,j] /= div
/home/hlappal/projects/compphys/2021/Lecture 3 - Linear and Nonlinear Equations/banded.py:62: RuntimeWarning: invalid value encountered in double_scalars
  A[i+k,j] -= A[up+k,m]*A[i,j]<jupyter_text>Exercise 4 - QR decompositionThe matrix $$\mathbf{A}=\left(\begin{array}{ccc}2 & 6 & 7  \\ 1 & 0 & -1  \\ 2 & 3 & -2    \end{array} \right)$$ has the following QR decomposition: $$\begin{equation*}	\mathbf{Q}=\frac{1}{3}\left(\begin{array}{cccc}2 & 2 & 1  \\ 1 & -2 & 2  \\ 2 & -1 & -2   \end{array} \right) \quad	\mathbf{R}=3\left(\begin{array}{cccc}1 & 2 & 1  \\ 0 & 1 & 2  \\ 0 & 0 & 1   \end{array} \right)\end{equation*} $$1. Verify that $\mathbf{Q}^T\mathbf{Q}=\mathbb{1}$.2. Check that $\mathbf{QR}=\mathbf{A}$. Talking points1. What do you observe?  Exercise 5 - Eigenvalues and eigenvectorsDiagonalize $$\mathbf{A}=\left(\begin{array}{ccc}2 & 1 & 2  \\ 1 & 0 & 3  \\ 2 & 3 & -2    \end{array} \right)$$1. Calculate the eigenvalues and eigenvectors of $\mathbf{A}$ using the $\texttt{numpy}$ function $\texttt{eigh}$.2. Verify that $\mathbf{V}^T\mathbf{AV}=\mathbb{1}$. Talking points1. What do you observe?<jupyter_code>import numpy as np

A = array([[2,6,7],[1,0,-1],[2,3,-2]],float)
Q = array([[2,2,1],[1,-2,2],[2,-1,-2]],float)/3
R = array([[1,2,1],[0,1,2],[0,0,1]],float)*3
np.matmul(Q.T,Q)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18245944039501527, 'fraction_numerical': 0.06689395720667764, 'mean_word_length': 2.6810904370402424, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '143583', 'license': 'MIT', 'n_tokens_mistral': 3561, 'n_tokens_neox': 3283, 'n_words': 1108, 'path': 'Lecture 3 - Linear and Nonlinear Equations/CompPhys_L3_1_in_class_exercises.ipynb', 'repo_name': 'hlappal/comp-phys', 'size': '8506'}