Final Exam Review¶
Reference Books¶
- LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
- L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
- TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.
Topics¶
The final exam will focus on material not tested on the midterm. However, this class is cumulative, and it is necessary to understand earlier material as well. In particular, it is expected that students have developed a working knowledge of Python.
The following is a list of topics that may be explicitly tested on the final exam.
Numerical Differentiation following TAK § 3.2. Notebook: Numerical Differentiation [download] [view] [azure] |
Numerical Integration following TAK § 3.3. |
Numerical Integration following TAK § 3.4. |
Matrices and vectors in Numpy following LL § 2.3. |
Gaussian Elimination following TAK § 4.2. |
LU Factorization and Differential Equations (TAK § 4.3). |
Least Squares (TAK § 4.5). |
Eigenvalues (TAK § 4.6). Notebook: Eigenvalues and Eigenvectors [download] [view] [azure] |
Dictionaries (L § 6.2). |
Begin discussing classes following L Chapter 7. |
Random Numbers following L Chapter 8. |
Monte Carlo Integration following L Chapter 8. |
Sources of review¶
- The homework assignments and solutions.
- The textbooks listed above and exercises on the topics above.
- The notebooks available on the Class' Calendar page
- Class videos, available from Blackboard Collaborate.
Problems¶
Below some problems are included that are aimed at reviewing the topics above. They do not necessarily represent problems that would be assigned on the final exam.
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
from scipy import linalg
from numpy import random
import random as random_number
1. Differentiation:¶
Write a funtion derivative(f, h) which takes as input:
- A function
f(x), which takes as input a floating point real number and produces a floating point real number as its output. - A float $h>0$.
The derivative function should return a new function df(x). The function df(x) should take as input a floating point real number x and return the symmetric two-point approximation of the derivative of $f$ at $x$ using points at distance $h$ from $x$.
The derivative function should make use of currying.
def derivative(f, h):
def df(x):
return (f(x+h) - f(x-h)) / (2*h)
return df
# Test
f = lambda x: np.sin(x)
df = derivative(f, 10**-8)
for x in np.linspace(0, 3, 201):
assert abs( df(x) - np.cos(x) ) < 10**-8, \
"The value df({:0.2f}) is incorrect.".format(x)
2. Speed of a crash¶
A malfunctioning drone crashes destroying itself. The dictionary altitude below maps time (in seconds) to measurements of a drone's altitude (in meters). The drone crashes shortly after the $5$ second mark.
- Approximately fast was the drone falling when it crashed into the ground? Store the answer (which would be in meters per second) in the variable
crash_speed. - Before the drone crashed, it accellerated upward for some unknown reason. What was the fastest speed it was moving upward over the provided time interval? At approximately what time was it moving at that speed? Store the respective answers in
up_speedandup_time.
altitude = {0.00: 2.0438, 0.10: 2.1013, 0.20: 2.1577, 0.30: 2.2136, 0.40: 2.2699,
0.50: 2.3271, 0.60: 2.3859, 0.70: 2.4472, 0.80: 2.5117, 0.90: 2.5799,
1.00: 2.6528, 1.10: 2.7311, 1.20: 2.8154, 1.30: 2.9066, 1.40: 3.0053,
1.50: 3.1123, 1.60: 3.2281, 1.70: 3.3535, 1.80: 3.4888, 1.90: 3.6347,
2.00: 3.7915, 2.10: 3.9593, 2.20: 4.1384, 2.30: 4.3285, 2.40: 4.5294,
2.50: 4.7406, 2.60: 4.9613, 2.70: 5.1904, 2.80: 5.4263, 2.90: 5.6672,
3.00: 5.9109, 3.10: 6.1543, 3.20: 6.3941, 3.30: 6.6263, 3.40: 6.8462,
3.50: 7.0480, 3.60: 7.2256, 3.70: 7.3716, 3.80: 7.4775, 3.90: 7.5341,
4.00: 7.5307, 4.10: 7.4553, 4.20: 7.2945, 4.30: 7.0336, 4.40: 6.6560,
4.50: 6.1436, 4.60: 5.4763, 4.70: 4.6323, 4.80: 3.5874, 4.90: 2.3156,
5.00: 0.7886, }
First to organize the input times. (Dictionaries are fickle, and floating point arithmetic is not exact. Probably floats are not great inputs for dictionaries.)
times = []
for t in altitude:
times.append(t)
times
crash_speed = (altitude[5] - altitude[4.9]) / 0.1
crash_speed
up_speed = 0
up_time = 0
for i in range(len(times)-1):
t = times[i]
next_t = times[i+1]
s = (altitude[next_t] - altitude[t]) / (0.1)
if s > up_speed:
up_speed = s
up_time = (t + next_t)/2
print("up_time = {}".format(up_time))
print("up_speed = {}".format(up_speed))
3. Weird Quadrature rule¶
There was an error in the question below which was now corrected.
This problem concerns the following quadrature rule, defined for any interval $[a, b]$: $$Q_{[a,b]}(f) = (b-a) \left(\frac{1}{4} f\left(\frac{2a+b}{3}\right) + \frac{1}{2} f\left(\frac{a+b}{2}\right) + \frac{1}{4} f\left(\frac{a+2b}{3}\right)\right).$$
Define a function
q(a , b, f)which takes as input floating point real numbers $a$ and $b$ with $a<b$ and a functionfrepresenting a real-valued function $f:{\mathbb R} \to {\mathbb R}$. The function should return $Q_{[a,b]}(f)$.What is the degree of of precision of $Q_{[a,b]}$? Store the answer in the variable
q_degree.
def q(a, b, f):
return (b-a) * (0.25 * f((2*a+b)/3) + 0.5 * f((a+b)/2) + 0.25 * f((a+2*b)/3))
# Just a check
f = lambda x: x**2
print(q(0, 1, f)) # computed value
print(0.25*(1/3)**2 + 0.5*(1/2)**2+0.25*(2/3)**2) # by hand
The above looks okay. Lets check the degree of precision. We'll use the interval $[0,1]$.
q(0, 1, lambda x: 1)
The above answer of $1.0$ agrees with $\int_0^1 1~dx$, so the degree of precision is at least $0$.
q(0, 1, lambda x: x)
The above answer of $0.5$ agrees with $\int_0^1 x~dx$, so the degree of precision is at least $1$.
q(0, 1, lambda x: x**2)
The above answer of $0.26\ldots$ is not the same as $\int_0^1 x^2~dx=\frac{1}{3}$, so the degree of precision is exactly $1$.
q_degree = 1
q_degree
4. Compound rule¶
Write a function q_compound(n, f) which takes as input a positive integer $n$ and a function $f$ as above. The function should divide the interval $[0,1]$ into $n$ equal sized subintervals. It should return the sum of $Q_{\ast}(f)$ where $\ast$ varies over the $n$ subintervals.
Note that we can divide $[0,1]$ into $n$ equal sized intervals with the code
np.linspace(0, 1, n+1)
This is demonstrated below.
n = 3
np.linspace(0, 1, 3+1)
def q_compound(n, f):
total = 0.0
x = np.linspace(0, 1, n+1)
for i in range(n):
a = x[i]
b = x[i+1]
total += q(a, b, f)
return total
Recall $\int_0^1 \sin x~dx=1-\cos(1).$ We'll check that our q_compound is about right.
val = q_compound(100, lambda x: np.sin(x))
print("val = {} should be about 1-cos(1) = {}.".format(val, 1-np.cos(1)))
5. Matrix Manipulation 1¶
Write a function manipulate1(m) which takes as input a $2$-dimensional numpy array m representing a matrix $M$. The function should make a copy of m (so that m is not changed). Then the following operations should be performed, in order:
- All diagonal entries should be doubled.
- The top right and bottom left entries should be swapped.
- The entries which are strictly closer to the bottom right courner than any other corner should be negated.
For example, the if the following matrix is input for manipulate(m)
$$\left(\begin{matrix}
2 & 3 & 4 \\
1 & -1 & 3 \\
8 & 5 & 11
\end{matrix}\right)$$
then the following matrix should be output:
$$\left(\begin{matrix}
4 & 3 & 8 \\
1 & -2 & 3 \\
4 & 5 & -22
\end{matrix}\right).$$
def manipulate1(m):
mat = m.copy()
rows,cols = mat.shape
# double all diagonal entries
for i in range( min(rows, cols) ):
mat[i,i] *= 2
# Swap top right and bottom left entry
top_right = mat[0, cols-1]
bot_left = mat[rows-1, 0]
mat[0, cols-1] = bot_left
mat[rows-1, 0] = top_right
# Find entries strictly closest to the bottom right corner, and negate them.
# A row is strictly closer to the last row than to the first row if
# (rows-1)-row < row
# A column is strictly closer to the last column than the first if
# (cols-1)-col < col
for row in range(rows):
for col in range(cols):
if (rows-1)-row < row and (cols-1)-col < col:
mat[row,col] *= -1
return mat
# Test suggested in the problem
m_in = np.array( [ [ 2, 3, 4 ],
[ 1, -1, 3 ],
[ 8, 5, 11] ] )
m_out = np.array( [ [ 4, 3, 8 ],
[ 1, -2, 3 ],
[ 4, 5, -22] ] )
assert (manipulate1(m_in) == m_out).all()
# Another test suggested in the problem
m_in = np.array( [ [ 2, 3, 4 ],
[ 1, -1, 3 ],
[ 1, -1, 3 ],
[ 1, -1, 3 ] ] )
# Done by hand:
m_out = np.array( [ [ 4, 3, 1 ],
[ 1, -2, 3 ],
[ 1, -1, -6 ],
[ 4, -1, -3 ] ] )
assert (manipulate1(m_in) == m_out).all()
6. Matrix Manipulation 2¶
Write a function manipulate2(m) which takes as input a $2$-dimensional numpy array m representing a matrix $M$. The function should return the matrix
$$\left(\begin{array}{r|r}
4I & 3 M \\
\hline
M^T & M^T M
\end{array}\right).$$
Here, $M^T$ represents the transpose of $M$. The matrix $I$ represents the identity matrix whose sizes are chosen to make the above block matrix form well-defined. The matrix should be returned as a $2$-dimensional numpy array.
This can be done with a block matrix command in numpy or done by hand. We'll denonstrate both. Here is a test matrix:
m = np.array( [ [ 1, 2 ],
[ 3, 4 ],
[ 5, 6 ] ] )
Typing np.blo<tab> shows that np.block() is a method of making a numpy array. We'll just try it.
out = np.block( [ [ 4*np.identity(3), 3*m ],
[ m.transpose(), m.transpose() @ m ] ] )
out
Here we make it a function.
def manipulate2(m):
rows, cols = m.shape
return np.block( [ [ 4*np.identity(rows), 3*m ],
[ m.transpose(), m.transpose() @ m ] ] )
# Test using the above example
assert (out == manipulate2(m)).all()
# done by hand
def manipulate2(m):
rows, cols = m.shape
r = np.zeros((rows+cols, rows+cols))
# fill in 4*I
for i in range(rows):
r[i, i] = 4
# fill in 3*m
for row in range(rows):
for col in range(cols):
r[row, rows+col] = 3*m[row, col]
# fill in m.transpose
for row in range(rows):
for col in range(cols):
r[rows+col, row] = m[row, col]
# compute M^T * M, which is a cols x cols matrix:
prod = m.transpose() @ m
for row in range(cols):
for col in range(cols):
r[rows+row,rows+col] = prod[row,col]
return r
# Test using the above example
assert (out == manipulate2(m)).all()
7. Column version of Gaussian elimination¶
Write a function elimination(m) which takes as input a $2$-dimensional numpy array m of floats representing a matrix $M$ with at least as many columns as rows (added some additional requirements for m). The function should make a copy of m so that this matrix is unchanged. The function should then perform operations of adding multiples of one column to another until the matrix has an upper triangular form. Other operations are not permissible. The upper triangular matrix should be returned.
# As an example 1 take:
m1 = np.array([[ 1, 2, 3],
[ 4, 5, 6]], dtype = float)
# We'd want to subtract 4/5 times column 1 to column 0.
# This can be done with the following:
m1[:,0] -= 4/5 * m1[:,1]
m1
Zeroing out entries to the left of the diagonal will work well unless the diagonal entry is zero. Assuming a diagonal entry is zero, we can make it non-zero (typically) by adding a column to the left. (If there are only zero entries to the left, then we don't need to do anything.) So, here is a slightly more challenging example.
# As a second example take:
m2 = np.array([[ 1, 2, 3],
[ 4, 0, 6]], dtype = float)
# Make entry in position (1, 1) non-zero
m2[:,1] += m2[:,0]
m2
# Now zero out:
m2[:,0] -= 4/4 * m2[:,1]
m2
# As a bigger example take:
m3 = np.array([[ 1, 2, 3],
[ 4, 0, 6],
[ 7, 3, 0]], dtype = float)
# First we'll work on row 1.
# Add column 0 to column 1.
m3[:,1] += m3[:,0]
m3
# Zero out the entry to the left of (1,1)
m3[:,0] -= m3[:,1]
m3
# Now we work on row 2.
# Add column 1 to column 2.
m3[:,2] += m3[:,1]
m3
# Zero out entry (2, 0)
m3[:,0] -= -3/10 * m3[:,2]
m3
Oops we messed up entry (1, 0). This means we need to work from the bottom row to the top row instead.
def elimination(m):
m = m.copy()
rows, cols = m.shape
assert cols >= rows
for i in range(rows-1, 0, -1): # Move through the diagonal entries.
zero_out_left = True # Will turn to false if we don't need to zero
# out entries to the left of position (i,i)
if m[i,i] == 0:
# Find a non-zero entry to the left
j = np.argmax(abs(m[i,:i]))
if m[i,j] == 0:
# All entries to the left are already zero
zero_out_left = False
else:
# We'll add column j to column i
m[:,i] += m[:,j]
if zero_out_left:
for j in range(i):
m[:,j] -= m[i,j]/m[i,i] * m[:,i]
return m
# First example:
m = np.array([[ 1, 2, 3],
[ 4, 5, 6]], dtype = float)
elimination(m)
# Second example:
m = np.array([[ 1, 2, 3],
[ 4, 0, 6]], dtype = float)
elimination(m)
# As a bigger example take:
m3 = np.array([[ 1, 2, 3],
[ 4, 0, 6],
[ 7, 3, 0]], dtype = float)
elimination(m3)
# Random 4x4 example
# found with np.random.random_sample(16).reshape((4,4))
m4 = np.array( [ [ 0.09978885, 0.40719259, 0.65440626, 0.44289351],
[ 0.57510288, 0.24356203, 0.05080939, 0.4282749 ],
[ 0.19026836, 0.66561491, 0.6450414, 0.01432258],
[ 0.49572319, 0.3316222, 0.61250047, 0.45688589]] )
print(m4)
print(elimination(m4))
Notice that in the above example, entries are not exactly zero. We can fix that by defining the entry that we are zeroing out to be exactly zero.
def elimination(m):
m = m.copy()
rows, cols = m.shape
assert cols >= rows
for i in range(rows-1, 0, -1): # Move through the diagonal entries.
zero_out_left = True # Will turn to false if we don't need to zero
# out entries to the left of position (i,i)
if m[i,i] == 0:
# Find a non-zero entry to the left
j = np.argmax(abs(m[i,:i]))
if m[i,j] == 0:
# All entries to the left are already zero
zero_out_left = False
else:
# We'll add column j to column i
m[:,i] += m[:,j]
if zero_out_left:
for j in range(i):
m[:,j] -= m[i,j]/m[i,i] * m[:,i]
m[i,j] = 0.0
return m
This new version fixes the above:
print(m4)
print(elimination(m4))
8. Differential Equations¶
Consider the differential equation $$y''(x) + \sin(x) y(x) = x^2$$ on the interval $[0, \pi]$ subject to the boundary conditions $y(0)=1$ and $y(\pi)=0$.
Given an integer $n \geq 3$, consider cutting the interval into $n$ equal-sized subintervals. You are interested in approximating the solution $y$ at the cut points $$x_0 = \frac{\pi}{n}, x_1 = \frac{\pi}{n}, \ldots, x_{n-1}=\frac{(n-1)\pi}{n}.$$ Set $y_i = y(x_i)$. Using numerical derivatives, find an $n-1 \times n-1$ matrix $M$ and a vector $v \in {\mathbb R}^{n-1}$ such that solving $M y= v$ leads to an approximate solution $y=(y_0, \ldots, y_{n-1})$ for the differential equation.
Write a function linearization(n) which returns a pair of numpy arrays (m,v). Here m should represent the matrix $M$ from the previous paragraph and v should reprensent the vector $v$.
The $x$ values can be obtained from $n$ with:
n = 4
x = np.linspace(np.pi/n, (n-1)/n * np.pi, n-1)
x
Recall that the second derivative of y is approximately $$y''(x) = \frac{y(x+h) - 2(x) + y(x-h)}{h^2}.$$ Setting $y_i=y(x_i)$ and $h=\frac{\pi}{n}$, this becomes $$y''(x_i) \approx \frac{n^2}{\pi^2} (y_{i+1} - 2 y_i + y_{i-1}).$$ We take $x_{-1}=0$ and $y_{-1}=y(0)=1$ and take $x_n=\pi$ and $y_n=y(\pi)=0$, but we'll deal with these special cases later.
The left side of our diferential equation is $y''(x)+\sin(x) y(x)$ which for $x=x_i$ becomes approximately $$\frac{n^2}{\pi^2} (y_{i+1} - 2 y_i + y_{i-1}) + \sin(x_i)y_i = \frac{n^2}{\pi^2} y_{i+1} + \big(sin(x_i)-2\frac{n^2}{\pi^2}\big) y_i + \frac{n^2}{\pi^2}y_{i-1}.$$
The right side of our differential equation is $x^2$ which at $x=x_i$ simply becomes $x_i^2$. So this gives rise to the linear equation: $$\frac{n^2}{\pi^2}y_{i+1} + \big(sin(x_i)-2\frac{n^2}{\pi^2}\big) y_i + \frac{n^2}{\pi^2}y_{i-1}=x_i^2.$$
Now we turn our attention to the special cases. When $i=0$, we have $x_{-1}=0$ and $y_{-1}=1$ and thus the equation becomes $$\frac{n^2}{\pi^2} y_{1} + \big(sin(x_0)-2\frac{n^2}{\pi^2}\big) y_0 + \frac{n^2}{\pi^2} y_{-1}=x_0^2$$ or $$\frac{n^2}{\pi^2} y_{1} + \big(sin(x_0)-2\frac{n^2}{\pi^2}\big) y_0 = x_0^2 - \frac{n^2}{\pi^2}$$ When $i=n-1$, we have to consider that $x_n=\pi$ and $y_n=0$. Thus we get $$\big(sin(x_{n-1})-2\frac{n^2}{\pi^2}\big) y_{n-1} + \frac{n^2}{\pi^2}y_{n-2}=x_{n-1}^2.$$
The matrix form of these linear equations are then: $$\left(\begin{matrix} sin(x_0)-2\frac{n^2}{\pi^2} & \frac{n^2}{\pi^2} & 0 & \ldots \\ \frac{n^2}{\pi^2} & sin(x_1)-2\frac{n^2}{\pi^2} & \frac{n^2}{\pi^2} & 0 & \ldots \\ 0 & \frac{n^2}{\pi^2} & sin(x_2)-2\frac{n^2}{\pi^2} & \frac{n^2}{\pi^2} & 0 & \ldots \\ \vdots & \ddots & \ddots & \ddots & \ddots \\ 0 & \ldots & 0 & \frac{n^2}{\pi^2} & sin(x_{n-2})-2\frac{n^2}{\pi^2} & \frac{n^2}{\pi^2} \\ 0 & \ldots & \ldots & 0 & \frac{n^2}{\pi^2} & sin(x_{n-1})-2\frac{n^2}{\pi^2} \end{matrix}\right) \left(\begin{matrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-2} \\ y_{n-1} \end{matrix}\right)= \left(\begin{matrix} x_0^2-\frac{n^2}{\pi^2} \\ x_1^2 \\ x_2^2 \\ \vdots \\ x_{n-2}^2 \\ x_{n-1}^2 \end{matrix}\right) $$
def linearization(n):
# Build and n-1 x n-1 matrix of all zeros
m = np.zeros((n-1, n-1))
# Build the x list
x = np.linspace(np.pi/n, (n-1)/n * np.pi, n-1)
# handle the diagonal entries.
for i in range(n-1):
m[i,i] = np.sin(x[i]) - 2 * n**2 / np.pi**2
# handle the off diagonal entries:
for i in range(n-2):
m[i,i+1] = n**2 / np.pi**2
m[i+1,i] = n**2 / np.pi**2
# Most of the vector entries are x[i]^2
v = x**2
# We can fix the first entry:
v[0] -= n**2 / np.pi**2
return (m, v)
To check our work, we can actually solve the equation.
n = 1000
x = np.linspace(np.pi/n, (n-1)/n * np.pi, n-1)
m,v = linearization(n)
y = np.linalg.solve(m, v)
plt.plot(x, y)
# Check the differential equation
# Padded version of y, with y_-1=1 and y_n=0 added.
yy = np.array([1.0] + list(y) + [0.0])
# Second derivative:
ypp = (yy[0:n-1] - 2*yy[1:n] + yy[2:n+1]) / (np.pi**2 / n**2)
# Left side minus right side:
# y'' + sin(x) * y - x**2
difference = ypp + np.sin(x)*y - x**2
max(abs(difference))
The above number is small, and indicates the solution works. It indicates that the maximum of $$| y''(x) + \sin(x) y(x) - x^2|$$ over $[0,\pi)$ is small.
9. Orthogonality¶
Consider the $L^2$ inner product on $[0, 1]$ given by
$$\langle f, g \rangle = \int_0^1 f(x)g(x)~dx.$$
This inner product has been implemented using the quad function from scipy.integrate
below.
from scipy.integrate import quad
def ip(f, g):
return quad(lambda x:f(x)*g(x), 0, 1)[0]
Given continuous functions $f$ and $g$ (which are not zero on $[0,1]$), there is always a unique constant $c$ such that $f$ and $g+cf$ are orthogonal with respect to $\langle, \rangle$. Here $g+cf$ denotes the function sending $x$ to $g(x)+c f(x)$.
Write a function orth(f, g) which returns this constant c.
The projection of $g$ onto $\mathit{span} \{f\}$ is $$\frac{\langle g, f\rangle}{\langle f, f\rangle} f.$$ Therefore $$g - \frac{\langle g, f\rangle}{\langle f, f\rangle} f.$$ is perpendicular to $f$. Therefore, we want $c=- \frac{\langle g, f\rangle}{\langle f, f\rangle} f.$
def orth(f, g):
return - ip(g, f) / ip(f, f)
# Check our work with an example.
f = lambda x: np.sin(x) # f(x) = sin x
g = lambda x: x * np.exp(x) # g(x) = x * e^x
c = orth(f, g)
h = lambda x: g(x) + c*f(x)
ip(f, h)
The above number is very small, so f and h are orthogonal and our code worked.
10. Least Squares¶
Write a function projection(f, g, h) which takes as input three functions representing functions ${\mathbb R} \to {\mathbb R}$. The function should return a pair of real numbers $(a,b)$ so that the function $a f + b g$ is the orthogonal projection of $h$ to the span of $f$ and $g$. Orthogonality should be measured with respect to the inner product from problem 9.
To compute the orthogonal projection of $h$ onto $\mathit{span}\{f, g\}$, we first need to compute an orthogonal basis for the span. We can define a replacement for $g$ following the previous problem, setting $$g' = g - \frac{\langle g, f\rangle}{\langle f, f\rangle} f.$$ Then the projection of $h$ onto $\mathit{span}\{f, g\}$ is the sum of the projections on $\mathit{span}\{f\}$ and $\mathit{span}\{g'\}$: $$\frac{\langle h, f\rangle}{\langle f, f\rangle} f+\frac{\langle h, g'\rangle}{\langle g', g'\rangle} g'.$$
def projection(f, g, h):
c = orth(f, g)
gp = lambda x: g(x) + c*f(x) # This is g' as above
coef_f = ip(h, f) / ip(f, f)
coef_gp = ip(h, gp) / ip(gp, gp)
return lambda x: coef_f*f(x) + coef_gp*gp(x)
# Lets test it.
f = lambda x: np.sin(x) # f(x) = sin(x)
g = lambda x: 1 - x**2 # g(x) = 1 - x^2
h = lambda x: np.exp(x) # h(x) = e^x
p = projection(f, g, h)
x = np.linspace(0, 1, 200)
plt.plot(x, h(x), "r") # Plot y=e^x in red
plt.plot(x, p(x)) # Plot the projection in blue
plt.show()
It looks sort of reasonable. I guestion if there is a better linear combination of $f$ and $g$ for fitting $h$. Thats hard to see with just a graph.
11. Eigenvalues and Eigenvectors¶
The matrix $A$ below has a real eigenvalue $\lambda$ which is very close to $-3$.
A = np.array([[ 4, 1, 2, 4],
[-2, 4, -1, 2],
[ 1, -2, 3, 4],
[ 3, 4, 2, 1]
])
We define $B = (A+3I)^{-1}$. This matrix is approximately given below.
B = np.array([[-0.36363636, -0.54545455, -0.27272727, 0.90909091],
[-0.81818182, -0.72727273, -0.36363636, 1.54545455],
[-1.40909091, -1.36363636, -0.43181818, 2.52272727],
[ 1.79545455, 1.81818182, 0.78409091, -3.23863636]])
Use powers of $B$ to compute the eigenvalue lambda_A of $A$ which is closest to $-3$. Store the associated eigenvector in eigenvector_A. Both quantities should be correct to six decimal places.
Suppose $A v = \lambda v$ with $\lambda$ very close to $-3$. Then $(A +3 I)v=(\lambda+3) v$, and $\lambda+3$ will be very close to zero. We will also have $(A +3 I)^{-1} v=\frac{1}{\lambda+3} v$ and $\frac{1}{\lambda+3}$ will be very large. Thinking carefully about this, we are looking for the dominant eigenvalue eigenvector pair for $B$. So, we can use the power method.
x = 2*np.random.random_sample(4) - 1
k = 0
x_k = x / np.linalg.norm(x)
eigenvalue_B = x_k @ B @ x_k
eigenvalue_old = float("+inf") # infinity
while abs(eigenvalue_B - eigenvalue_old) > 10**-8:
temp = B @ x_k
k += 1
x_k = temp / np.linalg.norm(temp)
eigenvalue_old = eigenvalue_B
eigenvalue_B = x_k @ B @ x_k
print("we found an eigenvector {} with eigenvalue {}.". \
format(x_k, eigenvalue_B))
# Check eigenvalue equation
B @ x_k - eigenvalue_B*x_k
That looks good. Now we have $$\frac{1}{\lambda_A+3} = \lambda_B,$$ Thus $\lambda_A = \frac{1}{\lambda_B} - 3.$$
lambda_A = 1/eigenvalue_B - 3
eigenvector_A = x_k
# Check
A @ eigenvector_A - lambda_A * eigenvector_A
Looks good!
12. The last eigenvalue¶
The matrix $C$ below has four distinct eigenvalues. Three of them are given below (correct to several digits of accuracy). Find the fourth and store the result in eigvalue_3 to approximately the same accuracy.
C = np.array([[ 5, 1, 2, 4],
[-2, 4, -1, 2],
[ 1, -2, 3, 4],
[ 3, 4, 2, 10]
])
eigvalue_0 = 12.92391499
eigvalue_1 = 1.23951179
eigvalue_2 = 1.60202972
The sum of the eigenvalues is the trace, and the trace is the sum of the diagonal entries.
trace = 0.0
for i in range(4):
trace += C[i,i]
trace
eigvalue_3 = trace - eigvalue_0 - eigvalue_1 - eigvalue_2
eigvalue_3
12. Estimating the area of the deltoid¶
The deltoid is the region in the plane satisfying $$(x^2+y^2)^2 + 18(x^2+y^2)-27 < 8(x^3-3xy^2).$$ The deltoid is contained in the rectangle $$R=[-2, 3] \times [-3, 3] = \{(x,y):~-2 \leq x \leq 3 \text{ and } -3 \leq y \leq 3\}.$$
Write a function deltoid_area(n) which takes as input a positive integer $n$. The function should choose $n$ uniformly random points from $R$ and check if they are contained in the deltoid. The number of points in the deltoid should be used to estimate its area using Monte Carlo Integration, and this area estimate should be returned.
def deltoid(x, y):
"""Return True if (x,y) is in the deltoid, False if not."""
return (x**2 + y**2)**2 + 18*(x**2 + y**2) - 27 < 8*(x**3 - 3*x*y**2)
def random_point():
"""Return a random point in the rectangle R=[-2, 3]x[-3,3]"""
x = -2 + 5*random.random_sample()
y = -3 + 6*random.random_sample()
return (x,y)
def deltoid_area(n):
"""Use Monte Carlo Integration with n points to estimate the area of the deltoid."""
total = 0
for i in range(n):
x,y = random_point()
if deltoid(x,y):
total += 1
return (total / n) * 30
deltoid_area(100000)
Remark: The actual area is $2 \pi$:
2 * np.pi
13. Classes¶
Write a class QField such that objects of the class represent numbers of the form $a+b\sqrt{2}$ where $a$ and $b$ are integers. For example, Qfield(3, 5) should return a QField object representing $3 + 5 \sqrt{2}$.
If x and y are two QField objects representing $a_x+b_x \sqrt{2}$ and $a_y+b_y \sqrt{2}$ respectively, the following should be possible:
x.float()should return a floating point approximation for $a_x+b_x \sqrt{2}$.x.plus(y)should return theQFieldobject representing $a_x+a_y+(b_x+b_y) \sqrt{2}$.x.times(y)should similarly return theQFieldobject representing the product (which also can be expressed as an integer plus an integer times the square root of two).
The formula for addition was given above but not multiplication. Observe $$(a_x + b_x \sqrt{2}) (a_y + b_y \sqrt{2}) = a_x a_y + 2 b_x by + (a_x b_y + a_y b_x) \sqrt{2}.$$
class QField:
def __init__(self, a, b):
self._a = a
self._b = b
def float(self):
return self._a + self._b * 2**0.5
def plus(self, y):
return QField(self._a + y._a, self._b + y._b)
def times(self, y):
return QField(self._a*y._a + 2*self._b*y._b, self._a*y._b + self._b*y._a)
# Testing. Construct the square root of 2:
sqrt2 = QField(0, 1)
print(sqrt2.float())
sqrt2.float() ** 2
# Test addition
x = QField(5, -2)
y = QField(-4, 3)
z = x.plus(y)
xf = x.float()
yf = y.float()
zf = z.float()
error = xf + yf - zf
error
# Test multiplication
x = QField(5, -2)
y = QField(-4, 3)
z = x.times(y)
xf = x.float()
yf = y.float()
zf = z.float()
error = xf*yf - zf
error
14. Shifting polynomials¶
Recall that a polynomial
$$p(x) = \sum_{k=0}^d c_k x^k$$
can be represented by a dictionary with one entry for every non-zero coefficient. The dictionary representing $p$ would then map a degree $k$ with a non-zero $c_k$ to the floating point number $c_k$. For example {9:3.0, 7:2.0} would represent $3 x^9 + 2x^7$.
Write a function shift(p, t) which takes as input a dictionary p representing a polynomial $p(x)$ and a floating point real number $t$. The function should return the dictionary representing the polynomial
$$q(x) = p(x+t).$$
Help: To do this, you need to know the coefficients of $(x+t)^k$. We have
$$(x+t)^k = \sum_{i=0}^k \binom{k}{i} t^{k-i} x^i.$$
The binomial coefficient $\binom{k}{i}$ can be computed using the binom function from the scipy.special package. For example, $\binom{4}{2}=6$ and this can be seen using the following python code:
from scipy.special import binom
binom(4, 2)First consider the math. Suppose $p(x) = \sum_{k=0}^d c_k x^k$. Then $$q(x) = p(x+t) = \sum_{k=0}^d c_k (x+t)^k.$$ Using the formula for $(x+t)^k$ we see that $$q(x) = \sum_{k=0}^d c_k \sum_{i=0}^k \binom{k}{i} t^{k-i} x^i = \sum_{k=0}^d \sum_{i=0}^k c_k \binom{k}{i} t^{k-i} x^i.$$
from scipy.special import binom
def shift(p, t):
q = {} # intitially an empty dictionary
for k in p:
ck = p[k]
for i in range(k+1):
coef = ck * binom(k, i) * t**(k-i)
# The following if statement adds coef*x^i to q.
if i in q:
q[i] += coef
else:
q[i] = coef
# Now if we have any zero coefficients, we remove them
# We can't change the keys of a dictionary we are interating through,
# so we convert the keys to a list first.
l = list(q.keys())
for k in l:
if q[k] == 0.0:
del q[k]
return q
p = {1: 1.0} # The polynomial p(x) = x
shift(p, 7) # Should be x+7
p = {2: 1.0, 1: 2.0, 0: 1.0} # The polynomial p(x) = x^2 + 2x + 1
shift(p, -1) # Should be x^2