Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

compmech-project04/notebooks/01_Linear-Algebra.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
1028 lines (1028 sloc) 37.5 KB
Raw Blame
Open in GitHub Desktop
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Linear Algebra Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Linear algebra is a representation of linear equations. We can separate the known constants from unknown variables to represent our system of equations. Take the following example, there are three unknowns and three equations\n",
"\n",
"1. $5x_{1}+3x_{2} =1$\n",
"\n",
"2. $x_{1}+2x_{2}+3x_{3} =2$\n",
"\n",
"3. $x_{1}+x_{2}+x_{3} =3$\n",
"\n",
"If we can represent our equations in a matrix-vector format, then we can use standard linear algebra routines to solve for the unknown variables, $x_1,~x_2,~and~x_3.$\n",
"\n",
"Consider the matrix form of equations 1-3 above:\n",
"\n",
"$\\left[ \\begin{array}{ccc}\n",
"5 & 3 & 0 \\\\\n",
"1 & 2 & 3 \\\\\n",
"1 & 1 & 1 \\end{array} \\right]\n",
"\\left[\\begin{array}{c} \n",
"x_{1} \\\\ \n",
"x_{2} \\\\\n",
"x_{3}\\end{array}\\right]=\\left[\\begin{array}{c} \n",
"1 \\\\\n",
"2 \\\\\n",
"3\\end{array}\\right]$\n",
"\n",
"$\\mathbf{Ax}=\\mathbf{b}$\n",
"\n",
"The __matrix__, $\\mathbf{A}$ contains all of the constants that are multiplied by our unknown variables $x_1,~x_2,~and~x_3.$ The __vector__, $\\mathbf{y}$ contains all of the known constants that are not multiplied by our unknown variables $x_1,~x_2,~and~x_3.$ Finally, the __vector__, $\\mathbf{x}=[x_1,~x_2,~x_3]$ contains our unknown values. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"Solve for $x_1,~x_2,~and~x_3.$ How did you approach the problem? Did you use Linear algebra?\n",
"\n",
"Highlight for answers: <font color=\"white\">x1 = 20, x2 = -33, x3 = 16</font>"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ -7.0 3.0 0.0] * [1.2] = [-20.0]\n",
"[ 7.0 -19.0 12.0] * [1.1] = [ 0.0]\n",
"[ 0.0 4.0 -12.0] * [1.0] = [ -8.0]\n"
]
}
],
"source": [
"A=np.array([[-7,3,0],[7,-19,12],[0,4,-12]])\n",
"b=np.array([-20,0,-8])\n",
"\n",
"np.linalg.solve(A,b)\n",
"for i in range(0,3):\n",
" print('[{:5.1f} {:5.1f} {:5.1f}] {} [{:3.1f}] {} [{:5.1f}]'.format(*A[i],'*',x[i],'=',b[i]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use linear algebra _operations_ to standardize the way we approach a set of equations. If we knew another matrix $\\mathbf{A}^{-1}$ that we could multiply both sides of our equation and get a new equation\n",
"\n",
"$\\mathbf{A^{-1} A x}= \\mathbf{A_n b}$\n",
"\n",
"where $\\mathbf{A^{-1} A x} = \\mathbf{x}$\n",
"\n",
"then $\\mathbf{x} = \\mathbf{A^{-1} b}$\n",
"\n",
"if this were a single equation with a single unknown, then we just use the inverse of A as such\n",
"\n",
"$12 x = 6$\n",
"\n",
"$\\frac{1}{12} 12 x = \\frac{1}{12} 6$\n",
"\n",
"$x=2$\n",
"\n",
"In this notebook, we will look at how to frame our problems as linear algebra problems and work on how to get $\\mathbf{A}^{-1}$ or $\\mathbf{x}$ in the next two notebooks. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example with Mixing Tanks\n",
"\n",
"![Mixing tank volume flow rates](../images/mixing_tanks.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the diagram above we have three tanks of water that are mixing two concentrations of salt water with $5~mg/m^3$ entering tank 1 and $1~mg/m^3$ entering tank three. The outlet is located on the middle tank 2, but the concentration is unknown. \n",
"\n",
"The volume flow rates of water are shown on each arrow to denote flow in/out of each tank. We want to know what the final concentration at the outlet is $c_2$, but in order to know that we need the concentrations of $c_1~and~c_3$. The mass flow of the salt is the concentration $\\times$ volume flow rate. The total mass flow in - mass flow out in each container is 0. We have three mass balance equations\n",
"\n",
"1. $(-7~m^3/s)~c_1 +(3~m^3/s)~c_2 +(4~m^3/s)~5~mg/m^3 = 0$\n",
"\n",
"2. $(7~m^3/s)~c_1 -(3~m^3/s)~c_2 -(4~m^3/s)~c_2 -(12~m^3/s)~c_2 +(12~m^3/s)~c_3 = 0$\n",
"\n",
"3. $(4~m^3/s)~c_2 -(12~m^3/s)~c_3 + (8~m^3/s)~1~mg/m^3 = 0$\n",
"\n",
"or rearranging our mass-balance equations we have\n",
"\n",
"1. $-7c_1+3c_2=-20$\n",
"\n",
"2. $7c_1-19c_2+12c_3=0$\n",
"\n",
"3. $4c_2-12c_3=-8$\n",
"\n",
"We can put this into the same form that we used above with the matrix $\\mathbf{A}$ and vectors $\\mathbf{x}$ and $\\mathbf{b}$\n",
"\n",
"$\\left[ \\begin{array}{ccc}\n",
"-7 & 3 & 0 \\\\\n",
"7 & -19 & 12 \\\\\n",
"0 & 4 & -12 \\end{array} \\right]\n",
"\\left[\\begin{array}{c} \n",
"c_{1} \\\\ \n",
"c_{2} \\\\\n",
"c_{3}\\end{array}\\right]=\\left[\\begin{array}{c} \n",
"-20 \\\\\n",
"0 \\\\\n",
"-8\\end{array}\\right]$\n",
"\n",
"$\\mathbf{Ax}=\\mathbf{b}$\n",
"\n",
"Now, let's use some numpy linear algebra to solve for $c_2$. First, define $\\mathbf{A}$ and $\\mathbf{b}$ our known constants. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"matrix A:\t vector b:\n",
"[-7 3 0] \t -5\n",
"[ 7 -19 12] \t 0\n",
"[ 0 4 -12] \t -8\n"
]
}
],
"source": [
"A=np.array([[-7,3,0],[7,-19,12],[0,4,-12]])\n",
"b=np.array([-5,0,-8])\n",
"print('matrix A:\\t vector b:')\n",
"for i in range(0,3):\n",
" print(A[i],'\\t',b[i])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can solve for $\\mathbf{x}$ with the function `np.linalg.solve`. This is an advanced linear algebra solver that incorporates some ideas we will explore in [Module 02](./02_Gauss_elimination.ipynb). For now, we just want to understand the inputs and outputs\n",
"\n",
"```python\n",
"x = np.linalg.solve(A,b)\n",
"```\n",
"\n",
"In the next cell, we run this line of code. The inputs are the matrix $\\mathbf{A}$ and vector $\\mathbf{b}$, as defined above as `A` and `b`. The output is our unknown vector $\\mathbf{x}=[c_1,~c_2,~c_3]$. If we plug in the values of `x` into our mass balance equations we will see that mass is conserved. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"c1 = 1.18 mg/m^3,\n",
"c2 = 1.08 mg/m^3,\n",
"c3 = 1.03 mg/mm^3\n"
]
}
],
"source": [
"x = np.linalg.solve(A,b)\n",
"print('c1 = {:.2f} mg/m^3,\\nc2 = {:.2f} mg/m^3,\\nc3 = {:.2f} mg/mm^3'.format(*x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"Show that $\\mathbf{Ax} = \\mathbf{b}$ in the previous mixing container example. Plug the values of `x` into the three equations and show that mass is conserved. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Vectors "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"We use vectors to represent unknown variables or known outputs. In numpy, a vector only has one dimension. \n",
"\n",
"```python\n",
"y = np.array([1,2,3])\n",
"```\n",
"\n",
"If you ask for the `shape` of `y`, you get an output of `(3,)`, which means it is a one dimensional vector. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(3,)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"y=np.array([1,2,3])\n",
"\n",
"y.shape"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### What's a vector?\n",
"\n",
"Vectors are everywhere: physics, engineering, mathematics, computer science, video games, and more. Each field's interpretation of what a vector *is* may be different, but vectors live a similar life in every space.\n",
"\n",
"The first episode in the wonderful video series, [_\"Essence of Linear Algebra\"_](http://3b1b.co/eola) tells you of three different ideas about vectors [1]:\n",
"\n",
"1. For physicists, a vector is an \"arrow\" of a given length (magnitude) and direction. It can represent directional quantities like velocity, force, acceleration.\n",
"2. For computer scientists, a vector is an ordered list of numbers. It can represent a set of variables or features stored in order.\n",
"3. For mathematicians, vectors are generic objects that behave a certain way when they are added or scaled: $\\mathbf{u}+\\mathbf{v}$, $\\alpha\\mathbf{v}$.\n",
"\n",
"<img src=\"../images/whatsavector.png\" style=\"width: 500px;\"/> "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### How you think of a vector depends on what you're doing\n",
"\n",
"In physics, vectors are almost always two- or three-dimensional (although in some fancy branches of physics they do go to higher dimensions). Vectors help physicists describe things like motion and electro-magnetic fields on a plane or in physical 3D space.\n",
"\n",
"But, as we saw in our first example of Linear algebra for a set of equations, the vector could be a set of known or unknown values. This is closer to how vectors are treated in computer science and data science, vectors are often multi-dimensional, that is, they have many components. They contain a set of ordered variables in a data model, like for example: the age, weight, daily hours of sleep, weekly hours of exercise, and blood pressure of an individual (five dimensions)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adding and subtracting scaled vectors\n",
"\n",
"In our first linear algebra problem, we had the vector, $\\mathbf{x}=[c_1,~c_2,~c_3]$ and solved for $\\mathbf{x}=[3.86,~2.33,~1.44]~mg/m^3$. We separated the flow rates out of the equation, but we could also have pulled out the flow rates in three vectors:\n",
"\n",
"$y=Ax=\\left[\\begin{array}{c} \n",
"-20 \\\\ \n",
"0 \\\\\n",
"-8\\end{array}\\right]\n",
"=\n",
"\\left[\\begin{array}{c} \n",
"-7 \\\\ \n",
"7 \\\\\n",
"0\\end{array}\\right] c_{1}+\n",
"\\left[\\begin{array}{c} \n",
"3 \\\\ \n",
"-19 \\\\\n",
"4\\end{array}\\right] c_{2}+\n",
"\\left[\\begin{array}{c} \n",
"0 \\\\ \n",
"12 \\\\\n",
"-12\\end{array}\\right] c_{3}$\n",
"\n",
"or \n",
"\n",
"$\\left[\\begin{array}{c} \n",
"-20 \\\\ \n",
"0 \\\\\n",
"-8\\end{array}\\right]\n",
"=\n",
"\\left[\\begin{array}{c} \n",
"-7 \\\\ \n",
"7 \\\\\n",
"0\\end{array}\\right] 3.86+\n",
"\\left[\\begin{array}{c} \n",
"3 \\\\ \n",
"-19 \\\\\n",
"4\\end{array}\\right] 2.33+\n",
"\\left[\\begin{array}{c} \n",
"0 \\\\ \n",
"12 \\\\\n",
"-12\\end{array}\\right] 1.44 = \n",
"\\left[\\begin{array}{c} \n",
"-20 \\\\ \n",
"0 \\\\\n",
"-8\\end{array}\\right]$\n",
"\n",
"When you multiply a vector by a scalar (something with a single magnitude, a number) the result is that each component of the vector is multiplied by that scalar. That's why we can separate the flow rates and multiply them by the individual concentrations and add the results. \n",
"\n",
"$y = a_{1}x_{1} + a_{2}x_{2} +...+a_{N}x_{N}$\n",
"\n",
"- $a_{i}$ is a column vector \n",
"\n",
"- $x_{i}$ is a scalar taken from the $i^{th}$ element of x.\n",
"\n",
"Multiplying a vector by a scalar is the same as multiplying each component of the vector by the scalar. So if we multiply a vector $\\mathbf{a}$ by 2, then it is the same as multiplying each component, $a_i$ by 2. \n",
"\n",
"$2\\mathbf{a}=\\left[\\begin{array}{c} \n",
"2a_{1} \\\\ \n",
"2a_{2} \\\\\n",
"\\vdots \\\\\n",
"2a_{n}\\end{array}\\right]$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"external inputs\n",
"flow in at 1: 20.0 g/s\n",
"flow in at 2: 0.0 g/s\n",
"flow in at 3: 8.0 g/s\n"
]
}
],
"source": [
"a1 = np.array([-7,7,0])\n",
"a2 = np.array([3,-19,4])\n",
"a3 = np.array([0,12,-12])\n",
"\n",
"out = x[0]*a1+x[1]*a2+x[2]*a3\n",
"\n",
"print('external inputs')\n",
"print('flow in at 1: {:.1f} g/s\\nflow in at 2: {:.1f} g/s\\nflow in at 3: {:.1f} g/s'.format(*out*-1))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Representation of a problem with Matrices and Vectors"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"There are two main types of problems that are solved with linear algebra. Let's use the following variables to describe our problems. \n",
"\n",
"- $\\mathbf{y}:$ a set of known outputs, $y_{1},~y_{2},~...y_{N}$ \n",
"\n",
"- $\\mathbf{A}:$ a set of known constants from equations, $A=\\left[ \\begin{array}{cccc}\n",
"A_{11} & A_{12} &...& A_{1N} \\\\\n",
"A_{21} & A_{22} &...& A_{2N} \\\\\n",
"\\vdots & \\vdots &\\ddots& \\vdots \\\\\n",
"A_{M1} & A_{M2} &...& A_{MN}\\end{array} \\right]$\n",
"\n",
"- $\\mathbf{x}:$a set of unknown inputs, $x_{1},~x_{2},~...x_{N}$\n",
"\n",
"- $\\lambda:$ an unknown constant\n",
"\n",
"Using the variables defined above we describe the two main types of linear algebra problems:\n",
"\n",
"1. linear system solutions where $\\mathbf{Ax} = \\mathbf{b}$\n",
"\n",
"2. eigenvalue soslutions where $\\mathbf{Ax} = \\lambda \\mathbf{x}$\n",
"\n",
"We saw an example of the first type, in the mixing example. Eigenvalue problems come up in a number of engineering and science applications. We will cover the application of eigenvalues in the last module when we look at boundary value problems. We will restrict our initial applications of linear algebra to linear systems. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. Another linear system solution\n",
"\n",
"![Example of Applying Newton's second law to obtain a linear system](../images/mass-pulley.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the above diagram, there are 4 masses, connected by 4 cords. Our goal is to create a system of equations that are solveable with Linear Algebra. We start with Newton's second law to determine 4 equations of motion\n",
"\n",
"1. $m_1 a_1 = T_1 + \\mu m_1g\\cos\\theta - m_1g\\sin\\theta$\n",
"\n",
"2. $m_2 a_2 = T_2 -T_1 + \\mu m_2g\\cos\\theta - m_2g\\sin\\theta$\n",
"\n",
"3. $m_3 a_3 = T_3 -T_2 + \\mu m_3g\\cos\\theta - m_3g\\sin\\theta$\n",
"\n",
"4. $m_4 a_4 = T_3 - m_4g$\n",
"\n",
"This gives us four equations with 7 unknowns $(a_1,~a_2,~a_3,~a_4,~T_1,~T_2,~and~T_3).$ We also have 3 constraints that relate the motion of masses 1-4:\n",
"\n",
"1. $a_1 = a_2$\n",
"\n",
"2. $a_2 = a_3$\n",
"\n",
"3. $-a_3 = a_4$\n",
"\n",
"So we can limit our description of acceleration to just $a$, because the system only has one overall degree of freedom (as long as the connecting cords remain taut). Now we have four equations and four unknowns $(a,~T_1,~T_2,~and~T_3).$\n",
"\n",
"1. $-m_1 a - T_1 = \\mu m_1g\\cos\\theta - m_1g\\sin\\theta$\n",
"\n",
"2. $-m_2 a - T_2 +T_1 = \\mu m_2g\\cos\\theta - m_2g\\sin\\theta$\n",
"\n",
"3. $-m_3 a - T_3 +T_2 = \\mu m_3g\\cos\\theta - m_3g\\sin\\theta$\n",
"\n",
"4. $m_4 a - T_3 = - m_4g$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or in the matrix-vector form:\n",
"\n",
"$\\left[ \\begin{array}{cccc}\n",
"-m_1 & -1 & 0 & 0 \\\\\n",
"-m_2 & 1 & -1 & 0\\\\\n",
"-m_3 & 0 & 1 & -1\\\\\n",
"m_4 & 0 & 0 & -1\\end{array} \\right]\n",
"\\left[\\begin{array}{c} \n",
"a \\\\ \n",
"T_1 \\\\\n",
"T_2 \\\\\n",
"T_{3}\\end{array}\\right]=\\left[\\begin{array}{c} \n",
"\\mu m_1g\\cos\\theta - m_1g\\sin\\theta\\\\\n",
"\\mu m_2g\\cos\\theta - m_2g\\sin\\theta\\\\\n",
" \\mu m_3g\\cos\\theta - m_3g\\sin\\theta\\\\\n",
"- m_4g\\end{array}\\right]$\n",
"\n",
"$\\mathbf{Ax}=\\mathbf{b}$\n",
"\n",
"Let's use the following constants to solve for acceleration and tensions:\n",
"\n",
"* $\\mu = 0.2$\n",
"* $m_1 = m_2 = m_3 = 2~kg$\n",
"* $m_4 = 1~kg$\n",
"* $\\theta=\\dfrac{\\pi}{4}=45^o$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A2 = \t\t y2=\n",
"[-2 -1 0 0] \t -12.486091542192055\n",
"[-2 1 -1 0] \t -12.486091542192055\n",
"[-2 0 1 -1] \t -12.486091542192055\n",
"[ 1 0 0 -1] \t -9.81\n"
]
}
],
"source": [
"mu = 0.1\n",
"m1=m2=m3=2\n",
"m4=1\n",
"g =9.81\n",
"a = np.pi/4\n",
"\n",
"A2 = np.array([[-m1,-1,0,0],[-m2,1,-1,0],[-m3,0,1,-1],[m4,0,0,-1]])\n",
"y2 = np.array([mu*m1*g*np.cos(a)-m1*g*np.sin(a),\\\n",
" mu*m2*g*np.cos(a)-m2*g*np.sin(a),\\\n",
" mu*m3*g*np.cos(a)-m3*g*np.sin(a),\\\n",
" -m4*g])\n",
"\n",
"print('A2 = \\t\\t y2=')\n",
"for i in range(0,4):\n",
" print(A2[i],'\\t',y2[i])"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a=3.95 m/s/s\n",
"T1=4.6 N\n",
"T2=9.2 N\n",
"T3=13.8 N\n"
]
}
],
"source": [
"x2 = np.linalg.solve(A2,y2)\n",
"\n",
"print('a={:.2f} m/s/s\\nT1={:.1f} N\\nT2={:.1f} N\\nT3={:.1f} N'.format(*x2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Exercise\n",
"\n",
"1. Plug in the values that we solved into the original 4 equations. Show that our values for accelerations and tensions satisfy our initial equations. \n",
"\n",
"2. Create a new vector `y3` where the coefficient of friction is $\\mu=0.5$ what is the new acceleration? Do the tensions increase or decrease when tension increases?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pause and ponder\n",
"\n",
"In this example, the unknown vector, $\\mathbf{x}=[a,~T_1,~T_2,~T_3],$ is a combination of acceleration and forces. This definition of a vector is less intuitive than a physics-based magnitude and direction, but it is __extremely__ useful in solving engineering and physics problems. Here are a few __key ideas__ from these two linear system exercises:\n",
"\n",
"* In order to solve for $n$ unknowns, we need $n$ independent equations \n",
"* A vector is a collection of numbers that we either know or want to know\n",
"* A matrix is a collection of known numbers _note: there are some cases where you might want to solve for a matrix, but for now let's restrict our use of linear algebra to known matrices_\n",
"\n",
"The specification of _independent equations_ is best illustrated using _dependent_ equations:\n",
"\n",
"1. $x+y = 3$\n",
"\n",
"2. $2x+2y=6$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"ename": "LinAlgError",
"evalue": "Singular matrix",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-5-b58617395425>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mA_sing\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0my_sing\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA_sing\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my_sing\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
"\u001b[0;32m/opt/miniconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m 401\u001b[0m \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'DD->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'dd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 402\u001b[0m \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 403\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 404\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 405\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/opt/miniconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m 95\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 96\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 97\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 98\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 99\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
]
}
],
"source": [
"A_sing = np.array([[1,1],[2,2]])\n",
"y_sing = np.array([3,6])\n",
"np.linalg.solve(A_sing,y_sing)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Singular and ill-conditioned matrices\n",
"\n",
"We receive the `LinAlgError: Singular matrix` because equation (2) is 2 $\\times$ (equation 1). A __singular matrix__ represents a system of equations where multiple solutions exist. \n",
"\n",
"For example, if x=1 and y=2, then\n",
"\n",
"1. $1+2 = 3$\n",
"\n",
"2. $2+4 = 6$\n",
"\n",
"But, we could also say that x=10 and y=-7\n",
"\n",
"1. $10-7 = 3$\n",
"\n",
"2. $20-14 = 6$\n",
"\n",
"Because the system of equations is __singular__, there are an infinite number of values for $x$ and $y$. In engineering, this usually means that we have used the same equation more than once to describe our system of equations, or that we have not properly constrained the problem.\n",
"\n",
"We can also run into matrices that are almost singular, the equations are independent, but they are very close to being dependent. Take the singular example, but let's add a small number, $\\delta$ to one of the constants. \n",
"\n",
"1. $\\delta x+y = 3$\n",
"\n",
"2. $2x+(2+\\delta)y=6$\n",
"\n",
"Now, the equations are independent so we should be able to solve for $x$ and $y$, but depending on the size of $\\delta$ and our machine epsilon (`np.floatinf('float').eps`) we may still have a singular matrix"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"d= 2e-16\n",
"x = 0.0, y = 3.0\n",
"d= 5e-16\n",
"x = -1.0, y = 4.0\n",
"d= 1e-15\n",
"x = -0.4, y = 3.4\n"
]
}
],
"source": [
"for d in [2e-16,5e-16,10e-16]:\n",
" A_ill = np.array([[1+d,1],[2,2]])\n",
" y_ill = np.array([3,6+d])\n",
" print('d=',d)\n",
" x=np.linalg.solve(A_ill,y_ill)\n",
" print('x = {:.1f}, y = {:.1f}'.format(*x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This problem is now __ill-conditioned__. For small changes (as small as the roundoff error in our floating point calculations) we can get different values of $x$ and $y.$ An __ill-conditioned__ matrix (system of equations) is _almost_ worse than the singular matrix because we still get a result, but we wouldn't know the result is so sensitive to our representation of floating point numbers unless we checked the sensitivity of our results. \n",
"\n",
"Luckily, __ill-conditioned__ matrices come up often enough that there is a standard way to recognize the __condition__ of a matrix. The __condition__ of a matrix is the ratio of the measure of its \"magnitude\" compared to the \"magnitude\" of its inverse. Here \"_magnitude_\" is in quotes, because a matrix does not have a single magnitude and direction, instead we call the measure of a matrix a __norm__. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Matrix norms\n",
"\n",
"The Euclidean norm of a vector is measure of the magnitude (in 3D this would be: $|x|=\\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$) in general the equation is:\n",
"\n",
"$||x||_{e}=\\sqrt{\\sum_{i=1}^{n}x_{i}^{2}}$\n",
"\n",
"For a matrix, A, the same norm is called the Frobenius norm:\n",
"\n",
"$||A||_{f}=\\sqrt{\\sum_{i=1}^{n}\\sum_{j=1}^{m}A_{i,j}^{2}}$\n",
"\n",
"In general we can calculate any $p$-norm where\n",
"\n",
"$||A||_{p}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{p}}$\n",
"\n",
"so the p=1, 1-norm is \n",
"\n",
"$||A||_{1}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{1}}=\\sum_{i=1}^{n}\\sum_{i=1}^{m}|A_{i,j}|$\n",
"\n",
"$||A||_{\\infty}=\\sqrt{\\sum_{i=1}^{n}\\sum_{i=1}^{m}A_{i,j}^{\\infty}}=\\max_{1\\le i \\le n}\\sum_{j=1}^{m}|A_{i,j}|$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Condition of a matrix \n",
"### *just checked in to see what condition my condition was in*\n",
"\n",
"The matrix condition is the product of \n",
"\n",
"$Cond(A) = ||A||\\cdot||A^{-1}||$ \n",
"\n",
"So each norm will have a different condition number, but the limit is $Cond(A)\\ge 1$\n",
"\n",
"An estimate of the rounding error is based on the condition of A:\n",
"\n",
"$\\frac{||\\Delta x||}{x} \\le Cond(A) \\frac{||\\Delta A||}{||A||}$\n",
"\n",
"So if the coefficients of A have accuracy to $10^{-t}$\n",
"\n",
"and the condition of A, $Cond(A)=10^{c}$\n",
"\n",
"then the solution for x can have rounding errors up to $10^{c-t}$\n",
"\n",
"The default norm using `np.linalg.norm` is the Frobenius norm, so let's look at the error that can affect our output, $x$ and $y$. \n",
"\n",
"Here, we can say that the accuracy in our constants is $10^{-15}$ so $t=15$ and our norm for $\\delta=10^{-15}$ is\n",
"\n",
"```python\n",
"[in]: print('{:e}'.format(np.linalg.cond(A_ill)))\n",
"[out]:\n",
" 3.045e+15\n",
"```\n",
"\n",
"so $c = 15$. The expected error in $x$ and $y$ is then \n",
"\n",
"$10^{15-15} = 1$\n",
"\n",
"and if we look at the values we calculated when we changed $\\delta$, the values were within a range of $\\approx 1$. \n",
"\n",
"The __key idea__ here is that the condition of the matrix is directly related to the accuracy of your solutions. "
]
},
{
"cell_type": "code",
"execution_count": 108,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.045100e+15\n"
]
}
],
"source": [
"print('{:e}'.format(np.linalg.cond(A_ill)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exercise\n",
"\n",
"One of the most well-known ill-conditioned matrices is the [Hilbert Matrix](https://en.wikipedia.org/wiki/Hilbert_matrix). It is created by placing smaller and smaller fractions in each successive row and column of a matrix. \n",
"\n",
"$H_{ij} = \\frac{1}{i+j-1}$\n",
"\n",
"$H = \\left[\\begin{array}{cccc} \n",
" 1 & \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} \\\\\n",
" \\frac{1}{2} & \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} \\\\\n",
" \\frac{1}{3} & \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} \\\\\n",
" \\frac{1}{4} & \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} \\\\\n",
" \\frac{1}{5} & \\frac{1}{6} & \\frac{1}{7} & \\frac{1}{8} & \\frac{1}{9}\n",
"\\end{array}\\right]$\n",
"\n",
"What is the condition of this $4 \\times 4$ matrix?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# What we've learned\n",
"\n",
"* How to solve a linear algebra problem with `np.linalg.solve`\n",
"* Creating a linear system of equations\n",
"* Identify constants in a linear system $\\mathbf{A}$ and $\\mathbf{b}$\n",
"* Identify unknown variables in a linear system $\\mathbf{x}$\n",
"* Identify a __singular__ or __ill-conditioned__ matrix\n",
"* Calculate the __condition__ of a matrix\n",
"* Estimate the error in the solution based upon the condition of a matrix\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# References\n",
"\n",
"1. Chapra, Steven _Applied Numerical Methods with Matlab for Engineers._ __ch 8.__ McGraw Hill. \n",
"\n",
"2. Kiusalaas, Jaan _Numerical Methods in Engineering with Python 3._ __ch 2.__ Cambridge University Press. \n",
"\n",
"3. [_\"Essence of Linear Algebra\"_](http://3b1b.co/eola) 3 Blue 1 Brown Linear algebra series. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Problems\n",
"\n",
"1. Consider 4 masses connected in series to 4 springs with K=1,000 N/m. What are the final positions of the masses i.e. when acceleration is 0? \n",
"\n",
"![Springs-masses](../images/mass_springs.png)\n",
"\n",
"The masses haves the following amounts, $m_1=1,~m_2=2,~m_3=3,~and~m_4=4 kg.$ Using a FBD for each mass:"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"$m_{1}g+k(x_{2}-x_{1})-kx_{1}=0$\n",
"\n",
"$m_{2}g+k(x_{3}-x_{2})-k(x_{2}-x_{1})=0$\n",
"\n",
"$m_{3}g+k(x_{4}-x_{3})-k(x_{3}-x_{2})=0$\n",
"\n",
"$m_{4}g-k(x_{4}-x_{3})=0$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"in matrix form:\n",
"\n",
"$\\left[ \\begin{array}{cccc}\n",
"2k & -k & 0 & 0 \\\\\n",
"-k & 2k & -k & 0 \\\\\n",
"0 & -k & 2k & -k \\\\\n",
"0 & 0 & -k & k \\end{array} \\right]\n",
"\\left[ \\begin{array}{c}\n",
"x_{1} \\\\\n",
"x_{2} \\\\\n",
"x_{3} \\\\\n",
"x_{4} \\end{array} \\right]=\n",
"\\left[ \\begin{array}{c}\n",
"m_{1}g \\\\\n",
"m_{2}g \\\\\n",
"m_{3}g \\\\\n",
"m_{4}g \\end{array} \\right]$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![HVAC diagram showing the flow rates and connections between floors](../images/hvac.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2. In the figure above you have an idealized Heating, Ventilation and Air conditioning (HVAC) system. In the current configuration, the three-room building is being cooled off by $15^oC$ air fed into the building at 0.1 kg/s. Our goal is to determine the steady-state temperatures of the rooms given the following information\n",
"\n",
"* $\\dot{m}_1=0.1~kg/s$\n",
"* $\\dot{m}_2=0.15~kg/s$\n",
"* $\\dot{m}_3=0.17~kg/s$\n",
"* $\\dot{m}_4=0.15~kg/s$\n",
"* $\\dot{m}_5=0.02~kg/s$\n",
"* $\\dot{m}_6=0.05~kg/s$\n",
"* $C_p=1000~\\frac{J}{kg-K}$\n",
"* $\\dot{Q}_{in} = 500~W$\n",
"* $T_{in} = 15^{o} C$\n",
"\n",
"The energy-balance equations for rooms 1-3 create three equations:\n",
"\n",
"1. $\\dot{m}_1 C_p T_{in}+\\dot{Q}_{in}-\\dot{m}_2 C_p T_{1}+\\dot{m}_6 C_p T_{2} = 0$\n",
"\n",
"2. $\\dot{m}_2 C_p T_{1}+\\dot{Q}_{in}+\\dot{m}_5 C_p T_{3}-\\dot{m}_3 C_p T_{2}-\\dot{m}_6 C_p T_{2} = 0$\n",
"\n",
"3. $\\dot{m}_3 C_p T_{2}+\\dot{Q}_{in}-\\dot{m}_5 C_p T_{3}-\\dot{m}_4 C_p T_{3} = 0$\n",
"\n",
"Identify the unknown variables and constants to create a linear algebra problem in the form of $\\mathbf{Ax}=\\mathbf{b}$.\n",
"\n",
"a. Create the matrix $\\mathbf{A}$\n",
"\n",
"b. Create the known vector $\\mathbf{b}$\n",
"\n",
"c. Solve for the unknown variables, $\\mathbf{x}$\n",
"\n",
"d. What are the warmest and coldest rooms? What are their temperatures?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3. The [Hilbert Matrix](https://en.wikipedia.org/wiki/Hilbert_matrix) has a high condition number and as the matrix increases dimensions, the condition number increases. Find the condition number of a \n",
"\n",
"a. $1 \\times 1$ Hilbert matrix\n",
"\n",
"b. $5 \\times 5$ Hilbert matrix\n",
"\n",
"c. $10 \\times 10$ Hilbert matrix\n",
"\n",
"d. $15 \\times 15$ Hilbert matrix\n",
"\n",
"e. $20 \\times 20$ Hilbert matrix\n",
"\n",
"If the accuracy of each matrix element is $\\approx 10^{-16}$, what is the expected rounding error in the solution $\\mathbf{Ax} = \\mathbf{b}$, where $\\mathbf{A}$ is the Hilbert matrix. \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}