In order to successfully complete this assignment, you must follow all the instructions in this notebook and upload your edited ipynb file to D2L with your answers on or before 11:59pm on Friday November 6th.
BIG HINT: Read the entire homework before starting.
Homework 4: Working with Linear Systems¶
In this homework we are going to explore many different problem examples. These examples are in the style of your bi-weekly quizzes and are provided as a way to practice.
Note, this is also a similar format to your final exam in this class.
Problem 1 - Traffic flow¶
Consider traffic going though W. Circle Dr. on the Michigan State University Campus. Note that the circle only allows "one way" traffic except the short section connecting Beal St. and Kalamazoo St. (which is two way). Also note that the roads connecting to the circle (Abbot Rd, Beal St., Kalamazoo St., Auditorium Rd and E. Circle Dr.) are all two way streets with traffic flowing both in and out of the Circle.
In the above figure, each intersection has been given a single letter variable corresponding where each road intersects with W. Circle Dr. These are as follows:
- U - Union (at Abbot Rd.)
- B - Beal St.
- K - Kalamazoo St.
- A - Auditorium Rd.
- E - E. Circle Dr.
On a busy game day, flow into and out of the circle at each of the intersections can be given the following values (units are vehicles per minute):
# Here are some libraries you may need to use
%matplotlib inline
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
sym.init_printing(use_unicode=True)
Uin = 20
Uout = 25
Bout = 33
Bin = 20
Kout = 18
Kin = 33
Aout = 55
Ain = 32
Eout = 7
Ein = 33
The unknowns in this system are the traffic flow values around the circle ($x_1, x_{2a}, x_{2b}, x_3, x_4, x_5$). We can assume that the total traffic flowing into an intersection is equal to the total traffic flowing out of an intersection. For example, The flow going into and out of the Union intersection (U) and the Beal Interscation (B) can be dented by the following equations:
$$ U_{in} + x_5 = U_{out} + x_1$$$$B_{in} + x_1 + x_{2a} = B_{out} + x_{1b}$$**Question 1.a:** Write out the remaining 3 linear equations (using latex syntax) that describe the traffic flow though each of the remaining intersections.
YOUR ANSWER HERE
**Question 1.b:** Transform all five equations into a linear system of the form $Ax=b$ and represent that system using an augmented matrix using the values in the python variables defined above (Uin, Uout, Bout, etc)
# Put your augmented matrix here.
from answercheck import checkanswer
checkanswer.detailedwarnings = False
checkanswer.eq_matrix(M,"64ed5d83acc4424c0bd5f3e9078a51b1");
**Question 1.c:** Calculate the reduced row echelon form of your equation and show the result below.
# YOUR CODE HERE
raise NotImplementedError()
**Question 2.d:** Rewrite the above answer to express each leading variable (first variable in each row) in terms of the remaining variables. For example, the first row's leading variable is $x_1$ which can be written out as $x_1 = x_5 - 5$. Do this same exercise for the remaining rows and write the answer using Latex syntax in markdown.
YOUR ANSWER HERE
**Question 1.e:** Assuming that $x_{2b}$ is 9, what is the minimum value of flow across $x_5$. (hint there can not be negative flow on a one way street). Given these values for $x_{2b}$ and $x_5$, what are the values for all the unknowns?
# Put your answer to the above question here.
from answercheck import checkanswer
checkanswer.vector([x1, x2a, x2b, x3, x4, x5],"a66bb881d0e6af98ce96b21ac090559d");
Problem 2- Wheatstone bridge¶
A common applications of linear algebra to electronics is to analyze electronic circuits. This problem is very similar to the traffic flow problem above. The goal is to calculate the current flowing in each branch of the circuit or to calculate the voltage at each node of the circuit. The following circuit represents a Wheatstone bridge which can be used to measure an unknown resistance ($R_G$).
Knowing the branch currents, the nodal voltages can easily be calculated, and knowing the nodal voltages, the branch currents can easily be calculated. Loop analysis finds the currents directly and nodal analysis find the voltages directly. Which method is simpler depends on the given circuit.
Lets assume the following values for the above circuit:
V1 = 9 #volts
R1 = 5 #Ohms
R2 = 2 #Ohms
Rg = 50 #Ohms
R3 = 10 #Ohms
R4 = 4 #Ohms
Like traffic flow, Kirchoff’s current law states that the sum of the currents flowing into a node equals the sum of the currents flowing out of the node. Here is an example:
Image from: commons.bcit.ca
Using the the Wheatstone bridge figure above we can apply Kirchoff's current law to three of the four nodes and get the following:
$$I_0 = I_1 + I_2$$$$I_1 = I_3 + I_g$$$$I_2 + I_g = I_4$$✅ **Question 2.a**: What is the equation you get if you apply Kirchoff's current law to the remaining (bottom) node?
YOUR ANSWER HERE
Ohm’s Law The voltage difference between the two ends of a resistor causes a current to flow through the resistor. For many substances the voltage and current are proportional. This is expressed by the formula:
$$V = I R$$This equation is called Ohm’s law and any device that obeys it is called a resistor.
$V$ is the difference in voltage between the two ends of the resistor measured in volts, $I$ is the current through the resistor measured in amperes, and the proportionality constant $R$ is the resistance of the resistor measured in ohms. Given any two of these quantities, Ohm’s law can be used to find the third.
Kirchoff’s Voltage Law states that, Around any closed path (loop) in an electric circuit, the sum of the voltage rises through the voltage sources equals the sum of the voltage drops through the resistors.
A closed path is a path through a circuit that ends where it starts.
For example: $I_0 \rightarrow I_1 \rightarrow I_3$ is a closed path. We can write this path using Kirchoff's Voltage law and Ohm's Law to come up with the following equation: $$V_1 = I_1R_1 + I_3R_3$$
Similarly, the equation for path through $I_1 \rightarrow I_g \rightarrow I_2$ is :
$$R_1I_1 + R_gI_g - R_2I_2 = 0$$✅ **Question 2.b**: What is the eaution for the path though $I_0 \rightarrow I_2 \rightarrow I_4$?
YOUR ANSWER HERE
If we ignore the answers to questions a and b and combine the three equations provided by Kirchoff's current law and the two equations provided by Kirchoff's Voltage Law, we get the following set of linear equations (represented as $Ax=b$ with $x = [I_0, I_1, I_2, I_3, I_4, I_g]^T$):
A = [[1, -1, -1, 0, 0, 0],
[0, 1, 0, -1, 0, -1],
[0, 0, 1, 0, -1, 1],
[0, R1, 0, R3, 0, 0],
[0, R1, -R2, 0, 0, Rg]]
sym.Matrix(A)
b = [[0], [0], [0], [V1], [0]]
sym.Matrix(b)
Unfortunately there are six unknowns and only five equations. We do not have enough information to find a unique solution to the problem.
✅ **Question 2.c**: Use your answers from question b to write a new system of linear equations such that there are now six (6) equations and (6) unknowns. Write out the the new values for $A$ and $b$ matrix in python.
# YOUR CODE HERE
raise NotImplementedError()
✅ **Question 2.d**: Using numpy, verify that the new matrix $A$ is invertible. Explain your answer in the space provided.
# YOUR CODE HERE
raise NotImplementedError()
Explain your answer to part d here.
✅ **Question 1.e**: Solve the system of equations, what are the current values for $I_0,I_1,I_2,I_3,I_4,I_g$?
#Put your answer here
from answercheck import checkanswer
checkanswer.vector(x,"2831d5d99336ac40955dad094df2f8b1");
Problem 3 - Large Scale Simulations¶
My types of large-scale simulations break a virtual world into small parts. This process is sometimes called meshing and involves creating a numerical representation of a real world structure. This type of simulation is used to understand weather, the formation of the universe, and how to build better airplanes and cars.
Each individual mesh element is simple compared to the entire simulation. Often there are only a few variables at each mesh element that need to be solved. However, the variables are dependent on variables in neighboring mesh elements.
One way to formulate the problem is to construct a very large matrix to solve all of the variables at once. The size of the matrix depends on the number of elements in a model but it is not uncommon to get matrices that are a million $\times$ million elements.
Solving these incredibly large-scale systems of equations can be done but requires a ton of computational time.
In this question we are going to discuss Big-O notation for comparing algorithm complexity. Watch the following video as an introduction to Big-O. If you do not understand, see if you can find other descriptions or tutorials on the web.
Cramer's Rule¶
✅ **Question 3.a:** In class we talked about using Cramer's method to calculate the solution to a system of equations. Explain why you would never use Cramer's method to solve such a large system. What specific calculation does Cramer's method use and why does that calculation cause a problem.
YOUR ANSWER HERE
Gauss-Jordon¶
Another option is to use the Gauss-Jordon method to solve the system of equations. The majority of computation in Gauss-Jordan occurs during the matrix multiply. In fact the Matrix Multiply dominates the calculation such that the Big-O complexity of Gauss-Jordan can be estimated as the complexity of matrix multiply. The following is a simple matrix multiply function defined in Python.
import numpy as np
A = [[3,4,5],[1,2,3]]
B = [[3, 3], [2,3], [1,3]]
C = [[0,0], [0,0]]
np.matrix(A)*B
#Simple Matrix Multiply
for i in range(len(A)):
for k in range(len(B[0])):
C[i][k] = 0;
for j in range(len(B)):
C[i][k] += A[i][j] * B[j][k];
C
✅ **Question 3.b:** Given the above, what is the Big-O complexity of this naive matrix multiply (and similarly Gauss-Jordon using this version of matrix multiply)?
YOUR ANSWER HERE
Although easy to code, the above matrix multiply is not very good. According to wikipedia (https://en.wikipedia.org/wiki/Matrix_multiplication_algorithm) The current $O(n^k)$ algorithm with the lowest known exponent $k$ is a generalization of the Coppersmith–Winograd algorithm that has an asymptotic complexity of $O(n^{2.3728639})$. This means there is a version of Gauss-Jordan that also has a Big-O complexity of $O(n^{2.3728639})$.
Jacobi¶
Remember from class that the Jacobi method is iterative. Each iteration of the algorithm takes $O(n)$ time. If the algorithm is guaranteed to converge after $m$ iterations then the Big-O complexity is $O(mn)$.
For some (not all problems) we can assume that the number of iterations is much less than the size of the matrix (i.e. $m<<n$), in these situations the Big-O complexity of Jacobi is $O(n)$. In other cases we say $m$ is approximately equal to $n$ so the Big-O complexity can be estimated as $O(n^2)$.
✅ **Question 3.c**: Given some case in which the number of iterations $m$ is approximately equal to $log(n)$. What is the Big-O complexity the Jacobi method in this case?
YOUR ANSWER HERE
Which is best?¶
Big-O notation compares algorithms when $n$ gets really big. Consider the algorithms depicted in the following made-up graph. In this example the linear Jacobi method is complex and adds a constant factor of 20 to it's calculations and the Gauss-Jordan method adds a constant factor 2. For big values of $n$, these constant factors are meaningless, but for small values of $n$ they can be important.
%matplotlib inline
import matplotlib.pylab as plt
import scipy as sp
from scipy.special import factorial
import numpy as np
mx_x = 10
n = np.linspace(1,mx_x)
plt.plot(n,n*20, label='Jacobi $O(n)$')
plt.plot(n,n*n, label='Jacobi $O(n^2)$')
plt.plot(n,(n*2)**2.3728639, label='Gauss-Jordan $O(n^{2.3728639})$')
plt.plot(n,factorial(n), label='Cramer\'s Rule $O(n!)$')
plt.axis([2,mx_x, 0,3000])
plt.xlabel('number of inputs (n)')
plt.ylabel('approx. calculation time')
plt.legend();
✅ **Question 3.d**: Given the above data, at what approximate value of $n$ is it faster to use Gauss-Jordan compared to Cramer's Rule?
YOUR ANSWER HERE
✅ **Question 3.e**: Explain why numerical libraries such as numpy.linalg solve equations of the form $Ax=b$ using Gauss-Jordan, even though Gauss-Jordan is slower than either of the Jacobi methods for most cases,
YOUR ANSWER HERE
Example motivated by (hint good resource if you are stuck): https://www.countbayesie.com/blog/2015/11/21/the-black-friday-puzzle-understanding-markov-chains
Consider the map of the book/gift store shown in the figure above. The store is shaped like a pentagon with each department on the outer edge. To the right of the map we see a Markov state transition model that describes the probabilities of customers moving between departments.
This probability model can be represented by the following matrix.
# Here are some libraries you may need to use
%matplotlib inline
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
sym.init_printing(use_unicode=True)
names = ['Books', 'Childrens', 'Toys', 'Puzzles', 'Music']
A = np.matrix([[ 0.5,0.1, 0.1, 0.05, 0.2],
[0.2, 0.3, 0.2, 0.15, 0.1],
[ 0.15, 0.2,0.3, 0.3, 0.1],
[ 0.1,0.3, 0.2 ,0.4 ,0.1],
[0.05,0.1 ,0.2,0.1,0.5]])
print(names)
sym.Matrix(A)
✅ **Question 4.a:** Assume the store opens and starts with exactly 50 people in each room and no one leaves or enters the store. How many people does the Markov model predict will be in each department in exactly two time steps?
#Put your answer to the above question here.
from answercheck import checkanswer
checkanswer.vector(b2,"b741dbee33e52d2fa3ea783aa8949d2c");
✅ **Question 4.b:** Use numpy to calculate the eigenvalues and eigenvectors for the matrix $A$. Note some of these eigenvectors may have imaginary values.
# YOUR CODE HERE
raise NotImplementedError()
✅ **Question 4.c:** The eigenvector associated with the largest eigenvalue (we can call it $v$) represents the long term steady state probably of the system. Rescale the vector $v$ into a new vector called $\bar{v}$ such that it's components represent probabilities that add to one. In other words, calculate a constant $c$ such as the following is true:
$$v = [v_1,v_2,v_3,v_4,v_5]$$$$\bar{v} = cv$$Where $c$ above fulfills the following equalities:
$$\begin{align} cv_1 + cv_2 + cv_3 + cv_4 + cv_5 &= \bar{v}_1 + \bar{v}_2 + \bar{v}_3 + \bar{v}_4 + \bar{v}_5\\ &= 1 \end{align} $$# YOUR CODE HERE
raise NotImplementedError()
✅ **Question 4.d:** Using $\bar{v}$ calculated above create a new $5 \times 5$ steady state matrix with $\bar{v}$ each of the five columns of the matrix.
$$ [ \bar{v}^T | \bar{v}^T | \bar{v}^T | \bar{v}^T | \bar{v}^T ] $$#Put your answer to the above question here.
from answercheck import checkanswer
checkanswer.matrix(M,"1052b1bc9525b87d6fac25520ae004c0");
✅ **Question 4.e:** Assume the store opens and starts with exactly 50 people in each room and no one leaves or enters the store. Using the Markov model, what is the long term steady state of the system. i.e. how many people are likely to be in each room as time approaches infinity?
#Put your answer to the above question here.
from answercheck import checkanswer
checkanswer.vector(b,"fe16c592cc06d695ee4db42562d7ee2a");
Congratulations, we're done!¶
Turn in your assignment using D2L no later than 11:59pm on the day of class. See links at the end of this document for access to the class timeline for your section.
Written by Dirk Colbry and Ming Yan, Michigan State University
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.