Skip to content

Instantly share code, notes, and snippets.

@tobydriscoll

tobydriscoll/TB-Lecture-20-julia.ipynb

Last active October 22, 2016 15:37
Show Gist options
  • Save tobydriscoll/e2f44e7756b864050b52c6773aa83b70 to your computer and use it in GitHub Desktop.
Save tobydriscoll/e2f44e7756b864050b52c6773aa83b70 to your computer and use it in GitHub Desktop.
Download ZIP
TB Lecture 20
Display the source blob
Display the rendered blob
Raw
TB-Lecture-20-julia.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
TB-Lecture-20-matlab.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lecture 20: Gaussian elimination\n",
"\n",
"To solve the system $Ax=b$, you are taught a method of \"triangular triangularization.\" For example..."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A =\n",
"\n",
" 17 2 3 13\n",
" 5 12 10 8\n",
" 9 7 7 12\n",
" 4 14 15 2\n"
]
}
],
"source": [
"A = magic(4) + eye(4)\n",
"x = (4:-1:1)'; % exact answer\n",
"b = A*x;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Start by introducing zeros below the diagonal in the first column, using row operations."
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AA =\n",
"\n",
" 17 2 3 13 93 \n",
" 0 194/17 155/17 71/17 963/17 \n",
" 0 101/17 92/17 87/17 574/17 \n",
" 0 230/17 243/17 -18/17 1158/17\n"
]
}
],
"source": [
"format rat\n",
"AA = [A b];\n",
"AA(2,:) = AA(2,:) - (5/17)*AA(1,:);\n",
"AA(3,:) = AA(3,:) - (9/17)*AA(1,:);\n",
"AA(4,:) = AA(4,:) - (4/17)*AA(1,:)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's pause to express these operations using matrix algebra. Each row operation is realized by performing the same operation to an identity matrix, resulting in\n",
"\n",
"$$L_{i1} = I + \\mu_i e_i e_1^*,$$\n",
"\n",
"where $\\mu_i$ is a multiplier and $e_j$ is the $j$th column of $I$."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 17 2 3 13 93 \n",
" 0 194/17 155/17 71/17 963/17 \n",
" 0 101/17 92/17 87/17 574/17 \n",
" 0 230/17 243/17 -18/17 1158/17\n"
]
}
],
"source": [
"I = eye(4);\n",
"L21 = I + (-5/17)*I(:,2)*I(:,1)';\n",
"L31 = I + (-9/17)*I(:,3)*I(:,1)';\n",
"L41 = I + (-4/17)*I(:,4)*I(:,1)';\n",
"\n",
"L41*(L31*(L21*[A b]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But it's clear from a little algebra that the product of the $L$ factors is very simple:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L1 =\n",
"\n",
" 1 0 0 0 \n",
" -5/17 1 0 0 \n",
" -9/17 0 1 0 \n",
" -4/17 0 0 1\n"
]
}
],
"source": [
"L1 = L41*L31*L21"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In words, we just accumulate all the mulitpliers for column 1 into a single matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's rewind to the start. This time, we'll ignore the $b$ vector, because we can always apply the $L$ factors to it later and get the same effect."
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A1 =\n",
"\n",
" 17 2 3 13 \n",
" 0 194/17 155/17 71/17 \n",
" 0 101/17 92/17 87/17 \n",
" 0 230/17 243/17 -18/17\n"
]
}
],
"source": [
"A1 = L1*A"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we eliminate in the second column, below the diagonal."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A2 =\n",
"\n",
" 17 2 3 13 \n",
" 0 194/17 155/17 71/17 \n",
" 0 0 129/194 571/194 \n",
" 0 0 338/97 -583/97\n"
]
}
],
"source": [
"L2 = I;\n",
"L2(3,2) = (-101/194);\n",
"L2(4,2) = (-230/194);\n",
"A2 = L2*A1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we have one last operation in column 3. (I've also reached the end of my patience with the rational format.)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A3 =\n",
"\n",
" 17.0000 2.0000 3.0000 13.0000\n",
" 0 11.4118 9.1176 4.1765\n",
" 0 0 0.6649 2.9433\n",
" 0 0 0 -21.4341\n"
]
}
],
"source": [
"format\n",
"L3 = I;\n",
"L3(4,3) = -A2(4,3)/A2(3,3);\n",
"A3 = L3*A2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In fact, we should have called this matrix $U$, for upper triangular. In summary:"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U =\n",
"\n",
" 17.0000 2.0000 3.0000 13.0000\n",
" 0 11.4118 9.1176 4.1765\n",
" 0 0 0.6649 2.9433\n",
" 0 0 0 -21.4341\n"
]
}
],
"source": [
"U = L3*(L2*(L1*A))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We resort to our favorite maneuver: reassociate matrix terms."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"M =\n",
"\n",
" 1.0000 0 0 0\n",
" -0.2941 1.0000 0 0\n",
" -0.3763 -0.5206 1.0000 0\n",
" 2.0853 1.5426 -5.2403 1.0000\n"
]
}
],
"source": [
"M = L3*L2*L1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a _unit_ lower triangular matrix, with ones on the diagonal. Those terms below the diagonal are not the same individual row multipliers from before:"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 1.0000 0 0\n",
" -0.2941 1.0000 0\n",
" -0.5294 -0.5206 1.0000\n",
" -0.2353 -1.1856 -5.2403\n"
]
}
],
"source": [
"[ L1(:,1) L2(:,2) L3(:,3) ]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, something amazing happens if we take the _inverse_ of the $M$ matrix above."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L =\n",
"\n",
" 1.0000 0 0 0\n",
" 0.2941 1.0000 -0.0000 -0.0000\n",
" 0.5294 0.5206 1.0000 -0.0000\n",
" 0.2353 1.1856 5.2403 1.0000\n"
]
}
],
"source": [
"L = inv(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The multipliers are back! So instead of having $U=MA$, we'll write $A=LU$ for a unit lower triangular $L$ whose elements are the negatives of the row multipliers."
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans =\n",
"\n",
" 1.0e-14 *\n",
"\n",
" 0 0 0 0\n",
" 0 0 0 0\n",
" -0.1776 -0.1776 -0.1776 -0.1776\n",
" -0.2665 -0.3553 -0.5329 -0.3553\n"
]
}
],
"source": [
"A - L*U"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Matlab",
"language": "matlab",
"name": "matlab"
},
"language_info": {
"codemirror_mode": "octave",
"file_extension": ".m",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
}
],
"mimetype": "text/x-octave",
"name": "matlab",
"version": "0.11.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment