{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(section:linear-equations-row-reduction)=\n",
"# Row Reduction (Gaussian Elimination)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Last revised on August 25, 2025**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**References:**\n",
"\n",
"- {cite}`Chasnov` Section 3.1, *Gaussian Elimination*.\n",
"\n",
"- {cite}`Sauer` Section 2.1.1, *Naive Gaussian Elimination*.\n",
"\n",
"- {cite}`Burden-Faires` Section 6.1, *Linear Systems of Equations*.\n",
"\n",
"- {cite}`Dionne` Section 4.1, *Gaussian Elimination with Backward Substitution*.\n",
"This includes partial pivoting, as seen here in the section {ref}`section:linear-equations-pivoting`.\n",
"\n",
"- {cite}`Chenney-Kincaid` Section 2.1, *Naive Gaussian Elimination*.\n",
"\n",
"- {cite}`Trefethen-Bau` Lecture 20, *Gaussian Elimination*. (This starts directly with the LU factorization approach as in section {ref}`section:linear-equations-lu-factorization` — that is not a bad idea!).\n",
"\n",
"- {cite}`Dahlquist-Bjorck-v2` Section 7.2.1, *Triangular Systems and Gaussian Elimination*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"This section introduces a basic algorithm for solving a system of $n$ linear equations in $n$ unknowns,\n",
"based on the strategy of row-reduction or Gaussian elimination.\n",
"You might have seen this strategy in a linear algebra course, but if not, do not worry:\n",
"this chapter introduces everything from the beginning, not assuming any previous knowledge of linear algebra beyond matrix notation.\n",
"\n",
"This section also takes the opportunity to introduce some general tactics for getting from mathematical facts to precise algorithms and then to code which uses these algorithms."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Strategy for getting from mathematical facts to a good algorithm and then to its implentation in [Julia] code\n",
"\n",
"Here I take the opportunity to illustrate some useful strategies for getting from mathematical facts and ideas to good algorithms and working code for solving a numerical problem.\n",
"The pattern we will see here, and often later, is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(subsubsection:algorithm-strategy-basic-algorithm)=\n",
"### Step 1. Get a basic algorithm:\n",
"1. Start with mathematical facts (like the equations $\\sum_{j=1}^n a_{ij}x_j = b_i$).\n",
"2. Solve to get an equation for each unknown — or for an updated aproximation of each unknown — in terms of other quantitities.\n",
"3. Specify an order of evaluation in which all the quantities at right are evaluated earlier.\n",
"\n",
"In this, it is often best to start with a verbal description before specifying the details in more precise and detailed mathematical form."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Step 2. Refine to get a more **robust** algorithm:\n",
"1. Identify cases that can lead to failure due to division by zero and such, and revise to avoid them.\n",
"2. Avoid inaccuracy due to problems like severe rounding error. One rule of thumb is that anywhere that a zero value is a fatal flaw (in particular, division by zero), a very small value is also a hazard when rounding error is present.\n",
"So **avoid very small denominators**. (We will soon examine this through the phenomenon of **loss of significance** and it extreme case **catastrophic cancellation**.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(subsubsection:algorithm-strategy-robustness)=\n",
"### Step 2. Refine to get a more **robust** algorithm:\n",
"1. Identify cases that can lead to failure due to division by zero and such, and revise to avoid them.\n",
"2. Avoid inaccuracy due to problems like severe rounding error. One rule of thumb is that anywhere that a zero value is a fatal flaw (in particular, division by zero), a very small value is also a hazard when rounding error is present.\n",
"So *avoid very small denominators*. (The section on\n",
"[Machine numbers, rounding error and Error Propagation](section:machine-numbers-rounding-error-error-propagation)\n",
"examines this through the phenomenon of **loss of significance**.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(subsubsection:algorithm-strategy-efficiency)=\n",
"### Step 3. Refine to get a more **efficient** algorithm\n",
"For example,\n",
"- Avoid repeated evaluation of exactly the same quantity.\n",
"- Avoid redundant calculations, such as ones whose value can be determnied in advance;\n",
"for example, values that can be shown in advance to be zero.\n",
"- Compare and choose between alternative algorithms."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Row Reduction, a.k.a. Gaussian Elimination\n",
"\n",
"We start by considering the most basic algorithm, based on ideas seen in a linear algebra course.\n",
"\n",
"The problem is best stated as a collection of equations for individual numerical values:\n",
"\n",
"Given coefficients $a_{i,j} 1 \\leq i \\leq n,\\, 1 \\leq j \\leq n$ and right-hand side values $b_i,\\, 1 \\leq i \\leq n$,\n",
"solve for the $n$ unknowns $x_j,\\, 1 \\leq j \\leq n$ in the equations\n",
"\n",
"$$\n",
"\\sum_{j=1}^n a_{i,j} x_j = b_i,\\, 1 \\leq i \\leq n.\n",
"$$\n",
"\n",
"In verbal form, the basic strategy of *row reduction* or *Gaussian elimination* is this:\n",
"\n",
"- **Choose** one equation and use it to eliminate one **chosen** unknown from all the other equations, leaving that chosen equation plus $n-1$ equations in $n-1$ unknowns.\n",
"- Repeat recursively, at each stage using one of the remaining equations to eliminate one of the remaining unknowns from all the other equations.\n",
"- This gives a final equation in just one unknown, preceeded by an equation in that unknown plus one other, and so on: solve them in this order, from last to first."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Determining those choices, to produce a first algorithm: \"Naive Row Reduction\"\n",
"\n",
"A precise algorithm must include rules specifying all the choices indicated above.\n",
"The simplest \"naive\" choice, which works in most but not all cases, is to eliminate from the top to bottom and left to right:\n",
"\n",
"- Use the first equation to eliminate the first unknown from all other equations.\n",
"- Repeat recursively, at each stage using the first remaining equation to eliminate the first remaining unknown. Thus, at step $k$, equation $k$ is used to eliminate unknown $x_k$.\n",
"- This gives one equation in just the last unknown $x_n$; another equation in the last two unknowns $x_{n-1}$ and $x_n$, and so on: solve them in this reverse order, evaluating the unknowns from last to first."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This usually works, but can fail because at some stage the (updated) $k$-th equation might not include the $k$-th unknown: that is, its coefficient might be zero, leading to division by zero.\n",
"\n",
"We will refine the algorithm to deal with that in the later section {ref}`section:linear-equations-pivoting`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The general case of solving $Ax = b$\n",
"\n",
"The problem of solving $Ax = b$ in general, when all you know is that $A$ is an $n \\times n$ matrix and $b$ is an $n$-vector, can in most cases be handled well by using standard software rather than by writing your own code. Here is an example in Julia, solving\n",
"\n",
"$$\n",
"\\left[ \\begin{array}{rrr} 4 & 2 & 7 \\\\ 3 & 5 & -6 \\\\ 1 & -3 & 2 \\end{array} \\right]\n",
"\\left[ \\begin{array}{r} x_1 \\\\ x_2 \\\\ x_3 \\end{array} \\right]\n",
"= \\left[ \\begin{array}{r} 2 \\\\ 3 \\\\ 4 \\end{array} \\right]\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A =\n",
"[4.0 2.0 7.0; 3.0 5.0 -6.0; 1.0 -3.0 2.0]\n",
"b = [2.0, 3.0, 4.0]\n"
]
}
],
"source": [
"A = [4.0 2.0 7.0; 3.0 5.0 -6.0; 1.0 -3.0 2.0]\n",
"println(\"A =\\n$(A)\")\n",
"b = [2.0; 3.0; 4.0]\n",
"println(\"b = $(b)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} On Julia\n",
":label: remark-julia-arrays\n",
"\n",
"- See the notes on {ref}`arrays` in {ref}`section:julia-language-notes`.\n",
"\n",
"- Julia mimics Matlab's notation for \"dividing from the left\":\n",
"the solution of $Ax = b$ is $x = A^{-1} b$ and given by `A\\b`; it is not $b A^{-1}$ which is what you get from the usual \"divide from the right\" notation of `b/A`.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia says that the solution of Ax=b is x=[1.8116883116883116, -1.0324675324675323, -0.45454545454545453]\n"
]
}
],
"source": [
"x = A\\b;\n",
"println(\"Julia says that the solution of Ax=b is x=$(x)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check the **residual** $b - Ax$, a measure of *backward error*:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"r = b-A*x;"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"As a check, the residual is\n",
" r = b - Ax = [0.0, 0.0, 8.881784197001252e-16]\n",
"and its infinity (or 'maximum') norm is\n",
" ||r|| = 8.881784197001252e-16\n"
]
}
],
"source": [
"println()\n",
"println(\"As a check, the residual is\")\n",
"println(\" r = b - Ax = $(r)\")\n",
"println(\"and its infinity (or 'maximum') norm is\")\n",
"println(\" ||r|| = $(maximum(abs.(r)))\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Not quite zero values and rounding\n",
":label: remark-1-not-quite-zero-values-and-rounding\n",
"\n",
"Some values here that you might hope to be zero are instead very small non-zero numbers, with exponent $10^{-16}$,\n",
"due to rounding error in computer arithmetic.\n",
"For details on this (like why \"-16\" in particular) see {ref}`section:machine-numbers-rounding-error-error-propagation`.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The naive row reduction algorithm, in pseudo-code\n",
"\n",
"Here the elements of the transformed matrix and vector after step $k$ are named $a_{i,j}^{(k)}$ and $b_{k}^{(k)}$, so that the original values are $a_{i,j}^{(0)} = a_{i,j}$ and $b_{i}^{(0)} = b_{i}$.\n",
"\n",
"The name $l_{i,k}$ is given to the multiple of row $k$ that is subtracted from row $i$ at step $k$.\n",
"This naming might seem redundant, but it becomes very useful later,\n",
"in the section on [LU factorization](section:linear-equations-lu-factorization)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:algorithm} naive row reduction\n",
":label: naive-row-reduction \n",
"\n",
"for k from 1 to n-1 $\\qquad$ *Step k: get zeros in column k below row k:*\n",
"
\n",
"$\\quad$ for i from k+1 to n\n",
"
\n",
"$\\qquad$ *Evaluate the multiple of row k to subtract from row i:*\n",
"
\n",
"$\\quad\\quad l_{i,k} = a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}$ $\\qquad$ **If** $a_{k,k}^{(k-1)} \\neq 0$!\n",
"
\n",
"$\\qquad$ *Subtract $(l_{i,k}$ times row k) from row i in matrix A ...:*\n",
"
\n",
"$\\quad\\quad$ for j from 1 to n\n",
"
\n",
"$\\quad\\quad\\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)} - l_{i,k} a_{k,j}^{(k-1)}$\n",
"
\n",
"$\\quad\\quad$ end\n",
"
\n",
"$\\qquad$ ... and at right, subtract $(l_{i,k}$ times $b_k)$ from $b_i$:\n",
"
\n",
"$\\quad\\quad b_i^{(k)} = b_i^{(k-1)} - l_{i,k} b_{k}^{(k-1)}$ \n",
"
\n",
"$\\quad$ end\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The rows before $i=k$ are unchanged, so they are ommited from the update;\n",
"however, in a situation where we need to complete the definitions of $A^{(k)}$ and $b^{(k)}$ we would also need the following inside the `for k` loop:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:algorithm} Inserting the zeros below the main diagonal\n",
":label: row-reduction -inserting-zeros\n",
"\n",
"$\\quad$ for i from 1 to k\n",
"
\n",
"$\\quad\\quad$ for j from 1 to n\n",
"
\n",
"$\\quad\\quad\\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)}$\n",
"
\n",
"$\\quad\\quad$ end\n",
"
\n",
"$\\quad\\quad b_i^{(k)} = b_i^{(k-1)}$\n",
"
\n",
"$\\quad$ end\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, the algorithm will usually be implemented by overwriting the previous values in an array with new ones, and then this part is redundant."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next improvement in efficiency: the updates in the first $k$ columns at step $k$ give zero values (that is the key idea of the algorithm!),\n",
"so there is no need to compute or store those zeros,\n",
"and thus the only calculations needed in the above `for j from 1 to n` loop are covered by `for j from k+1 to n`.\n",
"Thus from now on we use only the latter—except when, for demonstration purposes, we need those zeros.\n",
"\n",
"Thus, the standard algorithm looks like this:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:algorithm} basic row-reduction \n",
":label: row-reduction \n",
"\n",
"for k from 1 to n-1 $\\qquad$ *Step k: Get zeros in column k below row k:*\n",
"
\n",
"$\\quad$ for i from k+1 to n $\\qquad$ *Update only the rows that change: from k+1 on:*\n",
"
\n",
"$\\qquad$ *Evaluate the multiple of row k to subtract from row i:*\n",
"
\n",
"$\\quad\\quad l_{i,k} = a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}$ $\\qquad$ **If** $a_{k,k}^{(k-1)} \\neq 0$!\n",
"
\n",
"$\\qquad$ *Subtract $(l_{i,k}$ times row k) from row i in matrix A, in the columns that are not automaticaly zero:*\n",
"
\n",
"$\\quad\\quad$ for j from k+1 to n\n",
"
\n",
"$\\quad\\quad\\quad a_{i,j}^{(k)} = a_{i,j}^{(k-1)} - l_{i,k} a_{k,j}^{(k-1)}$\n",
"
\n",
"$\\quad\\quad$ end\n",
"
\n",
"$\\qquad$ *and at right, subtract $(l_{i,k}$ times $b_k)$ from $b_i$:*\n",
"
\n",
"$\\quad\\quad b_i^{(k)} = b_i^{(k-1)} - l_{i,k} b_{k}^{(k-1)}$ \n",
"
\n",
"$\\quad$ end\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The naive row reduction elimination algorithm, in Julia\n",
"\n",
"Conversion to actual Julia code is now quite straightforward; there is litle more to be done than:\n",
"\n",
"- Change the way that indices are described, from $b_i$ to `b[i]` and from $a_{i,j}$ to `A[i,j]`.\n",
"\n",
"- Use case consistently in array names, since the quirk in mathematical notation of using upper-case letters for matrix names but lower case letters for their elements is gone!\n",
"In these notes, matrix names will be upper-case and vector names will be lower-case (even when a vector is considered as 1-column matrix).\n",
"\n",
"- Rather than create a new array for each matrix $A^{(0)}$, $A^{(1)}$, etc. and each vector $b^{(0)}$, $b^{(1)}$,\n",
"we overwite each in the same array."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" for k in 1:n\n",
" for i in k+1:n\n",
" L[i,k] = A[i,k] / A[k,k]\n",
" for j in k+1:n\n",
" A[i,j] -= L[i,k] * A[k,j]\n",
" end\n",
" b[i] -= L[i,k] * b[k]\n",
" end\n",
" end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To demonstrate this, some additions are needed:\n",
"- Putting this algorithm into a function.\n",
"- Getting the value $n$ needed for the loop, using the fact that it is the length of vector `b`.\n",
"- Creating the array $L$.\n",
"- Copying the input arrays `A` and `b` into new ones, `U` and `c`, so that the original arrays are not changed. That is, when the row reduction is completed, `U` contains $A^{(n-1)}$ and `c` contains $b^{(n-1)}$.\n",
"\n",
"Also, for some demonstrations, the zero values below the main diagonal of `U` are inserted, though usually they would not be needed."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"function row_reduce(A, b)\n",
" # To avoid modifying the matrix and vector specified as input,\n",
" # they are copied to new arrays, with the function copy().\n",
" # Warning: it does not work to say \"U = A\" and \"c = b\";\n",
" # this makes these names synonyms, referring to the same stored data.\n",
"\n",
" U = copy(A) # not \"U=A\", which makes U and A synonyms\n",
" c = copy(b)\n",
" n = length(b)\n",
" L = zeros(n, n)\n",
" for k in 1:n-1\n",
" for i in k+1:n\n",
" # compute all the L values for column k:\n",
" L[i,k] = U[i,k] / U[k,k] # Beware the case where U[k,k] is 0\n",
" for j in k+1:n\n",
" U[i,j] -= L[i,k] * U[k,j]\n",
" end\n",
" # Put in the zeros below the main diagonal in column k of U;\n",
" # this is not important for calculations,\n",
" # since those elements of U are not used in backward substitution,\n",
" # but it helps for displaying results and for checking the results via residuals.\n",
" U[i,k] = 0.\n",
" \n",
" c[i] -= L[i,k] * c[k]\n",
" end\n",
" end\n",
" for i in 2:n\n",
" for j in 1:i-1\n",
" U[i,j] = 0.\n",
" end\n",
" end\n",
" return (U, c)\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is a helper function for displaying matrices.\n",
"Since it will be used in several sections, it is also in the module `NumericalMethods`, along with the above function;\n",
"see {ref}`module:NumericalMethods`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Julia Modules\n",
":label: julia-modules\n",
"\n",
"Modules and their usage are introduced in the notes on [Using modules and packages](modules).\n",
"As explained there, these two functions can be obtained with\n",
"\n",
" include(\"NumericalMethods.jl\")\n",
" using .NumericalMethods:row_reduce\n",
" using .NumericalMethods:print_matrix\n",
"\n",
"Details about creating your own modules will be given in [creating modules](creating-modules-future) once those notes are written;\n",
"meanwhile, examining {ref}`module:NumericalMethods` could help:\n",
"the actual module definition file is just the Julia code from there with the interspersed explanatory comments removed.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" Some helper functions to cleanup output"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"function print_matrix(A)\n",
" # A helper function to \"pretty print\" matrices\n",
" \n",
" (rows, cols) = size(A)\n",
" print(\"[ \")\n",
" for row in 1:rows\n",
" if row > 1\n",
" print(\" \")\n",
" end\n",
" for col in 1:cols\n",
" print(A[row,col], \" \")\n",
" end\n",
" if row < rows;\n",
" println()\n",
" else\n",
" println(\"]\")\n",
" end\n",
" end\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [],
"source": [
"# A shorthand for rounding off to n significant digits.\n",
"import Base: round\n",
"round(x, n::Integer) = round(x, sigdigits=n);"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A is\n",
"[ 4.0 2.0 7.0 \n",
" 3.0 5.0 -6.0 \n",
" 1.0 -3.0 2.0 ]\n",
"b = [2.0, 3.0, 4.0]\n"
]
}
],
"source": [
"println(\"A is\")\n",
"print_matrix(A)\n",
"println(\"b = $(b)\")"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"(U, c) =row_reduce(A, b);"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Row reduction gives\n",
"U =\n",
"[ 4.0 2.0 7.0 \n",
" 0.0 3.5 -11.25 \n",
" 0.0 0.0 -11.0 ]\n",
"c = [2.0, 1.5, 5.0]\n"
]
}
],
"source": [
"println(\"Row reduction gives\")\n",
"println(\"U =\")\n",
"print_matrix(U)\n",
"println(\"c = $(c)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take advantage of the fact that we have used Julia's built-in linear algebra command `b\\A` to get a very accurate approximation of the solution $x$ to $Ax=b$; this should also solve $Ux=c$, so check the backward error, a.k.a. the *residual*:"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The residual (backward error) r = c-Ux is [0.0, -2.22e-16, 0.0], with maximum norm 2.22e-16\n"
]
}
],
"source": [
"r = c - U * x\n",
"println(\"The residual (backward error) r = c-Ux is $(round.(r,4)),\" *\n",
" \" with maximum norm $(round(maximum(abs.(r)),4))\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Array slicing in Julia\n",
":label: julia-array-slicing\n",
"\n",
"Many operations in linear algebra can be expressed more concisely using\n",
"[array slicing](array-indexing-slicing) and\n",
"[vectorization](vectorization-and-broadcasting) as described in {ref}`section:julia-language-notes`\n",
"allowing the loops above to be expressed as\n",
"\n",
" for k in 1:n\n",
" L[k+1:end,end] = A[k+1:end,k] / A[k,k]\n",
" A[k+1:end,k+1:end] -= L[k+1:end,k] * A[[k],k+1:end]\n",
" b[k+1:end] -= L[k+1:end,k] * b[k]\n",
" end\n",
" end\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Matrix slicing in Julia\n",
":label: julia-array-slicing-2\n",
"\n",
"One subtlety here, as mentioned in the notes on [slicing](array-indexing-slicing):\n",
"that row slice at the end of line 3 has to be done as `A[[k],k+1:end]` with brackets around the row index,\n",
"in order to make it a 1-row matrix;\n",
"using `A[k,k+1:end]` instead would give a vector, and then the product would fail.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I will break my usual guideline of non-repetition by redefining `row_reduce`,\n",
"since this is just a restatement of exactly the same algorithm with different Julia notation.\n",
"\n",
"While I am about it, I add a `demo_mode`, for display of intermediate results."
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {},
"outputs": [],
"source": [
"function row_reduce(A, b; demo_mode=false)\n",
" # To avoid modifying the matrix and vector specified as input,\n",
" # they are copied to new arrays, with the method .copy()\n",
" # Warning: it does not work to say \"U = A\" and \"c = b\";\n",
" # this makes these names synonyms, referring to the same stored data.\n",
" \n",
" # This version vectorizes the inner loops, and all of the \"i, j\" loop for updating U.\n",
"\n",
" U = copy(A) # not \"U=A\", which makes U and A synonyms\n",
" c = copy(b)\n",
" n = length(b)\n",
" L = zeros(n, n)\n",
" for k in 1:n-1\n",
" if demo_mode; println(\"Step $k:\"); end\n",
" # compute all the L values for column k:\n",
" L[k+1:end,k] = U[k+1:end,k] / U[k,k] # Beware the case where U[k,k] is 0\n",
" if demo_mode\n",
" println(\"The multipliers in column $(k) are $(L[k+1:end,k])\")\n",
" end\n",
" U[k+1:end,k+1:end] -= L[k+1:end,k] * U[[k],k+1:end]\n",
" c[k+1:end] -= L[k+1:end,k] * c[k]\n",
" \n",
" # Insert the below-diagonal zeros in column k;\n",
" # this is not important for calculations,\n",
" # since those elements of U are not used in backward substitution,\n",
" # but it helps for displaying results and for checking the results via residuals.\n",
" U[k+1:end,k] .= 0.0\n",
"\n",
" if demo_mode\n",
" println(\"The updated matrix is\")\n",
" print_matrix(U)\n",
" println(\"The updated right-hand side is $c\")\n",
" end\n",
" end\n",
" return (U, c)\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Repeating the above testing:"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Step 1:\n",
"The multipliers in column 1 are [0.75, 0.25]\n",
"The updated matrix is\n",
"[ 4.0 2.0 7.0 \n",
" 0.0 3.5 -11.25 \n",
" 0.0 -3.5 0.25 ]\n",
"The updated right-hand side is [2.0, 1.5, 3.5]\n",
"Step 2:\n",
"The multipliers in column 2 are [-1.0]\n",
"The updated matrix is\n",
"[ 4.0 2.0 7.0 \n",
" 0.0 3.5 -11.25 \n",
" 0.0 0.0 -11.0 ]\n",
"The updated right-hand side is [2.0, 1.5, 5.0]\n",
"Row reduction gives U=\n",
"[ 4.0 2.0 7.0 \n",
" 0.0 3.5 -11.25 \n",
" 0.0 0.0 -11.0 ]\n",
"and right-hand side [2.0, 1.5, 5.0]\n"
]
}
],
"source": [
"(U, c) =row_reduce(A, b,demo_mode=true);\n",
"println(\"Row reduction gives U=\")\n",
"print_matrix(U)\n",
"println(\"and right-hand side $c\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Backward substitution with an upper triangular matrix\n",
"\n",
"The transformed equations have the form\n",
"\n",
"$$\n",
"\\begin{split}\n",
"u_{1,1} x_1 + u_{1,2} x_2 + u_{1,3} x_3 + \\cdots + u_{1,n} x_n &= c_1 \\\\\n",
"\\vdots \\\\\n",
"u_{i,i} x_i + u_{i+1,i+1} x_{i+1} + \\cdots + u_{i,n} x_n &= c_i \\\\\n",
"\\vdots \\\\\n",
"u_{n-1,n-1} x_{n-1} + u_{n-1,n} x_{n} &= c_{n-1} \\\\\n",
"u_{nn} x_n &= c_n \\\\\n",
"\\end{split}\n",
"$$\n",
"\n",
"and can be solved from bottom up, starting with $x_n = c_n/u_{n,n}$.\n",
"\n",
"All but the last equation can be written as\n",
"\n",
"$$\n",
"u_{i,i}x_i + \\sum_{j=i+1}^{n} u_{i,j} x_j = c_i, \\; 1 \\leq i \\leq n-1\n",
"$$\n",
"\n",
"and so solved as\n",
"\n",
"$$\n",
"x_i = \\frac{c_i - \\sum_{j=i+1}^{n} u_{i,j} x_j}{u_{i,i}},\n",
"\\qquad \\textbf{ If } u_{i,i} \\neq 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This procedure is **backward substitution**, giving the algorithm\n",
"```{prf:algorithm} Backward substitution\n",
":label: backward-substitution\n",
"\n",
"$x_n = c_n/u_{n,n}$\n",
"
for i from n-1 down to 1\n",
"
$\\displaystyle \\quad x_i = \\frac{c_i - \\sum_{j=i+1}^{n} u_{i,j} x_j}{u_{i,i}}$\n",
"
end\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This works so long as none of the main diagonal terms $u_{i,i}$ is zero, because when done in this order, everything on the right hand side is known by the time it is evaluated.\n",
"\n",
"For future reference, note that the elements $u_{k,k}$ that must be non-zero here, the ones on the **main diagonal** of $U$, are the same as the elements $a_{k,k}^{(k)}$ that must be non-zero in the row reduction stage above, because after stage $k$, the elements of row $k$ do not change any more: $a_{k,k}^{(k)} = a_{k,k}^{(n-1)} = u_{k,k}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Summing in Julia\n",
":label: julia-summing\n",
"\n",
"For an $n$-element single-index array `v`, the sum of its elements $\\sum_{i=1}^{n} v_i$ is given by `sum(v)`.\n",
"\n",
"Thus $\\sum_{i=a}^{b} v_i$, the sum over a subset of indices $[a,b]$, is given by `sum(v[a:b])`.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The backward substitution algorithm in Julia\n",
"\n",
"With all the above Julia details, the core code for backward substitution is:\n",
"\n",
" x[end] = c[end]/U[end,end]\n",
" for i in n-1:-1:1\n",
" x[i] = (c[i] - U[i,i+1:end] * x[i+1:end])) / U[i,i]\n",
" end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:observation} \n",
":label: mathematically-correct-notation\n",
"\n",
"Note that the backward substitution algorithm and its Julia coding have a nice mathematical advantage over the row reduction algorithm above:\n",
"the precise mathematical statement of the algorithm does not need any intermediate quantities distinguished by superscripts ${}^{(k)}$\n",
"and correspondingly, all variables in the code have fixed meanings, rather than changing at each step.\n",
"\n",
"In other words, all uses of the equal sign are mathematically correct as equations!\n",
"\n",
"This can be advantageous in creating algorithms and code that is more understandable and more readily verified to be correct,\n",
"and is an aspect of the *functional programming* approach.\n",
"We will soon go part way to that *functional* ideal, by rephrasing row-reduction in a form where all variables have clear,\n",
"fixed meanings, corresponding to the natural mathematical description of the process:\n",
"the method of **LU factorization** introduced in {ref}`section:linear-equations-lu-factorization`.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a final demonstration, we put this second version of the code into a complete working Julia function and test it:"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [],
"source": [
"function backward_substitution(U, c;demo_mode=false)\n",
" n = length(c)\n",
" x = zeros(n)\n",
" x[end] = c[end]/U[end,end]\n",
" if demo_mode\n",
" println(\"x_$n = $(x[n])\")\n",
" end\n",
" for i in n-1:-1:1\n",
" x[i] = ( c[i] - sum(U[i,i+1:end] .* x[i+1:end]) ) / U[i,i]\n",
" if demo_mode\n",
" println(\"x_$i = $(x[i])\")\n",
" end\n",
" end\n",
" return x\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which, as usual, is also available via\n",
"\n",
" using .NumericalMethods:backward_substitution"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x_3 = -0.45454545454545453\n",
"x_2 = -1.0324675324675323\n",
"x_1 = 1.8116883116883116\n",
" \n",
"x = [1.8116883116883116, -1.0324675324675323, -0.45454545454545453]\n",
"The residual b - Ax = [0.0, 0.0, 8.882e-16], with maximum norm 8.882e-16)\n"
]
}
],
"source": [
"x = backward_substitution(U, c,demo_mode=true)\n",
"println(\" \")\n",
"println(\"x = $x\")\n",
"r = b - A*x\n",
"println(\"The residual b - Ax = $(round.(r,4)), with maximum norm $(round(maximum(abs.(r)),4)))\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since one is often just interested in the solution given by the two steps of row reduction and then backward substitution,\n",
"they can be combined in a single function by composition:"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [],
"source": [
"solve_linear_system(A, b) = backward_substitution(row_reduce(A, b)...);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} The Julia \"splat\" operation\n",
":label: julia-splat\n",
"\n",
"The **splat** notation \"...\" takes the item to its left (here the tuple `(U, c)` returned by `row_reduce`) and unpacks it into separate items [here `U` and `c`, as needed for input to `backwardsubstitution`.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x = [1.8116883116883116, -1.0324675324675323, -0.45454545454545453]\n"
]
}
],
"source": [
"x = solve_linear_system(A, b)\n",
"println(\"x = $x\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Two code testing hacks: starting from a known solution, and using randomly generated examples\n",
"\n",
"An often useful strategy in developing and testing code is to create a test case with a known solution;\n",
"another is to use random numbers to avoid accidently using a test case that in unusually easy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} function rand! from module Random\n",
":label: julia-Random-rand\n",
"\n",
"- The preferred style is to have all `import` and `using` statements near the top,\n",
"but since this is the first time we've heard of module `Random` I did not want it to be mentioned mysteriously above.\n",
"\n",
"- The exclamation point in `rand!` indicates that the function modifies its input.\n",
"This will be discussed in [Functions part 2](functions2-future) when that section gets written.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [],
"source": [
"using Random: rand!"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x_random = [0.09041138353921396, -0.5531448066870839, 0.4495100318153302]\n"
]
}
],
"source": [
"n = length(b)\n",
"x_random = zeros(n)\n",
"rand!(x_random) # fill with random values, uniform in [0, 1)\n",
"(x_random *= 2.0) .-= 1.0 # double and then subtract 1: now uniform in [-1, 1)\n",
"println(\"x_random = $x_random\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a right-hand side b that automatically makes `x_random` the correct solution:"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"b_random = A * x_random;"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A is\n",
"[ 4.0 2.0 7.0 \n",
" 3.0 5.0 -6.0 \n",
" 1.0 -3.0 2.0 ]\n",
"\n",
"brandom is [2.40192614349, -5.1915500737097595, 2.648865867231126]\n",
"\n",
"U is\n",
"[ 4.0 2.0 7.0 \n",
" 0.0 3.5 -11.25 \n",
" 0.0 0.0 -11.0 ]\n",
"\n",
"Residualc_random - U*x_random = [0.0, 0.0, 0.0]\n",
"\n",
"x_computed is [0.0904113835392, -0.553144806687, 0.449510031815]\n",
"\n",
"Residual b_random-A*x_computed is [4.441e-16, 0.0, -4.441e-16]\n",
"\n",
"Backward error is 4.441e-16\n",
"\n",
"Error x_random - x_computed is [0.0, 1.11e-16, 0.0]\n",
"\n",
"Absolute error |x_random-x_computed| is 1.11e-16\n"
]
}
],
"source": [
"println(\"A is\");print_matrix(A)\n",
"println(\"\\nbrandom is $b_random\")\n",
"(U,c_random) =row_reduce(A,b_random)\n",
"println(\"\\nU is\");print_matrix(U)\n",
"println(\"\\nResidualc_random - U*x_random = $(round.(c_random - U*x_random,4))\")\n",
"x_computed = backward_substitution(U,c_random)\n",
"println(\"\\nx_computed is $(round.(x_computed,12))\")\n",
"r = b_random - A*x_computed\n",
"println(\"\\nResidual b_random-A*x_computed is $(round.(r,4))\")\n",
"println(\"\\nBackward error is $(round(maximum(abs.(r)),4))\")\n",
"x_error =x_random - x_computed\n",
"println(\"\\nError x_random - x_computed is $(round.(x_error,4))\")\n",
"println(\"\\nAbsolute error |x_random-x_computed| is $(round(maximum(abs.(x_error)),4))\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What can go wrong? Some examples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} An obvious division by zero problem\n",
":label: example-obvious-division-by-zero\n",
"\n",
"Consider the system of two equations\n",
"\n",
"\\begin{align*}\n",
"x_2 &= 1\n",
"\\\\\n",
"x_1 + x_2 &= 2\n",
"\\end{align*}\n",
"\n",
"It is easy to see that this has the solution $x_1 = x_2 = 1$;\n",
"in fact it is already in \"reduced form\".\n",
"However when put into matrix form\n",
"\n",
"$$\n",
"\\left[\\begin{array}{rr} 0 & 1 \\\\ 1 & 1 \\end{array}\\right]\n",
"\\left[\\begin{array}{r} x_1 \\\\ x_2 \\end{array}\\right] = \\left[\\begin{array}{r} 1 \\\\ 2 \\end{array}\\right]\n",
"$$\n",
"\n",
"the above algorithm fails, because the fist *pivot element* $a_{11}$ is zero:\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A1 is)\n",
"[ 0.0 1.0 \n",
" 1.0 1.0 ]\n",
"b1 is [1.0, 1.0]\n"
]
}
],
"source": [
"A1 = [0.0 1.0 ; 1.0 1.0]\n",
"b1 = [1.0 ; 1.0]\n",
"\n",
"println(\"A1 is)\"); print_matrix(A1)\n",
"println(\"b1 is $(b1)\")"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U1 is\n",
"[ 0.0 1.0 \n",
" 0.0 -Inf ]\n",
"c1 is [1.0, -Inf]x1 is [NaN, NaN]"
]
}
],
"source": [
"(U1, c1) =row_reduce(A1, b1)\n",
"x1 = backward_substitution(U1, c1)\n",
"\n",
"println(\"U1 is\");print_matrix(U1)\n",
"print(\"c1 is $c1\")\n",
"print(\"x1 is $x1\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} IEEE fake numbers Inf and NaN\n",
":label: inf-nan\n",
"\n",
"- `Inf`, meaning \"infinity\", is a special value given as the result of calculations like division by zero.\n",
"Surprisingly, it can have a sign!\n",
"- `NaN`, meaning \"Not a Number\", is a special value given as the result of a calculation like `0/0`.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} A less obvious division by zero problem\n",
":label: example-less-obvious-division-by-zero\n",
"\n",
"Next consider this system\n",
"\n",
"$$\n",
"\\left[\\begin{array}{rrr} 1 & 1 & 1 \\\\ 1 & 1 & 2 \\\\ 1 & 2 & 2 \\end{array}\\right]\n",
"\\left[\\begin{array}{r} x_1 \\\\ x_2 \\\\ x_3 \\end{array}\\right] = \\left[\\begin{array}{r} 3 \\\\ 4 \\\\ 5 \\end{array}\\right]\n",
"$$\n",
"\n",
"The solution is $x_1 = x_2 = x_3 = 1$, and this time none of th diagonal elements is zero,\n",
"so it is not so obvious that a division by zero problem will occur, but:\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A2 is\n",
"[ 1.0 1.0 1.0 \n",
" 1.0 1.0 2.0 \n",
" 1.0 2.0 2.0 ]\n",
"b2 is [3.0, 4.0, 5.0]\n"
]
}
],
"source": [
"A2 = [1.0 1.0 1.0 ; 1.0 1.0 2.0 ; 1.0 2.0 2.]\n",
"b2 = [3.0 ; 4.0 ; 5.]\n",
"\n",
"println(\"A2 is\");print_matrix(A2)\n",
"println(\"b2 is $b2\")"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U2 is\n",
"[ 1.0 1.0 1.0 \n",
" 0.0 0.0 1.0 \n",
" 0.0 0.0 -Inf ]\n",
"c2 is [3.0, 1.0, -Inf]\n",
"x2 is [NaN, NaN, NaN]\n"
]
}
],
"source": [
"(U2, c2) =row_reduce(A2, b2)\n",
"x2 = backward_substitution(U2, c2)\n",
"\n",
"println(\"U2 is\");print_matrix(U2)\n",
"println(\"c2 is $c2\")\n",
"println(\"x2 is $x2\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What happens here is that the first stage subtracts the first row from each of the others ..."
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [],
"source": [
"A2[2,:] -= A2[1,:]\n",
"b2[2] -= b2[1]\n",
"A2[3,:] -= A2[1,:]\n",
"b2[3] -= b2[1];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"... and the new matrix has the same problem as above at the next stage:"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Now A2 is\n",
"[ 1.0 1.0 1.0 \n",
" 0.0 0.0 1.0 \n",
" 0.0 1.0 1.0 ]\n",
"and b2 is [3.0, 1.0, 2.0]\n"
]
}
],
"source": [
"println(\"Now A2 is\");print_matrix(A2)\n",
"println(\"and b2 is $(b2)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus, the second and third equations are\n",
"\n",
"$$\n",
"\\left[\\begin{array}{rr} 0 & 1 \\\\ 1 & 1 \\end{array}\\right]\n",
"\\left[\\begin{array}{r} x_2 \\\\ x_3 \\end{array}\\right] = \\left[\\begin{array}{r} 1 \\\\ 2 \\end{array}\\right]\n",
"$$\n",
"\n",
"with the same problem as in {prf:ref}`example-obvious-division-by-zero`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} Problems caused by inexact arithmetic: \"divison by almost zero\"\n",
":label: example-almost-division-by-zero\n",
"\n",
"The equations\n",
"\n",
"$$\n",
"\\left[\\begin{array}{rr} 1 & 10^{16} \\\\ 1 & 1 \\end{array}\\right]\n",
"\\left[\\begin{array}{r} x_1 \\\\ x_2 \\end{array}\\right] = \\left[\\begin{array}{r} 1+10^{16} \\\\ 2 \\end{array}\\right]\n",
"$$\n",
"\n",
"again have the solution $x_1 = x_2 = 1$, and the only division that happens in the above algorithm for row reduction is by that pivot element\n",
"$a_{11} = 1, \\neq 0$, so with exact arithmetic, all would be well. But:\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A3 is\n",
"[ 1.0 1.0e16 \n",
" 1.0 1.0 ]\n",
"b3 is [1.0e16, 2.0]\n"
]
}
],
"source": [
"A3 = [1.0 1e16 ; 1.0 1.0]\n",
"b3 = [1.0 + 1e16 ; 2.0]\n",
"\n",
"println(\"A3 is\");print_matrix(A3)\n",
"println(\"b3 is $b3\")"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U3 is\n",
"[ 1.0 1.0e16 \n",
" 0.0 -1.0e16 ]\n",
"c3 is [1.0e16, -9.999999999999998e15]\n",
"x3 is [2.0, 0.9999999999999998]\n"
]
}
],
"source": [
"(U3, c3) =row_reduce(A3, b3)\n",
"x3 = backward_substitution(U3, c3)\n",
"\n",
"println(\"U3 is\");print_matrix(U3)\n",
"println(\"c3 is $c3\")\n",
"println(\"x3 is $x3\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This gets $x_2 = 1$ fairly accurately, but $x_1$ is completely wrong!\n",
"\n",
"One hint is that $b_1$, which should be $1 + 10^{16} = 1000000000000001$, is instead just given as $10^{16}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On the other hand, all is well with less large values, like $10^{15}$:"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A3a is\n",
"[ 1.0 1.0e15 \n",
" 1.0 1.0 ]\n",
"b3a is [1.000000000000001e15, 2.0]\n"
]
}
],
"source": [
"A3a = [1.0 1e15 ; 1.0 1.0]\n",
"b3a = [1.0 + 1e15 ; 2.0]\n",
"\n",
"println(\"A3a is\");print_matrix(A3a)\n",
"println(\"b3a is $b3a\")"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U3a is\n",
"[ 1.0 1.0e15 \n",
" 0.0 -9.99999999999999e14 ]\n",
"c3a is [1.000000000000001e15, -9.99999999999999e14]\n",
"x3a is [1.0, 1.0]\n"
]
}
],
"source": [
"(U3a, c3a) =row_reduce(A3a, b3a)\n",
"x3a = backward_substitution(U3a, c3a)\n",
"\n",
"println(\"U3a is\");print_matrix(U3a)\n",
"println(\"c3a is $c3a\")\n",
"println(\"x3a is $x3a\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:example} Avoiding small denominators\n",
":label: example-avoiding-small-denominators\n",
"\n",
"The first equation in {prf:ref}`example-almost-division-by-zero` can be divided by $10^{16}$ to get an equivalent system with the same problem:\n",
"\n",
"$$\n",
"\\left[\\begin{array}{rr} 10^{-16} & 1 \\\\ 1 & 1 \\end{array}\\right]\n",
"\\left[\\begin{array}{r} x_1 \\\\ x_2 \\end{array}\\right] = \\left[\\begin{array}{r} 1+10^{-16} \\\\ 2 \\end{array}\\right]\n",
"$$\n",
"\n",
"Now the problem is more obvious: this system differs from the system in {prf:ref}`example-obvious-division-by-zero`\n",
"just by a tiny change of $10^{-16}$ in that pivot elements $a_{11}$, and the problem is *division by a value very close to zero*.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"A4 is\n",
"[ 1.0e-16 1.0 \n",
" 1.0 1.0 ]\n",
"b4 is [1.0, 2.0]\n"
]
}
],
"source": [
"A4 = [1e-16 1.0 ; 1.0 1.0]\n",
"b4 = [1.0 + 1e-16 ; 2.0]\n",
"\n",
"println(\"A4 is\");print_matrix(A4)\n",
"println(\"b4 is $b4\")"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"U4 is\n",
"[ 1.0e-16 1.0 \n",
" 0.0 -1.0e16 ]\n",
"c4 is [1.0, -9.999999999999998e15]\n",
"x4 is [2.220446049250313, 0.9999999999999998]\n"
]
}
],
"source": [
"(U4, c4) =row_reduce(A4, b4)\n",
"x4 = backward_substitution(U4, c4)\n",
"\n",
"println(\"U4 is\");print_matrix(U4)\n",
"println(\"c4 is $c4\")\n",
"println(\"x4 is $x4\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One might think that there is no such small denominator in {prf:ref}`example-almost-division-by-zero`,\n",
"but what counts for being \"small\" is magnitude relative to other values — 1 is very small compared to $10^{16}$.\n",
"\n",
"To understand these problems more (and how to avoid them) we will explore {ref}`section:machine-numbers-rounding-error-error-propagation` in the next section."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## When naive row reduction is safe: diagonal dominance\n",
"\n",
"There are several important cases when we can guarantee that these problem do not occur.\n",
"One obvious case is when the matrix $A$ is diagonal and non-singular (so with all non-zero elements);\n",
"then it is already row-reduced and with all denominators in backward substitution being non-zero.\n",
"\n",
"A useful measure of being \"close to diagonal\" is *diagonal dominance*:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:definition} Strict Diagonal Dominance\n",
":label: definition-strictly-diagonally-dominant\n",
"\n",
"A matrix $A$ is **row-wise strictly diagonally dominant**,\n",
"sometimes abbreviated as just **strictly diagonally dominant** or **SDD**,\n",
"if \n",
"\n",
"$$\\sum_{1 \\leq k \\leq n, k \\neq i}|a_{i,k}| < |a_{i,i}|$$\n",
"\n",
"Loosely, each main diagonal \"dominates\" in size over all other elements in its row.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:definition} Column-wise Strict Diagonal Dominance\n",
":label: definition-columnwise-strictly-diagonally-dominant\n",
"\n",
"If instead\n",
"\n",
"$$\\sum_{1 \\leq k \\leq n, k \\neq i}|a_{k,i}| < |a_{i,i}|$$\n",
"(so that each main diagonal element \"dominates its column\")\n",
"the matrix is called **column-wise strictly diagonally dominant**.\n",
"\n",
"Note that this is the same as saying that the transpose $A^T$ is SDD.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark}\n",
":label: non-strict-diagonal-dominance\n",
"\n",
"**Aside:** If only the corresponding non-strict inequality holds, the matrix is called *diagonally dominant*.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:theorem}\n",
":label: theorem-row-reduction-preserves-sdd\n",
"\n",
"For any strictly diagonally dominant matrix $A$, each of the intermediate matrices $A^{(k)}$ given by the naive Gaussan elimination algorithm is also strictly diagonally dominant, and so the final upper triangular matrix $U$ is.\n",
"In particular, all the diagonal elements $a_{i,i}^{(k)}$ and $u_{i,i}$ are non-zero, so no division by zero occurs in any of these algorithms, including the backward substitution solving for $x$ in $Ux = c$.\n",
"\n",
"The corresponding fact also true if the matrix is column-wise strictly diagonally dominant: that property is also preserved at each stage in naive row reduction.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus in each case the diagonal elements — the elements divided by in both row reduction and backward substitution — are in some sense safely away from zero.\n",
"We will have more to say about this in the sections on {ref}`section:linear-equations-pivoting`\n",
"and {ref}`section:linear-equations-lu-factorization`\n",
"\n",
"For a column-wise SDD matrix, more is true: at stage $k$, the diagonal dominance says that\n",
"the pivot element on the diagonal, $a_{k,k}^{(k-1)}$, is larger (in magnitude) than any of the elements $a_{i,k}^{(k-1)}$ below it, so the multipliers $l_{i,k}$ have\n",
"\n",
"$$|l_{i,k}| = |a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}| < 1.$$\n",
"\n",
"As we will see when we look at the effects of rounding error in {ref}`section:machine-numbers-rounding-error-error-propagation` and\n",
"{ref}`section:linear-equations-error-bounds`,\n",
"keeping intermediate values small is generally good for accuracy, so this is a nice feature."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{prf:remark} Positive definite matrices\n",
":label: remark-positive-definite-matrices-also-work\n",
"\n",
"Another class of matrices for which naive row-reduction works well is **positive definite matrices** which arise in any important situations; that property is in some sense more natural than diagonal dominance.\n",
"However that topic will be left for later.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{exercise-start}\n",
":label: linear-equations-row-reduction-exercise-1\n",
"```\n",
"Solve $A \\vec{x} = \\vec b$\n",
"with\n",
"\n",
"$$\n",
"A = \\left[ \\begin{array}{ccc} 4. & 2. & 1. \\\\ 9. & 3. & 1. \\\\ 25. & 5. & 1. \\end{array} \\right],\n",
"\\quad\n",
"\\vec b = \\left[ \\begin{array}{c} 0.693147 \\\\ 1.098612 \\\\ 1.609438 \\end{array} \\right],\n",
"$$\n",
"using *naive row-reduction ,* computing all intermediate values as decimal approximations rounded to four significant digits -- no fractions!\n",
"Call this approximate solution $\\vec{x}_1$.\n",
"\n",
"Then compute the residual $\\vec{r}_1 = \\vec{b} - A \\vec{x}_1$ and its norm $\\| \\vec{r}_1 \\|_\\infty$.\n",
"This should be done with high accuracy (not rounding to four decimal places) and could be done using a Julia notebook as a \"calculator\".\n",
"\n",
"Next, use Julia software to compute the *condition number* (see {prf:ref}`condition-number`)of $A$, $\\kappa(A) = \\kappa_\\infty(A)$,\n",
"(see section {ref}`section:linear-equations-error-bounds`),\n",
"and use this to get a bound on the absolute error in each of the above approximations.\n",
"\n",
"Finally some \"cheating\": use Julia software to compute the \"exact\" solution and use this to compute the actual absolute error;\n",
"compare this to the bound computed above.\n",
"```{exercise-end}\n",
"```"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.11.4",
"language": "julia",
"name": "julia-1.11"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}