2018-02-02.ipynb

Kernel: SageMath 8.1

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

Homework 3: due February 2, 2018

Please enter all answers within this notebook unless otherwise specified. As usual, don't forget to cite sources and collaborators.

Through this problem set, use the SageMath 8.1 kernel unless otherwise specified.

Problem 1: Gradient vector fields

Grading criterion: correctness of code.

1a. Compute the gradient of $f(x,y) = 3\sin(x) - 2\cos(2y) - x - y$.

1b. Plot the 2-dimensional vector field defined by the gradient of $f$ in the rectangle $(-2,-2) \leq (x,y) \leq (2,2)$.

Problem 2: 3D plotting

Grading criterion: correctness of code.

2a. Draw a 3D plot of a torus.

2b. Draw a single 3D plot containing the five regular polytopes in it: tetrahedron, cube, octahedron, dodecahedron, icosahedron. All five must be visible.

2c. Draw a 3D plot of the Mexican hat function. Try to choose the parameter $\sigma$ to get an authentic "sombrero" look.

Problem 3: MATLAB (and Octave) vs. Sage

Grading criterion: correctness of code.

This exercise refers to the Math 18 MATLAB exercise set.

3a. Do MATLAB exercise 4.5 twice: once using Octave, and a second time using Sage. (For this problem, you will need to switch the kernel between Octave and SageMath.)

your Octave answer goes here

your SageMath answer goes here

3b. Repeat with MATLAB exercise 5.6, using numpy to obtain the analogue of MATLAB's least squares functionality.

your Octave answer goes here

your SageMath answer goes here

Problem 4: Eigenvectors and eigenvalues

Grading criteria: correctness of code and explanations.

Let $M$ be the following matrix with rational entries: $$M = \left(\begin{array}{rrrr} \frac{1}{2} & 0 & -1 & 0 \\ 1 & \frac{1}{2} & 1 & 1 \\ 0 & 1 & -1 & -1 \\ 2 & 1 & -2 & 1 \end{array}\right)$$

4a. State the Cayley-Hamilton theorem, then use Sage to verify that it holds for $M$.

4b. Compute the eigenvalues and eigenvectors of $M$, and verify (using numerical approximations) that each eigenvector is indeed an eigenvector with the specified eigenvalue.

4c. State the relationship between the characteristic polynomial, the determinant, and the eigenvalues of a matrix; then verify numerically that this holds for $M$.

Problem 5: Hilbert matrices and numerical stability

Grading criterion: correctness of code and thoroughness of explanation.

5a. Define a Python function f, which on the input of a positive integer $n$, returns the $n \times n$ Sage matrix over the rational numbers $$H_{ij} = \frac{1}{i+j-1} \qquad (i,j=1,\dots,n).$$ See below for a sample output. Hint: remember that Sage, like Python, starts indexing with 0 instead of 1.

def f(n):
    # your code goes here

  File "<ipython-input-29-a3ab5dd50729>", line 2
    # your code goes here
                         ^
IndentationError: expected an indented block

## Sample output of f(3):
[  1 1/2 1/3]
[1/2 1/3 1/4]
[1/3 1/4 1/5]

5b. Write Sage code to compute the inverse of f(25), print out the top left entry (not the whole matrix), and verify that the whole answer agrees with the formula in Wikipedia.

5c. Use the change_ring method to redo the computation of the inverse of f(25) over the ring RR (i.e., using floating-point real numbers with double precision == 53 bits), agian printing out the top left entry.

5d. In as much detail as you can, explain what you have just observed. You may want to compute the determinant of f(25) and/or read about numerical stability.

Problem 6: Linear algebra over Q

Grading criteria: correctness of code and thoroughness of explanation.

6a. Explain in detail what the following code does and how it is being done. Assume that the input M is a matrix with entries in QQ.

def f(M):
    M = copy(M) # replace M with a copy to avoid clobbering the original
    n = M.dimensions()[0]
    if n != M.dimensions()[1]: raise ValueError
    for i in range(n):
        j = i
        while j<n and M[j,i] == 0: j += 1
        if j==n: return(0)
        if i<j:
            M.swap_rows(i,j)
            M.rescale_row(i, -1)
        for k in range(i+1,n):
            M.add_multiple_of_row(k,i,-M[k,i]/M[i,i])
    return(prod(M[i,i] for i in range(n)))

6b. Define a Python function g that, on input a positive integer n, returns an $n \times n$ matrix in which each entry is a Sage integer is chosen uniformly at random among 0 and 1. (Hint: use the randint function, but read the docstring carefully!)

Then check that on the output of g(25), the function f agrees with the native Sage way of doing the same computation.

def g(n):
    # your code goes here

6c. Define the binary height of a nonzero rational number $\frac{r}{s}$, written in lowest terms, to be the quantity $$\left\lfloor\log_2 r \right\rfloor + 1 + \left\lfloor \log_2 s \right\rfloor + 1.$$ Use the function binary_height defined below to compute this.

Take the code defining f, recopy it below, then modify it so that at the last step, instead of returning the product of the diagonal entries of M, it returns the maximum of the heights of the diagonal entries.

def binary_height(x):
    return(x.numer().nbits() + x.denom().nbits())

def f_modified(M):
    # your modified code goes here

6d. Explain, in words, the output of the following code. (It should run without errors!)

M = g(25).change_ring(QQ)
print(binary_height(f(M)))
print(f_modified(M))

6e. Sage does not use the method illustrated by the definition of f for matrices of integers or rational numbers; instead, it uses an approach based on the Chinese remainder theorem, in which one does the analogous computation modulo $p$ for a few primes $p$. Explain why you think this decision was made.