Files
cds-numerical-methods/Week 3/7 Linear Equation Systems.ipynb

670 lines
19 KiB
Plaintext

{
"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
}