Homework 9 Key

CHE 116: Numerical Methods and Statistics

2/21/2019



In [30]:

    
import numpy as np

1. Creating Matrices (12 Points)

Create the following matrices using the given constraints. 4 Points each.

In 3 lines of python (not including prints/imports), create a 6x12 matrix whose second column (where we count "first", "second", etc.) is the powers of 2 (e.g., 1, 2, 4, 8, 16). Its third is all 4's. All other elements are zero
In 5 lines of python, accomplish the following: create a random 10x10 matrix using np.random.normal centered at 10 with standard deviation of 5, replace all negative values with 0 and then modify it so its rows sum to 1 and 0, use np.round to round to one decimal place.
In at most 3 lines of python, create a 9x9 matrix where all elements are 0 except in the diagonal, which is 1's

1.1



In [2]:

    
import numpy as np

m = np.zeros((6, 12))
m[:,1] = 2**np.arange(6)
m[:,2] = 4
print(m)









    



[[ 0.  1.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  2.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  4.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  8.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. 16.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. 32.  4.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

1.2



In [23]:

    
data = np.random.normal(size=(10,10), scale=5, loc=10)


data[data < 0] = 0

for i in range(10):
    data[i,:] /= np.sum(data[i,:])
print(np.round(data, 1))









    



[[0.1 0.1 0.1 0.1 0.1 0.  0.1 0.1 0.1 0.1]
 [0.1 0.1 0.  0.1 0.  0.2 0.  0.2 0.2 0.1]
 [0.  0.1 0.2 0.1 0.1 0.2 0.1 0.1 0.1 0.1]
 [0.1 0.1 0.  0.1 0.2 0.1 0.1 0.  0.1 0.1]
 [0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.  0.1 0.1]
 [0.2 0.2 0.1 0.1 0.1 0.  0.  0.1 0.1 0. ]
 [0.1 0.1 0.1 0.1 0.  0.1 0.1 0.1 0.1 0.1]
 [0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
 [0.1 0.1 0.1 0.2 0.1 0.1 0.2 0.1 0.1 0. ]
 [0.1 0.1 0.1 0.2 0.1 0.1 0.  0.1 0.1 0.1]]

1.3



In [22]:

    
data = np.zeros((9,9))
for i in range(9):
    data[i,i] = 1
print(data)









    



[[1. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 1.]]

2. Solving Systems of Equations (8 Points)

Solve the following systems of equations. 4 Points each. Write out your answer in Markdown

1.

$$\begin{array}{ll} 4x - 2y + z &= 0\\ 2x - 4y + z &= 1 \\ 2x + y + 3z &= 3\\ \end{array}$$

2.

$$\begin{array}{ll} e^{x} + y + 4z &= 4\\ 4e^{x} - 2y - z &= -1 \\ 3e^{x} + y + z &= 2\\ \end{array}$$

2.1



In [24]:

    
import numpy.linalg as  lin

A = np.array([[4, -2, 1], [2, -4, 1], [2, 1, 3]])
Ainv=  lin.inv(A)
B = np.array([0, 1, 3])
x = Ainv.dot(B)
print(x)









    



[-0.38235294 -0.11764706  1.29411765]

$$ x = -0.38\, y = -0.12\, z = 1.29\, $$

2.2



In [27]:

    
A = np.array([[1, 1, 4], [4, -2, -1], [3, 1, 1]])
Ainv=  lin.inv(A)
B = np.array([4, -1, 2])
x = Ainv.dot(B)
print(np.log(x[0]), x[1], x[2])









    



-1.5198257537444133 0.53125 0.8125

$$ x = -1.52\, y = 0.53\, z = 0.81 $$

3. Eigenvalue/Eigenvector Problems (8 Points)

Calculate the eigenvalues and eigenvectors for the following matrices. Solve in Python and then write out the eigenvalues/eigenvectors in LaTeX. When writing decimals, only report two significant figures.

[3 Points] $$A = \left[\begin{array}{lcr} 4 & 2 & 2\\ 4 & -8 & 4\\ 8 & 6 & -10\\ \end{array}\right]$$

[3 Points] $$A = \left[\begin{array}{lcr} 1 & 5 & 1\\ 2 & -1 & 2\\ 1 & 2 & -3\\ \end{array}\right]$$

[2 Points] Why would you use eigh over eig?

3.1



In [28]:

    
A = np.array([[4, 2, 2], [4, -8, 4], [8, 6, -10]])
print(lin.eig(A))









    



(array([  6.26181582,  -6.16107966, -14.10073615]), array([[ 0.77639311,  0.26691702, -0.03470413],
       [ 0.36237847, -0.74614174, -0.53196789],
       [ 0.51565064, -0.60994083,  0.84605307]]))

$$\lambda_1 = 6.3$$$$ v_1 = \left[0.78, 0.36, 0.52\right] $$$$\lambda_2 = -6.2$$$$ v_2 = \left[0.27, -0.75, -0.61\right] $$$$\lambda_3 = -14.1$$$$ v_3 = \left[-0.034, -0.53, 0.85\right] $$

3.2



In [52]:

    
A = np.array([[1, 5, 1], [2, -1, 2], [1, 2, -3]])
print(lin.eig(A))









    



(array([ 3.92715706, -2.57344229, -4.35371476]), array([[-0.85603357, -0.78017368,  0.4909431 ],
       [-0.45042445,  0.47652482, -0.64319425],
       [-0.25362244,  0.40528153,  0.58760193]]))

$$\lambda_1 = 3.9$$$$ v_1 = \left[-0.85, -0.78, 0.49\right] $$$$\lambda_2 = -2.6$$$$ v_2 = \left[-0.45, 0.48, -0.64\right] $$$$\lambda_3 = -4.4$$$$ v_3 = \left[-0.25, 0.40, 0.59\right] $$

3.3

If you had physical insight about the problem, like your eigenvalues provided solutions to an ODE, that disqualified complex eigenvalues you can use eigh. Otherwise, use eig.

4. Slicing Practice (6 Points)

Using numpy, create a sum or difference of array slices that yields the requested quantity. Consider the following example:

To create this sequence: $$x_0 , x_2, x_4, \ldots $$

Use this slice

x[::2]

Use this particular array for this: x = np.arange(15), but use len(x) when you need to refer to the length of the array. 2 Points each.

$x_1 \cdot x_0, x_2 \cdot x_1, \ldots x_N\cdot x_{N - 1}$
$ x_0 + x_N, x_1 + x_{N - 1}, \ldots, x_0 + x_N$
$x_N - x_{N - 1}, x_{N - 1} - x_{N - 2},\ldots ,x_1 - x_0$

4.1



In [20]:

    
x = np.arange(15)
x[1:] * x[:-1]









    Out[20]:





array([  0,   2,   6,  12,  20,  30,  42,  56,  72,  90, 110, 132, 156,
       182])

4.2



In [21]:

    
x + x[::-1]









    Out[21]:





array([14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14])

4.3



In [22]:

    
x[:0:-1] - x[-2::-1]









    Out[22]:





array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

5. Numerical Differentation/Integration Methods (5 Points)

Given the following problems, what is the correct method to use? Do not solve the problems, just state the best method. 1 Point each

Compute $\int_0^1 e^{-x^2}\,dx$
You are given $f(x)$ evaluated at the following x values: $[0,0.1, 0.5, 0.9, 1.2]$. Compute the derivative at $f(0.5)$.
Using the example from 5.2, compute $\int_{0.1}^{1.2} f(x)\,dx$
$g(x) = x^3$. What is the derivative $g$ at $g'(x = 0.5)$
Compute $\int_0^1 e^{-x}\,dx$

5.1

Gaussian quadrature (quad)

5.2

Central difference

5.3

Trapezoidal

5.4

Analytic

5.5

Analytic or quad

6. Numerical Differentation/Integration Methods (22 Points)

[4 Points] Compute $\int_{-\infty}^{\infty} x^2 e^{-x^2}\,dx$. Use np.inf to refer to infinity and use a lambda function. Make sure it's clear what is the value of the integral in your print.

[4 Points] Compute the numerical derivative of the following data:

x = [0, 1, 2, 3, 4, 6, 7, 9]
fx = [0.0, 0.84, 0.91, 0.14, -0.76, -0.28, 0.66, 0.41]

[2 Points] Compute the numerical derivative of the following data:

x2 = [0.0, 0.42, 0.83, 1.25, 1.67, 2.08, 2.5, 2.92, 3.33, 3.75, 4.17, 4.58, 5.0, 5.42, 5.83, 6.25, 6.67, 7.08, 7.5, 7.92, 8.33, 8.75, 9.17, 9.58, 10.0]
fx2 = [0.0, 0.4, 0.74, 0.95, 1.0, 0.87, 0.6, 0.22, -0.19, -0.57, -0.85, -0.99, -0.96, -0.76, -0.43, -0.03, 0.37, 0.72, 0.94, 1.0, 0.89, 0.62, 0.26, -0.16, -0.54]

[6 Points] Plot your data from 6.2 and 6.3 against $\cos(x)$, which is the correct derivative. Does the numerical derivative work even with the non-uniform coarse data in part 2?
[6 Points] Integrate the data from 6.2 and compare against the true integral $\int_0^{9} \sin(x)\,dx$. How accurate is it?

6.1



In [57]:

    
from scipy.integrate import quad
quad(lambda x: np.exp(-x**2) * x**2, -np.inf, np.inf)[0]









    Out[57]:





0.8862269254527599

6.2



In [58]:

    
x = np.array([0, 1, 2, 3, 4, 6, 7, 9])
fx = np.array([0.0, 0.84, 0.91, 0.14, -0.76, -0.28, 0.66, 0.41])
diff = (fx[1:] - fx[:-1]) / (x[1:] - x[:-1])
cdiff = (diff[1:] + diff[:-1]) / 2
print(cdiff)









    



[ 0.455  -0.35   -0.835  -0.33    0.59    0.4075]

6.3



In [59]:

    
x2 = np.array([0.0, 0.42, 0.83, 1.25, 1.67, 2.08, 2.5, 2.92, 3.33, 3.75, 4.17, 4.58, 5.0, 5.42, 5.83, 6.25, 6.67, 7.08, 7.5, 7.92, 8.33, 8.75, 9.17, 9.58, 10.0])
fx2 = np.array([0.0, 0.4, 0.74, 0.95, 1.0, 0.87, 0.6, 0.22, -0.19, -0.57, -0.85, -0.99, -0.96, -0.76, -0.43, -0.03, 0.37, 0.72, 0.94, 1.0, 0.89, 0.62, 0.26, -0.16, -0.54])
diff2 = (fx2[1:] - fx2[:-1]) / (x2[1:] - x2[:-1])
cdiff2 = (diff2[1:] + diff2[:-1]) / 2
print(cdiff2)









    



[ 0.89082462  0.66463415  0.30952381 -0.09901278 -0.47996516 -0.77380952
 -0.95238095 -0.95238095 -0.78571429 -0.50406504 -0.13501742  0.27380952
  0.64053426  0.8786295   0.95238095  0.90301974  0.68873403  0.33333333
 -0.06271777 -0.45557491 -0.75       -0.94076655 -0.96457607]

6.4



In [62]:

    
import matplotlib.pyplot as plt
%matplotlib inline
np.concatenate((fx[:1],cdiff,fx[-1:]))









    Out[62]:





array([ 0.    ,  0.455 , -0.35  , -0.835 , -0.33  ,  0.59  ,  0.4075,
        0.41  ])



In [65]:

    
plt.plot(x[1:-1], cdiff, label='6.2 data')
plt.plot(x2[1:-1], cdiff2, label='6.3 data')
plt.plot(x2, np.cos(x2), label='True derivative')
plt.legend()
plt.show()

The numerical derivative is quite accurate, even though it only exists at 4 points. It may not be accurate though at the end points, which have no central difference derivative

6.5



In [78]:

    
print(np.sum((fx[1:] + fx[:-1]) / 2 * (x[1:] - x[:-1])), quad(np.sin, 0, 9)[0])









    



1.7299999999999998 1.9111302618846777

The answer is off by only 0.2, which is not bad since there are only 6 data points.