\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: gaussian-elimination-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!), so there is no need to compute or store those zeros, 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: gaussian-elimination\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", "$end\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} Syntax for for loops and 0-based array indexing\n", ":label: remark-python-for-0-based\n", "\n", "Since array indices in Python (and in Java, C, C++, C#, Swift, etc.) start from zero, not from one, it will be convenient to express linear algebra algorithms in a form compatible with this.\n", "\n", "- Every index is one less than in the above!\n", "Thus in an array with $n$ elements, the index values $i$ are $0 \\leq i < n$, **excluding n**, which is the half-open interval of integers $[0, n)$.\n", "\n", "- In the indexing of an array, one can refer to the part the array with indices $a \\leq i < b$, **excluding b**, with the **slice** notation `a:b`.\n", "\n", "- Similarly, when specifiying the range of consecutive integers $i$, $a \\leq i < b$ in a `for` loop, one can use the expression `range(a,b)`.\n", "\n", "Also, when indices are processed in order (from low to high), these notes will abuse notation slightly, refering to the values as a set — specifically, a semi-open interval of integers.\n", "\n", "For example, the above loop\n", "\n", " for j from k+1 to n:\n", "\n", "first gets all indices lowered by one, to\n", "\n", " for j from k to n-1:\n", "\n", "and then this will sometimes be described in terms of the set of `j` values:\n", "\n", " for j in [k,n):\n", "\n", "which in Python becomes\n", "\n", " for j in range(k, n):\n", " ```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This new notation needs care initially, but helps with clarity in the long run.\n", "For one thing, it means that the indices of an $n$-element array, $[0,n-1)$, are described by `range(0,n)` and by `0:n`.\n", "In fact, the case of \"starting at the beginning\", with index zero, can be abbreviated: `range(n)` is the same as `range(0,n)`, and `:b` ia the same as `0:b`.\n", "\n", "Another advantage is that the index ranges `a:b` and `b:c` together cover the same indices as `a:c`, with no gap or duplication of `b`, and likewise `range(a,b)` and `range(b,c)` combine to cover `range(a,c)`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The naive row reduction algorithm, in Pythonic zero-based pseudo-code\n", "\n", "Here the above notational shift is made, along with eliminating the above-noted redundant formulas for values that are either zero or are unchanged from the previous step.\n", "It is also convenient for $k$ to be the index of the row being used to reduce subsequent rows, and so also the index of the column in which values below the main diagonal are being set to zero." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:algorithm}\n", ":label: gaussian-elimination-0-based\n", "\n", "for k in [0, n-1):\n", "

$\\quad$ for i in [k+1, n):\n", "

$\\quad\\quad l_{i,k} = a_{i,k}^{(k)}/a_{k,k}^{(k)}\\qquad$ **If** $a_{k,k}^{(k)} \\neq 0$!\n", "

$\\quad\\quad$ for j in [k+1, n):\n", "

$\\quad\\quad\\quad a_{i,j}^{(k+1)} = a_{i,j}^{(k)} - l_{i,k} a_{k,j}^{(k)}$\n", "

$\\quad\\quad$ end\n", "

$\\quad\\quad b_i^{(k+1)} = b_i^{(k)} - l_{i,k} b_{k}^{(k)}$\n", "

$\\quad$ end\n", "

end\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The naive row reduction algorithm, in Python\n", "\n", "Conversion to actual Python 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^{(0)}$, etc. and each vector $b^{(0)}$, $b^{(1)}$,\n", "we overwite each in the same array." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark}\n", ":label: remark-1-to-0-easy\n", "\n", "We will see that this simplicity in translation is quite common once algorithms have been expressed with zero-based indexing. The main ugliness is with loops that count backwards; see below.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " for k in range(n-1):\n", " for i in range(k+1, n):\n", " L[i,k] = A[i,k] / A[k,k]\n", " for j in range(k+1, n):\n", " A[i,j] -= L[i,k] * A[k,j]\n", " b[i] -= L[i,k] * b[k]" ] }, { "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": 32, "metadata": {}, "outputs": [], "source": [ "def rowReduce(A, b, demoMode=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", " U = A.copy()\n", " c = b.copy()\n", " n = len(b)\n", " # The function zeros_like() is used to create L with the same size and shape as A,\n", " # and with all its elements zero initially.\n", " L = np.zeros_like(A)\n", " for k in range(n-1):\n", " for i in range(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 range(k+1, n):\n", " U[i,j] -= L[i,k] * U[k,j]\n", " # Insert the below-diagonal zeros in column k;\n", " # this is not important for calculations, 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", " return (U, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** As usual, you could omit the above `def` and instead import this function with" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "from numericalMethods import rowReduce" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "though that actually gives the slightly revised version seen below." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Row reduction gives\n", "U =\n", "[[ 4. 2. 7. ]\n", " [ 0. 3.5 -11.25]\n", " [ 0. 0. -11. ]]\n", "c = [2. 1.5 5. ]\n" ] } ], "source": [ "(U, c) = rowReduce(A, b)\n", "print(f\"Row reduction gives\\nU =\\n{U}\")\n", "print(f\"c = {c}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take advantage of the fact that we have used `la.solve` to get a very accurate approximation of the solution x of $Ax=b$;\n", "this should also solve $Ux=c$, so check the backward error, a.k.a. the *residual*:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "The residual (backward error), b-Ax, is [0.0000000e+00 0.0000000e+00 8.8817842e-16], with maximum norm 8.881784197001252e-16.\n", "\n", "The residual for the row reduced form, c-Ux, is [ 0.00000000e+00 -2.22044605e-16 0.00000000e+00], with maximum norm 2.220446049250313e-16.\n" ] } ], "source": [ "r_Ab = b - A@x\n", "print(f\"\\nThe residual (backward error), b-Ax, is {r_Ab}, with maximum norm {max(abs(r_Ab))}.\")\n", "r_Uc = c - U@x\n", "print(f\"\\nThe residual for the row reduced form, c-Ux, is {r_Uc}, with maximum norm {max(abs(r_Uc))}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} Array slicing in Python\n", ":label: python-array-slicing\n", "\n", "Operations on a sequence of array indices, with \"slicing\": vectorization\n", "\n", "Python code can specify vector operations on a range of indices $[c,d)$, referred to withthe slice notaiton `c:d`.\n", "For example, the *slice* notation `A[c:d,j]` refers to the array containing the $d-c$ elements `A[i,j]` for $i$ in the semi-open interval $[c,d)$.\n", "\n", "Thus, each of the three arithmetic calculations above can be specified over a range of index values in a single command, eliminating all the inner-most `for` loops;\n", "this is somtimes called *vectorization*.\n", "Only `for` loops that contains other `for` loops remain.\n", "\n", "Apart from mathematical elegance, this usually allows far faster execution.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " for k in range(n-1):\n", " L[k+1:n,k] = U[k+1:n,k] / U[k,k] # compute all the L values for column k\n", " for i in range(k+1, n):\n", " U[i,k+1:n] -= L[i,k] * U[k,k+1:n] # Update row i\n", " c[k+1:n] -= L[k+1:n,k] * c[k] # update c values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I will break my usual guideline by redefining `rowReduce`, since this is just a different statement of the same algorithm\n", "(though with some added \"demoMode\" output)." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "def rowReduce(A, b, demomode=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", " U = A.copy()\n", " c = b.copy()\n", " n = len(b)\n", " # The function zeros_like() is used to create L with the same size and shape as A, and with all its elements zero initially.\n", " L = np.zeros_like(A)\n", " for k in range(n-1):\n", " if demomode: print(f\"Step {k=}\")\n", " # compute all the L values for column k:\n", " L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0\n", " if demomode:\n", " print(f\"The multipliers in column {k+1} are {L[k+1:,k]}\")\n", " for i in range(k+1, n):\n", " U[i,k+1:n] -= L[i,k] * U[k,k+1:n] # Update row i\n", " c[k+1:n] -= L[k+1:n,k] * c[k] # update c values\n", " if demomode:\n", " # Insert the below-diagonal zeros in column k;\n", " # this is not important for calculations, 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:, k] = 0.\n", " print(f\"The updated matrix is\\n{U}\")\n", " print(f\"The updated right-hand side is\\n{c}\")\n", " return (U, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} another way to select some rows or columns of a matrix\n", ":label: python-array-dicing\n", "\n", "As a variant on slicing, one can give a list of indices to select rows or columns a of a matrix; for example:\n", "\n", " A[[r1 r2 r3], :]\n", "\n", "gives a three row part of array `A` and\n", "\n", " A[2:, [c1 c2 c3 c4]]\n", "\n", "selects the indicated four columns — but only from row 2 onwards.\n", "\n", "This gives another way to describe the update of the lower-right block `U[k+1:n,k+1:n]` with a single matrix multiplication:\n", "it is the *outer product* of part of column `k` of `L` after row `k` by the part of row `k` of `U` after column `k`.\n", "\n", "To specify that the piecws of `L` nd `U` are identifies as a 1-column matrix and a 1-row matrix respectively, rather than as vectors,\n", "the above \"row/column list\" method must be used, with the list being just `[k]` in each case.\n", "```" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "def rowReduce(A, b, demomode=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", " U = A.copy()\n", " c = b.copy()\n", " n = len(b)\n", " # The function zeros_like() is used to create L with the same size and shape as A, and with all its elements zero initially.\n", " L = np.zeros_like(A)\n", " for k in range(n-1):\n", " if demomode: print(f\"Step {k=}\")\n", " # compute all the L values for column k:\n", " L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0\n", " if demomode:\n", " print(f\"The multipliers in column {k+1} are {L[k+1:,k]}\")\n", " U[k+1:n,k+1:n] -= L[k+1:n,[k]] @ U[[k],k+1:n] # The new \"outer product\" method.\n", " # Insert the below-diagonal zeros in column k;\n", " # this is not important for calculations, 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:n,k] = 0.0\n", " \n", " c[k+1:n] -= L[k+1:n,k] * c[k] # update c values\n", " if demomode:\n", " U[k+1:n, k] = 0. # insert zeros in column k of U:\n", " print(f\"The updated matrix is\\n{U}\")\n", " print(f\"The updated right-hand side is\\n{c}\")\n", " return (U, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Repeating the above testing:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(U, c) = rowReduce(A, b, demomode=True)\n", "print()\n", "print(f\"U =\\n{U}\")\n", "print(f\"c = {c}\")" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "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", "\n", "```{prf:algorithm}\n", ":label: backward-substitution-1\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": [ "## The backward substitution algorithm in zero-based pseudo-code\n", "\n", "Again, a zero-based version is more convenient for programming in Python (or Java, or C++):\n", "\n", "```{prf:algorithm}\n", ":label: backward-substitution-\n", "\n", "$x_{n-1} = c_{n-1}/u_{n-1,n-1}$\n", "

for i from n-2 down to 0\n", "

$\\displaystyle \\quad x_i = \\frac{c_i - \\sum_{j=i+1}^{n-1} u_{i,j} x_j}{u_{i,i}}$\n", "

end\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} Indexing from the end of an array and counting backwards\n", ":label: python-counting-backwards\n", "\n", "To express the above backwards counting in Python, we have to deal with the fact that `range(a,b)` counts upwards and excludes the \"end value\" `b`.\n", "The first part is easy: the extended form `range(a, b, step)` increments by `step` instead of by one, so that `range(a, b, 1)` is the same as `range(a,b)`, and `range(a, b, -1)` counts down: $a, a-1, \\dots, b+1$.\n", "\n", "But it still stops just before $b$, so getting the values from $n-1$ down to $0$ requires using $b= -1$, and so the slightly quirky expression `range(n-1, -1, -1)`.\n", "\n", "One more bit of Python: for an $n$-element single-index array `v`, the sum of its elements $\\sum_{i=0}^{n-1} v_i$ is given by `sum(v)`.\n", "Thus $\\sum_{i=a}^{b-1} v_i$, the sum over a subset of indices $[a,b)$, is given by `sum(v[a:b])`.\n", "\n", "And remember that multiplication of Numpy arrays with `*` is pointwise.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The backward substitution algorithm in Python\n", "\n", "With all the above Python details, the core code for backward substitution is:\n", "\n", " x[n-1] = c[n-1]/U[n-1,n-1]\n", " for i in range(n-2, -1, -1):\n", " x[i] = (c[i] - sum(U[i,i+1:] * x[i+1:])) / U[i,i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} \n", ":label: mathematically-correct-notation\n", "\n", "Note that the backward substitution algorithm and its Python coding have a nice mathematical advantage over the row reduction algorithm above: the precise mathematical statement of the algorithm does not need any intermediate quantities distinguished by superscripts ${}^{(k)}$, 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, 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, fixed meanings, corresponding to the natural mathematical description of the process:\n", "the method of **LU factorization** introduced in\n", "{ref}`section:linear-algebra-lu-factorization`.\n", "```" ] }, { "cell_type": "markdown", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "```{prf:remark} Another way to count backwards along an array\n", ":label: another-way-to-count-backwards\n", "\n", "On the other hand, there is an elegant way access array elements \"from the top down\".\n", "Firstly (or \"lastly\") `x[-1]` is the last element: the same as `x[n-1]` when `n = len(x)`, but without needing to know that length $n$.\n", "\n", "More generally, `x[-i]` is `x[n-i]`.\n", "\n", "Thus, one possibly more elegant way to describe backward substitution is to count with an increasing index, the \"distance from the bottom\":\n", "from `x[n-1]` which is `x[-1]` to `x[0]`, which is `x[-n]`.\n", "That is, index `-i` replaces index $n - i$: \n", "\n", " x[-1] = c[-1]/U[-1,-1]\n", " for i in range(2, n+1):\n", " x[-i] = (c[-i] - sum(U[-i,1-i:] * x[1-i:])) / U[-i,-i]\n", "\n", "There is still the quirk of having to \"overshoot\", referring to `n+1` in `range` to get to final index `-n`.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a final demonstration, we put this second version of the code into a complete working Python function and test it:" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "def backwardSubstitution(U, c, demomode=False):\n", " \"\"\"Solve U x = c for b.\"\"\"\n", " n = len(c)\n", " x = np.zeros(n)\n", " x[-1] = c[-1]/U[-1,-1]\n", " if demomode: print(f\"x_{n} = {x[-1]}\")\n", " for i in range(2, n+1):\n", " x[-i] = (c[-i] - sum(U[-i,1-i:] * x[1-i:])) / U[-i,-i]\n", " if demomode: print(f\"x_{n-i+1} = {x[-i]}\")\n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which as usual is also available via\n", "\n", " from numericalMethods import backwardSubstitution" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = [ 1.81168831 -1.03246753 -0.45454545]\n", "\n", "The residual b - Ax = [0.0000000e+00 0.0000000e+00 8.8817842e-16],\n", "with maximum norm 8.88e-16.\n" ] } ], "source": [ "x = backwardSubstitution(U, c)\n", "print(f\"x = {x}\")\n", "r = b - A@x\n", "print(f\"\\nThe residual b - Ax = {r},\")\n", "print(f\"with maximum norm {max(abs(r)):.3}.\")" ] }, { "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": 59, "metadata": {}, "outputs": [], "source": [ "def solveLinearSystem(A, b): return backwardSubstitution(*rowReduce(A, b));" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "```{prf:remark} On Python\n", ":label: python-splat-*\n", "\n", "The `*` here takes the value to its right (a single tuple with two elements `U` and `c`)\n", "and \"unpacks\" it to the two separate variables `U` and `c` needed as input to `backwardSubstitution`\n", "```" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1.81168831, -1.03246753, -0.45454545])" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solveLinearSystem(A, b)" ] }, { "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": [ "Prefered Python style is to have all `import` statements at the top,\n", "but since this is the first time we've heard of module `random`,\n", "I did not want it to be mentioned mysteriously above." ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "import random" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_random = [-0.27352765 0.31619792 -0.02981991]\n" ] } ], "source": [ "x_random = empty(len(b)) # An array the same length as b, with no values specified yet\n", "for i in range(len(x)):\n", " x_random[i] = random.uniform(-1, 1) # gives random real value, from uniform distribution in [-1, 1]\n", "print(f\"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": 67, "metadata": {}, "outputs": [], "source": [ "b_random = A @ x_random" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A =\n", "[[ 4. 2. 7.]\n", " [ 3. 5. -6.]\n", " [ 1. -3. 2.]]\n", "\n", "b_random = [-0.67045415 0.9393261 -1.28176123]\n", "\n", "U=\n", "[[ 4. 2. 7. ]\n", " [ 0. 3.5 -11.25]\n", " [ 0. 0. -11. ]]\n", "\n", "Residual c_random - U@x_random = [0. 0. 0.]\n", "\n", "x_computed = [-0.27352765 0.31619792 -0.02981991]\n", "\n", "Residual b_random - A@x_computed = [0. 0. 0.]\n", "\n", "Backward error |b_random - A@x_computed| = 0.0\n", "\n", "Error x_random - x_computed = [0. 0. 0.]\n", "\n", "Absolute error |x_random - x_computed| = 0.0\n" ] } ], "source": [ "print(f\"A =\\n{A}\")\n", "print(f\"\\nb_random = {b_random}\")\n", "(U, c_random) = rowReduce(A, b_random)\n", "print(f\"\\nU=\\n{U}\")\n", "print(f\"\\nResidual c_random - U@x_random = {c_random - U@x_random}\")\n", "x_computed = backwardSubstitution(U, c_random)\n", "print(f\"\\nx_computed = {x_computed}\")\n", "print(f\"\\nResidual b_random - A@x_computed = {b_random - A@x_computed}\")\n", "print(f\"\\nBackward error |b_random - A@x_computed| = {max(abs(b_random - A@x_computed))}\")\n", "print(f\"\\nError x_random - x_computed = {x_random - x_computed}\")\n", "print(f\"\\nAbsolute error |x_random - x_computed| = {max(abs(x_random - x_computed))}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What can go wrong? Three 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": 71, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U1 = \n", "[[ 0. 1.]\n", " [ 0. -inf]]\n", "c1 = [ 1. -inf]\n", "x1 = [nan nan]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/xv/1pdcq_w93rl5n_1hk34hgvz1qc6snk/T/ipykernel_79766/2412833099.py:15: RuntimeWarning: divide by zero encountered in divide\n", " L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0\n", "/var/folders/xv/1pdcq_w93rl5n_1hk34hgvz1qc6snk/T/ipykernel_79766/1659329062.py:5: RuntimeWarning: invalid value encountered in scalar divide\n", " x[-1] = c[-1]/U[-1,-1]\n" ] } ], "source": [ "A1 = array([[0., 1.], [1. , 1.]])\n", "b1 = array([1., 1.])\n", "(U1, c1) = rowReduce(A1, b1)\n", "print(f\"U1 = \\n{U1}\")\n", "print(f\"c1 = {c1}\")\n", "x1 = backwardSubstitution(U1, c1)\n", "print(f\"x1 = {x1}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{prf:remark} On Python \"Infinity\" and \"Not a Number\"\n", "\n", "- `inf`, meaning \"infinity\", is a special value given as the result of operations like division by zero.\n", "Surprisingly, it can have a sign!\n", "(This is available in Python from package Numpy as `numpy.inf`)\n", "- `nan`, meaning \"not a number\", is a special value given as the result of calculations like `0/0`.\n", "(This is available in Python from package Numpy as `numpy.nan`)\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 the 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": 74, "metadata": {}, "outputs": [], "source": [ "A2 = array([[1., 1., 1.], [1., 1., 2.],[1., 2., 2.]])\n", "b2 = array([3., 4., 5.])" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U2 = \n", "[[ 1. 1. 1.]\n", " [ 0. 0. 1.]\n", " [ 0. 0. -inf]]\n", "c2 = [ 3. 1. -inf]\n", "x2 = [nan nan nan]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/var/folders/xv/1pdcq_w93rl5n_1hk34hgvz1qc6snk/T/ipykernel_79766/2412833099.py:15: RuntimeWarning: divide by zero encountered in divide\n", " L[k+1:,k] = U[k+1:n,k] / U[k,k] # Beware the case where U[k,k] is 0\n", "/var/folders/xv/1pdcq_w93rl5n_1hk34hgvz1qc6snk/T/ipykernel_79766/1659329062.py:5: RuntimeWarning: invalid value encountered in scalar divide\n", " x[-1] = c[-1]/U[-1,-1]\n" ] } ], "source": [ "(U2, c2) = rowReduce(A2, b2)\n", "print(f\"U2 = \\n{U2}\")\n", "print(f\"c2 = {c2}\")\n", "x2 = backwardSubstitution(U2, c2)\n", "print(f\"x2 = {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": 77, "metadata": {}, "outputs": [], "source": [ "A2[1,:] -= A2[0,:]\n", "b2[1] -= b2[0]\n", "A2[2,:] -= A2[0,:]\n", "b2[2] -= b2[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... and the new matrix has the same problem as above at the next stage:" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Now A2 is \n", "[[1. 1. 1.]\n", " [0. 0. 1.]\n", " [0. 1. 1.]]\n", "and b2 is [3. 1. 2.]\n" ] } ], "source": [ "print(f\"Now A2 is \\n{A2}\")\n", "print(f\"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": 82, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A3 = \n", "[[1.e+00 1.e+16]\n", " [1.e+00 1.e+00]]\n", "b3 = [1.e+16 2.e+00]\n" ] } ], "source": [ "A3 = array([[1., 1e16], [1. , 1.]])\n", "b3 = array([1. + 1e16, 2.])\n", "print(f\"A3 = \\n{A3}\")\n", "print(f\"b3 = {b3}\")" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U3 = \n", "[[ 1.e+00 1.e+16]\n", " [ 0.e+00 -1.e+16]]\n", "c3 = [ 1.e+16 -1.e+16]\n", "x3 = [2. 1.]\n" ] } ], "source": [ "(U3, c3) = rowReduce(A3, b3)\n", "print(f\"U3 = \\n{U3}\")\n", "print(f\"c3 = {c3}\")\n", "x3 = backwardSubstitution(U3, c3)\n", "print(f\"x3 = {x3}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gets $x_2 = 1$ correct, 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": 86, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A3a = \n", "[[1.e+00 1.e+15]\n", " [1.e+00 1.e+00]]\n", "b3a = [1.e+15 2.e+00]\n" ] } ], "source": [ "A3a = array([[1., 1e15], [1. , 1.]])\n", "b3a = array([1. + 1e15, 2.])\n", "print(f\"A3a = \\n{A3a}\")\n", "print(f\"b3a = {b3a}\")" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U3a = \n", "[[ 1.e+00 1.e+15]\n", " [ 0.e+00 -1.e+15]]\n", "c3a = [ 1.e+15 -1.e+15]\n", "x3a = [1. 1.]\n" ] } ], "source": [ "(U3a, c3a) = rowReduce(A3a, b3a)\n", "print(f\"U3a = \\n{U3a}\")\n", "print(f\"c3a = {c3a}\")\n", "x3a = backwardSubstitution(U3a, c3a)\n", "print(f\"x3a = {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:\n", "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": 89, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A4 = \n", "[[1.e-16 1.e+00]\n", " [1.e+00 1.e+00]]\n", "b4 = [1. 2.]\n" ] } ], "source": [ "A4 = array([[1e-16, 1.], [1. , 1.]])\n", "b4 = array([1. + 1e-16, 2.])\n", "print(f\"A4 = \\n{A4}\")\n", "print(f\"b4 = {b4}\")" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "U4 = \n", "[[ 1.e-16 1.e+00]\n", " [ 0.e+00 -1.e+16]]\n", "c4 = [ 1.e+00 -1.e+16]\n", "x4 = [2.22044605 1. ]\n" ] } ], "source": [ "(U4, c4) = rowReduce(A4, b4)\n", "print(f\"U4 = \\n{U4}\")\n", "print(f\"c4 = {c4}\")\n", "x4 = backwardSubstitution(U4, c4)\n", "print(f\"x4 = {x4}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One might think that there is no such small denominator in\n", "{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\n", "[Machine Numbers, Rounding Error and Error Propagation](section:machine-numbers-rounding-error-error-propagation)\n", "in the next section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## When naive Guassian elimination 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", "$$\n", "\\sum_{1 \\leq k \\leq n, k \\neq i}|a_{i,k}| < |a_{i,i}|\n", "$$\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", "$$\n", "\\sum_{1 \\leq k \\leq n, k \\neq i}|a_{k,i}| < |a_{i,i}|\n", "$$\n", "\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", "````\n", "\n", "**Aside:** If only the corresponding non-strict inequality holds, the matrix is called *diagonally dominant*." ] }, { "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 Guassian elimination.\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\n", "[pivoting](section:linear-algebra-pivoting)\n", "and\n", "[LU factorization](section:linear-algebra-lu-factorization)\n", "\n", "For a column-wise SDD matrix, more is true: at stage $k$, the diagonal dominance says that\n", "the pivot elemet 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", "$$\n", "|l_{i,k}| = |a_{i,k}^{(k-1)}/a_{k,k}^{(k-1)}| < 1.\n", "$$\n", "\n", "As we will see when we look at the effects of rounding error in the sections on\n", "[Machine Numbers, Rounding Error and Error Propagation](section:machine-numbers-rounding-error-error-propagation)\n", "and\n", "[Error bounds for linear algebra](section:linear-algebra-errors)\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-algebra-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 Gaussian elimination,* 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 Python notebook as a \"calculator\".\n", "\n", "Next, use Python software to compute the condition number of $A$, $K(A) = K_\\infty(A)$,\n", "and use this to get a bound on the absolute error in each of the above approximations.\n", "\n", "Finally some \"cheating\": use Python 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": "Python 3 (ipykernel)", "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.12.4" } }, "nbformat": 4, "nbformat_minor": 4 }