Due on September 26 at 11:59pm electronically via email.
Use a word-processor like Microsoft Word or LibreOffice (free), or \( \LaTeX \) (also free) to create a single PDF file with all of your answers. If you want to include code, do not put them in this file. This file should include your graphs, equations, tables, and explanations. The name of this file should have the following format: LASTNAME_STUDENTNUMBER_A1.pdf. You can hand-write your answers, but you will have to scan your papers, or take a picture of them. In any case, send me one PDF file. I will not accept paper submissions.
If you want to include multiple code files (.m or .py files), put them all in one folder, zip the folder, and send me only one zip file. The name of this file should look like this: LASTNAME_STUDENTNUMBER_A1.zip
Attach both files to an email and send it to msamani[at]fields.utoronto.ca. Make sure the title of the email starts with MATH3090.
Question 1: Consider the number \(x=1.43\).
Questions from the textbook: 2.3, 2.5, 2.9, 2.10, 2.19.
The fucntions in part (a) are as follows:
import numpy as np import scipy.linalg # I could not find an implementation of magic in python, so we have to define it ourselves. # Code was adopted from https://stackoverflow.com/questions/47834140/numpy-equivalent-of-matlabs-magic and modified on Sep 10, 2018. def magic(n): n = int(n) if n < 3: raise ValueError("Size must be at least 3") if n % 2 == 1: p = np.arange(1, n + 1) return n * np.mod(p[:, None] + p - (n + 3) // 2, n) + np.mod(p[:, None] + 2 * p - 2, n) + 1 elif n % 4 == 0: J = np.mod(np.arange(1, n + 1), 4) // 2 K = J[:, None] == J M = np.arange(1, n * n + 1, n)[:, None] + np.arange(n) M[K] = n * n + 1 - M[K] else: p = n // 2 M = magic(p) M = np.block([[M, M + 2 * p * p], [M + 3 * p * p, M + p * p]]) i = np.arange(p) k = (n - 2) // 4 j = np.concatenate((np.arange(k), np.arange(n - k + 1, n))) M[np.ix_(np.concatenate((i, i + p)), j) ] = M[np.ix_(np.concatenate((i + p, i)), j)] M[np.ix_([k, k + p], [0, k])] = M[np.ix_([k + p, k], [0, k])] return M # scipy.linalg has a built-in command for Hilbert matrices H = scipy.linalg.hilbert(n) # this will make a Hilbert matrix of size 3 # scipy.linalg also has a built in package for Pascal matrices P = scipy.linalg.pascal(n) # numpy has an equivalent eye command E = np.eye(n, n)
I have translated the functions defined in section 2.7 into Python:
# Assumes L is a lower-triangular square matrix def forward(L, b): N = len(L) x = np.zeros(N) for k in range(N): j = slice(k) x[k] = (b[k] - L[k,j].dot(x[j]) ) / L[k,k] return x # Assumes U is an upper-triangular square matrix def backsubs(U, b): N = len(U) x = np.zeros(N) for k in range(N-1, -1, -1): j = slice(k,N) x[k] = (b[k] - U[k,j].dot(x[j])) / U[k,k] return x # LU factorization def LUtx(AA): A = np.copy(AA) N = len(A) p = np.array(range(N)).T for k in range(N-1): # Fint the largest element below diagonal in k-th column m = np.argmax(np.abs(A[k:N,k])) + k # Skip elimination if column is zero if A[m,k] != 0: # Swap pivot now if m != k: A[[k,m],:] = A[[m,k],:] p[[k,m]] = p[[m,k]] # Compute multipliers i = slice(k+1,N) A[i,k] = A[i,k]/A[k,k] # Update the remainder of the matrix j = slice(k+1,N) A[i,j] = A[i,j] - A[i,k:k+1] * A[k:k+1,j] L = np.tril(A, k=-1) + np.eye(N) U = np.triu(A, k=0) return L, U, p # It solves A.dot(x) = b and returns x def bslashtx(A,b): N = len(A) if np.all(np.triu(A, k=1) == np.zeros((N,N))): # Lower triangular return forward(A, b) elif np.all(np.tril(A, k=-1) == np.zeros((N,N))): # Upper triangular return backsubs(A, b) elif np.all(A == A.T): # symmetric L = np.linalg.cholesky(C) y = forward(L, b) x = backsubs(L.T,y) else: L,U,p = LUtx(A) y = forward(L, b[p]) x = backsubs(U,y) return x
Question 2.19 The sparse matrix libraries are in scipy.sparse