Homework 12

CHE 116: Numerical Methods and Statistics

Prof. Andrew White

Version 1.0 (4/9/2015)



In [1]:

    
%matplotlib inline
import random
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt, pi, erf
import seaborn
seaborn.set_context("talk")
seaborn.set_style("whitegrid")
import scipy.stats

1. Rearranging ODEs (4 Points)

Convert the following ODEs into standard form. If it is a second-order, make sure you specify both dimensions.

$\frac{dx}{dt} - t = 2$
$\frac{dx}{dt} = -rx$
$\frac{d^2x}{dt^2} + k x = \sin t$
$\frac{d^2x}{dt^2} + b \frac{dx}{dt} + k x = 0$

Answers

$\frac{dx}{dt} = t + 2$
$\frac{dx}{dt} = -rx$
$\frac{d^2x}{dt^2} + k x = \sin t$

$\frac{dx_1}{dt} = x_2$

$\frac{dx_2}{dt} + kx_1 = \sin t$

$\frac{dx_2}{dt} = \sin t - kx_1$
$\frac{dx_1}{dt} = x_2$

$\frac{dx_2}{dt} =- b x_2 -k x_1$

2. Solve This ODE (6 Points)

Solve the tank problem from lecture:

$$\frac{dV}{dt} - k_1 \sqrt{V} = k_2$$

and use $k_1 = 1$, $k_2 = 4$. Plot the solution from $t=0$ to $t=10$ with the initial condition that $V_0 = 10$



In [13]:

    
def deriv_2(V, t, k_1=1, k_2=4):
    return k_2 + k_1 * np.sqrt(V)

vinit = 10.0
t_points = np.linspace(0,25,100)
v_points = scipy.integrate.odeint(deriv_2, vinit, t_points)

plt.plot(t_points, v_points)

plt.show()

3. Solve This Other ODE (8 Points)

Solve a forced, dampened ODE governed by this equation:

$$\frac{d^2x}{dt^2} + \frac{1}{10}\frac{dx}{dt} + k x = 5e^{-c / t}$$

using the interact command to change $c$ and $k$

$$\frac{dx_1}{dt} = x_2$$$$\frac{dx_2}{dt} + \frac{x_2}{10} + kx_1 = 5e^{-c / t}$$$$\frac{dx_2}{dt} = e^{-c / t} - \frac{x_2}{10} - kx_1$$



In [11]:

    
def deriv_3(x, t, k=1, c=1):
    return np.array([x[1], 5 * np.exp(-c / t) - x[0] * k - 0.1 * x[1]])

xinit = np.array([0.0, 1.0])
t_points = np.linspace(0.001,10, 100)

def solve_and_plot(k, c):
    x_points = scipy.integrate.odeint(lambda x,t: deriv_3(x,t, k, c), xinit, t_points)
    plt.plot(t_points, x_points[:,0])

from IPython.html.widgets import interactive

interactive(solve_and_plot, k=(0.01, 5, 0.1), c=(0.001, 1, 0.01))

4. Linear Algebra (3 Points)

Given that:

$$\mathbf{A} = \left[\begin{array}{lcr} 7 & 0 & 1\\ 1 & 3 & 4\\ 4 & 5 & 2\\ \end{array}\right]$$$$\mathbf{b} = \left[\begin{array}{lcr} 4\\ 13\\ -6\\ \end{array}\right]$$

Answer the following problems:

What is $\mathbf{A}\mathbf{b}$?
What is the largest eigenvector of $\mathbf{A}$?
Solve $\mathbf{A}\mathbf{x} = \mathbf{b}$



In [4]:

    
a_list = [[7, 0, 1],\
          [1, 3, 4],\
          [4, 5, 2]]
np_a = np.array(a_list)
np_b = np.array([4,13,-6]).transpose()



In [5]:

    
print np_a.dot(np_b)



In [6]:

    
import numpy.linalg as linalg

e_values, e_vectors = linalg.eig(np_a)

#The eigenvectors are in order of eignevalue
print e_vectors[:, 0]









    



[ 0.43264513  0.57338348  0.69573672]



In [7]:

    
np_a_inv = linalg.inv(np_a)
np_x = np_a_inv.dot(np_b)

print np_x









    



[-0.25714286 -3.31428571  5.8       ]



In [8]:

    
#to check, make sure we get back b
print np_a.dot(np_x)









    



[  4.  13.  -6.]