From bc721d8bec252ae49ba31f24f94ab0bdd0581e93 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Thu, 17 Feb 2022 12:39:14 +0100 Subject: [PATCH] Fetch assignments week 3 --- Week 3/7 Linear Equation Systems.ipynb | 669 +++++++++++++++++++++++++ 1 file changed, 669 insertions(+) create mode 100644 Week 3/7 Linear Equation Systems.ipynb diff --git a/Week 3/7 Linear Equation Systems.ipynb b/Week 3/7 Linear Equation Systems.ipynb new file mode 100644 index 0000000..e58aa77 --- /dev/null +++ b/Week 3/7 Linear Equation Systems.ipynb @@ -0,0 +1,669 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4ec40081b048ce2f34f3f4fedbb0be10", + "grade": false, + "grade_id": "cell-98f724ece1aacb67", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "# CDS: Numerical Methods Assignments\n", + "\n", + "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n", + "\n", + "- Solutions must be submitted via the Jupyter Hub.\n", + "\n", + "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n", + "\n", + "## Submission\n", + "\n", + "1. Name all team members in the the cell below\n", + "2. make sure everything runs as expected\n", + "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n", + "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n", + "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n", + "6. submit your solutions **in time (before the deadline)**" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "team_members = \"\"" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "2900c663e4345c7f0707be39cc0cc7f4", + "grade": false, + "grade_id": "cell-b6b5c93e38567117", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "## Linear Equation Systems\n", + "\n", + "In the following you will implement the Gauss-Seidel (GS), Steepest Descent (SD) and the Conjugate Gradient (CG) algorithms to solve linear equation systems of the form \n", + "\n", + "$$A \\mathbf{x} = \\mathbf{b},$$ \n", + "\n", + "with $A$ being an $n \\times n$ matrix." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "55d777356d474cb7327bb087dfe5e644", + "grade": true, + "grade_id": "cell-bcf2dd8f9a194be8", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Import packages here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "9717dea99ca91d38963fc25b2ef95e03", + "grade": false, + "grade_id": "cell-595bb99f65f79ca7", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 1\n", + "First, you need to implement a Python function $\\text{diff(a,b)}$, which returns the difference $\\text{d}$ between two $n$-dimensional vectors $\\text{a}$ and $\\text{b}$ according to \n", + "\n", + "$$ d = || \\mathbf{a} - \\mathbf{b}||_\\infty = \\underset{i=1,2,\\dots,n}{\\operatorname{max}} |a_i - b_i|. $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "437477eff36ed8c083b4903a6921591c", + "grade": true, + "grade_id": "cell-bcd26d8a9f0e6447", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def diff(a, b):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "a10d859fa1b26c11ae00c7b846fdd1e8", + "grade": false, + "grade_id": "cell-19ae2e6fe5c9264c", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 2 \n", + "\n", + "The Gauss-Seidel iteration scheme to solve the linear equation system \n", + "\n", + "$$A \\mathbf{x} = \\mathbf{b}$$\n", + "\n", + "is defined by \n", + "\n", + "$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right].$$\n", + "\n", + "Note especially the difference in the sums: the first one involves $x_j^{(k)}$ and the second one $x_j^{(k-1)}$.\n", + "\n", + "\n", + "Give the outline of the derivation in LaTeX math notation in the markdown cell below. (Double click on \"YOUR ANSWER HERE\" to open the cell, and ctrl+enter to compile.) \n", + "\n", + "Hint: Similar to the Jacobi scheme, start by seperating the matrix $A$ into its diagonal ($D$), lower triangular ($L$) and upper triangular ($U$) forms, such that $A = D - L - U$." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "3c0c1bc83919a390739fa5802fc75743", + "grade": true, + "grade_id": "cell-3a8403fee85d9582", + "locked": false, + "points": 2, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "source": [ + "YOUR ANSWER HERE" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "dfbc2086b4294a7d942d35d9e1a22d5b", + "grade": false, + "grade_id": "cell-4ebd4c774509912d", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 3\n", + "\n", + "Implement the Gauss-Seidel iteration scheme derived above\n", + "\n", + "$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right],$$\n", + "\n", + "where $a_{ij}$ are the elements of the matrix $A$, and $x_i$ and $b_i$ the elements of vectors $\\mathbf{x}$ and $\\mathbf{b}$, respectively.\n", + "\n", + "Write a Python function $\\text{GS(A, b, eps)}$, where $\\text{A}$ represents the $n \\times n$ $A$ matrix, $\\text{b}$ represents the $n$-dimensional solution vector $\\mathbf{b}$, and $\\text{eps}$ is a scalar $\\varepsilon$ defining the accuracy up to which the iteration is performed. Your function should return both the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. \n", + "\n", + "Use an assertion to make sure the diagonal elements of $A$ are all non-zero. Initialize your iteration with $\\mathbf{x}^{(0)} = \\mathbf{0}$ (or with $\\mathbf{x}^{(1)} = D^{-1}\\mathbf{b}$, with $D$ the diagonal of $A$) and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "65b01c4e4c650dc6971b7d6dde62044b", + "grade": true, + "grade_id": "cell-ec2390c3cafa2f8e", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def GS(A, b, eps):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "badd9e6110882657c17c42e3d9ad8370", + "grade": false, + "grade_id": "cell-78d8cad45f5350b9", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 4\n", + "\n", + "Verify your implementation by comparing your approximate result to an exact solution. Use $\\text{numpy.linalg.solve()}$ to obtain the exact solution of the system\n", + "\n", + "$$\n", + "\\begin{align*}\n", + " \\begin{pmatrix}\n", + " 10 & -1 & 2 & 0 \\\\ \n", + " -1 & 11 &-1 & 3 \\\\\n", + " 2 & -1 & 10&-1 \\\\\n", + " 0 & 3 & -1& 8\n", + " \\end{pmatrix} \\mathbf{x}^*\n", + " =\n", + " \\begin{pmatrix}\n", + " 6 \\\\\n", + " 25 \\\\\n", + " -11\\\\\n", + " 15\n", + " \\end{pmatrix}\n", + "\\end{align*}\n", + "$$\n", + "\n", + "Then compare you approximate result $\\mathbf{\\tilde{x}}$ to the exact result $\\mathbf{x^*}$ by plotting $|| \\mathbf{x}^* - \\mathbf{\\tilde{x}}||_\\infty$ for different accuracies $\\varepsilon = 10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}$. \n", + "\n", + "Implement a unit test for your function using this system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "b54e0939053f9afee0f0ba46e490ba17", + "grade": true, + "grade_id": "cell-0585da22f1d304ab", + "locked": false, + "points": 4, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# Do your plotting here ...\n", + "\n", + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "d8e54a45261b013b9dd07a401d86401c", + "grade": true, + "grade_id": "cell-6524fb6d322a6ea0", + "locked": false, + "points": 0, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_GS():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_GS()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "f88ad91685d2bd73471acc5b600b5177", + "grade": false, + "grade_id": "cell-59514f17e337138f", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 5\n", + "\n", + "Next, implement the Steepest Descent algorithm in a similar Python function $\\text{SD(A, b, eps)}$, which calculates\n", + "\n", + "\\begin{align*}\n", + " \\mathbf{v}^{(k)} &= \\mathbf{b} - A \\mathbf{x}^{(k-1)} \\\\\n", + " t_k &= \\frac{ \\langle \\mathbf{v}^{(k)}, \\mathbf{v}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)}\\rangle } \\\\\n", + " \\mathbf{x}^{(k)} &= \\mathbf{x}^{(k-1)} + t_k \\mathbf{v}^{(k)} .\n", + "\\end{align*}\n", + "\n", + "Initialize your iteration again with $\\mathbf{x}^{(0)} = \\mathbf{0}$ and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n", + "Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "215218647aa6a6a314c093e82a4a24fb", + "grade": true, + "grade_id": "cell-a105ce7e1a34ce3c", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def SD(A, b, eps):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "f554d4f402c9e1a4cbff1f4e9018e9d0", + "grade": true, + "grade_id": "cell-62100bf55800145b", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_SD():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_SD()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "695745368fdf84d75691832839fb2834", + "grade": false, + "grade_id": "cell-3277aa3696a5c87f", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 6\n", + "\n", + "Finally, based on your steepest decent implementation from task 5, implement the Conjugate Gradient algorithm in a Python function $\\text{CG(A, b, eps)}$ in the following way: \n", + "\n", + "Initialize your procedure with:\n", + "\n", + "\\begin{align*}\n", + " \\mathbf{x}^{(0)} &= \\mathbf{0} \\\\\n", + " \\mathbf{r}^{(0)} &= \\mathbf{b} - A \\mathbf{x}^{(0)} \\\\\n", + " \\mathbf{v}^{(0)} &= \\mathbf{r}^{(0)}\n", + "\\end{align*}\n", + "\n", + "Then increase $k$ and repeat the following until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$.\n", + "\n", + "\\begin{align*}\n", + " t_k &= \\frac{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)} \\rangle } \\\\\n", + " \\mathbf{x}^{(k+1)} &= \\mathbf{x}^{(k)} + t_k \\mathbf{v}^{(k)} \\\\\n", + " \\mathbf{r}^{(k+1)} &= \\mathbf{r}^{(k)} - t_k A \\mathbf{v}^{(k)} \\\\\n", + " s_k &= \\frac{ \\langle \\mathbf{r}^{(k+1)}, \\mathbf{r}^{(k+1)} \\rangle }{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle } \\\\\n", + " \\mathbf{v}^{(k+1)} &= \\mathbf{r}^{(k+1)} + s_k \\mathbf{v}^{(k)}\n", + "\\end{align*}\n", + "\n", + "Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n", + "Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products.\n", + "\n", + "How do you expect the number of needed iteration steps to behave when changing the accuracy $\\epsilon$? What do you see?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "f709b7bfc8af2ba2d607f60ed03e1074", + "grade": true, + "grade_id": "cell-8ac3d07c9d237e47", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def CG(A, b, eps):\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "0d5b479bee6c0297974c003b4ebd99f3", + "grade": true, + "grade_id": "cell-4f5bc81f40ddb3fa", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "def test_CG():\n", + " # YOUR CODE HERE\n", + " raise NotImplementedError()\n", + " \n", + "test_CG()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "65dd782269d01314b02643d0610e3348", + "grade": false, + "grade_id": "cell-308595a088486388", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 7\n", + "\n", + "Apply all three methods to the following system\n", + "\n", + "\\begin{align*}\n", + "\\begin{pmatrix}\n", + "0.2& 0.1& 1.0& 1.0& 0.0 \\\\ \n", + "0.1& 4.0& -1.0& 1.0& -1.0 \\\\\n", + "1.0& -1.0& 60.0& 0.0& -2.0 \\\\\n", + "1.0& 1.0& 0.0& 8.0& 4.0 \\\\\n", + "0.0& -1.0& -2.0& 4.0& 700.0\n", + "\\end{pmatrix} \\mathbf{x}^*\n", + "=\n", + "\\begin{pmatrix}\n", + "1 \\\\\n", + "2 \\\\\n", + "3 \\\\\n", + "4 \\\\\n", + "5\n", + "\\end{pmatrix}.\n", + "\\end{align*}\n", + " \n", + "Plot the number of needed iterations for each method as a function of $\\varepsilon$, using $\\varepsilon = 10^{-1}, 10^{-2}, ..., 10^{-8}$.\n", + "\n", + "Explain the observed behavior with the help of the condition number (which you can calculate using $\\text{numpy.linalg.cond()}$). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "3d94dc67b052b82dae655b7ff562c7d2", + "grade": true, + "grade_id": "cell-490cd96dc58cbff8", + "locked": false, + "points": 4, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "deletable": false, + "editable": false, + "nbgrader": { + "cell_type": "markdown", + "checksum": "4d7421993a0fef53e4f0f83ffae975b1", + "grade": false, + "grade_id": "cell-c1b8b533da043a26", + "locked": true, + "schema_version": 3, + "solution": false, + "task": false + } + }, + "source": [ + "### Task 8\n", + "\n", + "Try to get a better convergence behavior by pre-conditioning your matrix $A$. Instead of $A$ use\n", + "\n", + "$$ \\tilde{A} = C A C,$$\n", + "\n", + "where $C = \\sqrt{D^{-1}}$. If you do so, you will need to replace $\\mathbf{b}$ by \n", + "\n", + "$$\\mathbf{\\tilde{b}} = C \\mathbf{b}$$\n", + "\n", + "and the vector $\\mathbf{\\tilde{x}}$ returned by your function will have to be transformed back via\n", + "\n", + "$$\\mathbf{x} = C \\mathbf{\\tilde{x}}.$$ \n", + "\n", + "What is the effect of $C$ on the condition number and why?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "deletable": false, + "nbgrader": { + "cell_type": "code", + "checksum": "213ac0f9b57d67bb90bbd75dd3808955", + "grade": true, + "grade_id": "cell-de3d1bd5ccfa8dfb", + "locked": false, + "points": 3, + "schema_version": 3, + "solution": true, + "task": false + } + }, + "outputs": [], + "source": [ + "# YOUR CODE HERE\n", + "raise NotImplementedError()" + ] + } + ], + "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.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}