1088 lines
83 KiB
Plaintext
1088 lines
83 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 = \"Koen Vendrig, Kees van Kempen\""
|
|
]
|
|
},
|
|
{
|
|
"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": 1,
|
|
"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 numpy as np\n",
|
|
"import numpy.linalg as linalg\n",
|
|
"from matplotlib import pyplot as plt"
|
|
]
|
|
},
|
|
{
|
|
"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": 2,
|
|
"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",
|
|
" return np.max(np.abs(a - b))"
|
|
]
|
|
},
|
|
{
|
|
"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": [
|
|
"---\n",
|
|
"\n",
|
|
"We start from our linear equations:\n",
|
|
"\n",
|
|
"$$Ax = b$$\n",
|
|
"\n",
|
|
"We separate A into different components (diagonal, strictly lower triangular and strictly upper triangular):\n",
|
|
"\n",
|
|
"$$A = D - L - U$$\n",
|
|
"\n",
|
|
"\n",
|
|
"We write $D - L$ as $L'$ to get:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$(L' - U)x = b$$\n",
|
|
"\n",
|
|
"\n",
|
|
"We take the iterative process of the Gauss-Seidel method to write:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$\n",
|
|
"L'x^k = b + Ux^{k-1}\\\\\n",
|
|
"x^k = L'^{-1}(b + Ux^{k-1})\\\\\n",
|
|
"$$\n",
|
|
"If we write every component of the matrix $A$ as $a_{ij}$, we can use forward substitution to rewrite our previous equation to:\n",
|
|
"\n",
|
|
"\n",
|
|
"$$x^k _i = \\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].$$"
|
|
]
|
|
},
|
|
{
|
|
"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": 3,
|
|
"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, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Gauss-Seidel algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" # First we decompose A = D - L - U.\n",
|
|
" D = np.diag(np.diag(A))\n",
|
|
" U = -np.triu(A) + D\n",
|
|
" L = -np.tril(A) + D\n",
|
|
" \n",
|
|
" # We need non-zero diagonals elements.\n",
|
|
" assert np.all(np.diag(D) != 0)\n",
|
|
" \n",
|
|
" x_prev = np.zeros(n)\n",
|
|
" x_cur = np.dot(linalg.inv(D), b)\n",
|
|
" \n",
|
|
" k = 1\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max:\n",
|
|
" k += 1\n",
|
|
" # We will have to copy, as the array elements will point to the same\n",
|
|
" # memory otherwise, and changes to one array will change the other aswell.\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" for i in range(n):\n",
|
|
" x_cur[i] = 1/A[i, i]*(-np.dot(A[i, :i], x_cur[:i]) - np.dot(A[i, i + 1:], x_prev[i + 1:]) + b[i])\n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"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": 4,
|
|
"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": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQxUlEQVR4nO3df2id133H8c9nipOoGZVL445YiWcXBdGwEAzCgXgEN1mmFKLEhKxLukFLTEwysrF/RG3S9b/WYWJjLQlkLg1uyvLDpEbYiVcxSF2HEDYr0aiTGXUmECI5zMkWqU25UNv97g9dx7IiyffqPI+ee+95v0D4Pue599xvOFw+Oec897mOCAEAkOL3qi4AAND+CBMAQDLCBACQjDABACQjTAAAyQgTAECyy6ouoApXX311bNy4seoyAKCtvPHGGx9GxLrFzmUZJhs3btT4+HjVZQBAW7H97lLnslrmsj1ke+/s7GzVpQBAR8kqTCLiUETs7OnpqboUAOgoWYUJAKAchAkAIFlWG/C2hyQN9fX1VV0KAKyq0YlpjYxN6tRMTevXdmt4sF/bN/cW1n9WMxP2TADkaHRiWrsPHNf0TE0haXqmpt0Hjmt0Yrqw98gqTAAgRyNjk6qdOXdRW+3MOY2MTRb2HoQJAHS4UzO1ptpXgjABgA63fm13U+0rkVWY8KVFADkaHuxX95qui9q613RpeLC/sPfIKkzYgAeQo+2be7Xn3hvVu7ZbltS7tlt77r2x0Ku5sro0GABytX1zb6HhsVBWMxMAQDkIEwBAMsIEAJCMMAEAJMsqTLg0GADKkVWYcGkwAJQjqzABAJSDMAEAJCNMAADJCBMAQDLCBACQjDABACTLKkz4ngkAlCOrMOF7JgBQjqzCBABQDsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAECyrMKE26kAQDmyChNupwIA5cgqTAAA5SBMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJCBMAQLKOCBPbX7L9lO0XbT9SdT0AkJvKw8T207ZP235rQfudtidtn7S9a7k+IuJERDws6auStpZZLwDg0yoPE0n7JN05v8F2l6QnJX1F0g2SHrB9g+0bbb+04O8L9dfcLellSYdXt3wAwGVVFxARR21vXNC8RdLJiHhHkmw/L+meiNgj6a4l+jko6aDtlyU9u/C87Z2SdkrShg0bivsPAABUHyZL6JX03rzjKUk3L/Vk29sk3SvpCi0xM4mIvZL2StLAwEAUVCcAQK0bJk2JiCOSjlRcBgBkqxX2TBYzLem6ecfX1tuS2B6yvXd2dja1KwDAPK0aJsckXW97k+3LJd0v6WBqpxFxKCJ29vT0JBcIALig8jCx/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE21XWCQBYWuV7JhHxwBLth8VlvgDQFiqfmawm9kwAoBxZhQl7JgBQjqzCBABQjqzChGUuAChHVmHCMhdQvtGJaW19/BVt2vWytj7+ikYnkr8ihjZQ+dVcADrH6MS0dh84rtqZc5Kk6Zmadh84Lknavrm3ytJQsqxmJgDKNTI2+UmQnFc7c04jY5MVVYTVQpgAKMypmVpT7egcWYUJG/BAudav7W6qHZ0jqzBhAx4o1/Bgv7rXdF3U1r2mS8OD/RVVhNXCBjyAwpzfZB8Zm9SpmZrWr+3W8GA/m+8ZIEwAFGr75l7CI0NZLXMBAMqRVZiwAQ8A5cgqTNiAB4ByZBUmAIByECYAgGSECQAgGWECAEiWVZhwNRcAlCOrMOFqLgAoR6FhYvuzRfYHAGgPhd1OxfYjks7ZvjUi/rKofgEAra/Ie3Odqf97tsA+AQBtoMhlrvcl9Uo6XWCfAIA2UGSY3CzpNUmbCuwTANAGCguTiPi2pP+T9FBRfQIA2kOhV3NFxJsRMVNkn0XieyYAUI6mN+Btb2jwqTMR8atm+y9TRBySdGhgYIDZEwAUaCVXc/2ogeeEpH2SnllB/wCANtN0mETEl8soBADQvrJa5gIAlINlLgBAsqbCxHYXy1wAgIWWvTTYdq/t79r+TL3p720/UT/3vdKrAwC0hUt9z+RGSTdJGqgffyTp3frjX9s+ZPsqSbI9aPu1csoEALSyZZe5IuKntr8eEUfrTVsk/Vv93Ldsf03SEdu/lfSxpF2lVgsAaElN7ZlExN22PydJtm/X3K1TfiPpGkkPRsRk8SUCAFpdI7dTeXP+QUR8VH/4mKS/i4htku6T9ILt24otDwDQDi45M4mIkSXab5v3+Ljtr0j6iaRbiiuvWLaHJA319fVVXQoAdJQi7xr8vqTbi+qvDPwGPACUo+i7BteK7A8A0B6Sw6S+dAQAyFgRM5PvFNAHAKCNFREmLqAPAEAbKyJMooA+AABtrNANeABAnggTAECyIsLkfwroAwDQxpLDJCLuKKIQAED7YpkLAJCMMAEAJCNMAADJmgqTeb+q+PvllAMAaEfNzkw+Z/tRSX9cRjEAgPbUbJjcLukbkr5o+wvFl7Nytq+yPW77rqprAYDcNBsm/yHpQUnvRsTpIgqw/bTt07bfWtB+p+1J2ydtN/Lb8t+UtL+ImgAAzWn2N+BP1B/+osAa9kl6QtIz5xtsd0l6UtIdkqYkHbN9UFKXpD0LXv+gpJsk/ZekKwusCwDQoKbCpAwRcdT2xgXNWySdjIh3JMn285LuiYg9kj61jGV7m6SrJN0gqWb7cET8rsy6AQAXNLzMZbvX9ndtf6Z+/D3bZd1+vlfSe/OOp+pti4qIxyLibyU9K+kHiwWJ7Z31PZXxDz74oOh6ASBrzeyZ3Ki55aSB+vGvJR2cd7nwoO3XCq6vKRGxLyJeWuLc3ogYiIiBdevWrXZpANDRGl7mioif2v56RBytH3/L9tckHbH9W0kfS2pko7wR05Kum3d8bb0NANCCVvwNeNu3S3pI0m8kXS3pbyLi1YLqOibpetubbF8u6X5JB1M7tT1ke+/s7GxygQCAC5oNkzfnPX5M0rcjYpuk+yS9YPu2Zguw/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE2832vVBEHIqInT09PaldAQDmcUQxv7pr+xpJP4mIWwrpsEQDAwMxPj5edRkA0FZsvxERA4uda/rSYNsbljm9Y975mYj4VbP9l8n2kKShvr6+qksBgI7S9MzE9s+WOR2SXP93X0Q8s8xzK8PMBACaV+jMJCK+nF4SAKCTFL3MNV/LLXMBAMqxktup/KiB54Tm7rnVUstc7JkAQDkKu5qrnbBn0jlGJ6Y1MjapUzM1rV/breHBfm3fvOSddwAkKHTPBGgVoxPT2n3guGpnzkmSpmdq2n3guCQRKMAq4zfg0bZGxiY/CZLzamfOaWRssqKKgHxlFSbcTqWznJqpNdUOoDxZhQm3U+ks69d2N9UOoDxZhQk6y/Bgv7rXdF3U1r2mS8OD/RVVBOSLDXi0rfOb7FzNBVSPMEFb2765l/AAWkBWy1xswANAObIKEzbgAaAcWYUJAKAchAkAIBlhAgBIRpgAAJIRJgCAZFmFCZcGA0A5sgoTLg0GgHJkFSYAgHIQJgCAZIQJACAZYQIASEaYAACSESYAgGRZhQnfMwGAcmQVJnzPBADKkVWYAADKQZgAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEiWVZhwOxUAKEdWYcLtVACgHFmFCQCgHIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASNYRYWJ7m+1XbT9le1vV9QBAbioPE9tP2z5t+60F7XfanrR90vauS3QTkj6WdKWkqbJqBQAs7rKqC5C0T9ITkp4532C7S9KTku7QXDgcs31QUpekPQte/6CkVyPi57b/QNI/SvqLVagbAFBXeZhExFHbGxc0b5F0MiLekSTbz0u6JyL2SLprme4+knTFYids75S0U5I2bNiQWjYAYJ7Kl7mW0CvpvXnHU/W2Rdm+1/Y/S/qx5mY5nxIReyNiICIG1q1bV2ixAJC7ymcmRYiIA5IOVF0HAOSqVWcm05Kum3d8bb0tie0h23tnZ2dTuwIAzNOqYXJM0vW2N9m+XNL9kg6mdhoRhyJiZ09PT3KBAIALKg8T289Jel1Sv+0p2zsi4qykRyWNSTohaX9EvF1lnQCApVW+ZxIRDyzRfljS4SLfy/aQpKG+vr4iuwWA7FU+M1lNLHMBQDmyChMAQDkIEwBAsqzChEuDAaAcWYUJeyYAUI6swgQAUA7CBACQjDABACTLKkzYgAeAcmQVJqkb8KMT09r6+CvatOtlbX38FY1OJN97EgA6QuW3U2kXoxPT2n3guGpnzkmSpmdq2n3guCRp++Ylf2oFALKQ1cwkxcjY5CdBcl7tzDmNjE1WVBEAtA7CpEGnZmpNtQNATrIKk5QN+PVru5tqB4CcZBUmKRvww4P96l7TdVFb95ouDQ/2F1UeALQtNuAbdH6TfWRsUqdmalq/tlvDg/1svgOACJOmbN/cS3gAwCKyWuYCAJSDMAEAJMsqTLidCgCUI6sw4fdMAKAcWYUJAKAcjoiqa1h1tj+Q9O4ip3okLVwDW6ztakkfllDapSxWy2r108hrLvWc5c4vda7Vx0QqZlzKGpNGnlfWuLT7mKy0n07+rPxhRKxb9ExE8Ff/k7S3wbbxVqlvtfpp5DWXes5y55c61+pjUtS4lDUmVY5Lu49JmePSiZ8VlrkudqjBtqoUVctK+mnkNZd6znLnlzrX6mMiFVNPWWPSyPM6cVz4rDReSyGyXOZKZXs8IgaqrgMXMCathzFpTWWNCzOTldlbdQH4FMak9TAmramUcWFmAgBIxswEAJCMMAEAJCNMAADJCJOC2b7K9rjtu6quBXNsf8n2U7ZftP1I1fVAsr3d9g9sv2D7T6uuB3Nsf9H2D22/2OxrCZM620/bPm37rQXtd9qetH3S9q4GuvqmpP3lVJmfIsYlIk5ExMOSvippa5n15qCgMRmNiIckPSzpz8usNxcFjcs7EbFjRe/P1VxzbN8q6WNJz0TEH9XbuiT9UtIdkqYkHZP0gKQuSXsWdPGgpJskfV7SlZI+jIiXVqf6zlXEuETEadt3S3pE0o8j4tnVqr8TFTUm9df9g6R/iYg3V6n8jlXwuLwYEfc18/780mJdRBy1vXFB8xZJJyPiHUmy/bykeyJij6RPLWPZ3ibpKkk3SKrZPhwRvyuz7k5XxLjU+zko6aDtlyURJgkK+qxY0uOS/pUgKUZRn5WVIkyW1yvpvXnHU5JuXurJEfGYJNn+huZmJgRJOZoal3rI3yvpCkmHyywsY02NiaS/lvQnknps90XEU2UWl7FmPyufl/QdSZtt766HTkMIkxJExL6qa8AFEXFE0pGKy8A8EfF9Sd+vug5cLCL+V3P7WE1jA35505Kum3d8bb0N1WJcWg9j0ppWbVwIk+Udk3S97U22L5d0v6SDFdcExqUVMSatadXGhTCps/2cpNcl9duesr0jIs5KelTSmKQTkvZHxNtV1pkbxqX1MCatqepx4dJgAEAyZiYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJkCLqP+I0X/W//7dNp9PtA2+AQ+0CNv/LenWiHi/6lqAZvF/PkDrOCzpF7b/qepCgGbxeyZAC7B9iyRLuqZ+cz6grTAzAVrDn0n6ZUSc9ZzPVl0Q0Az2TIAWYHuLpB9KCkk1SX8VEW9UWxXQOMIEAJCMZS4AQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMn+H8is0a3Xr+ppAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"A = np.array([[ 10, - 1, 2, 0],\n",
|
|
" [- 1, 11, - 1, 3],\n",
|
|
" [ 2, - 1, 10, - 1],\n",
|
|
" [ 0, 3, - 1, 8]])\n",
|
|
"b = np.array( [ 6, 25, -11, 15] )\n",
|
|
"x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
"eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n",
|
|
"diff_list = []\n",
|
|
"for eps in eps_list:\n",
|
|
" x, k = GS(A, b, eps)\n",
|
|
" diff_list.append(diff(x_exact, x))\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"\n",
|
|
"ax.scatter(eps_list, diff_list)\n",
|
|
"ax.set_xscale(\"log\")\n",
|
|
"ax.set_yscale(\"log\")\n",
|
|
"ax.set_xlabel(\"$\\epsilon$\")\n",
|
|
"ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n",
|
|
"\n",
|
|
"fig.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# As the three algorithm functions will have the same signature,\n",
|
|
"# it makes sense to only write the test function once.\n",
|
|
"\n",
|
|
"def test_alg(alg, alg_name):\n",
|
|
" \"\"\"\n",
|
|
" Check that function alg returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" A = np.array([[ 10, - 1, 2, 0],\n",
|
|
" [- 1, 11, - 1, 3],\n",
|
|
" [ 2, - 1, 10, - 1],\n",
|
|
" [ 0, 3, - 1, 8]])\n",
|
|
" b = np.array( [ 6, 25, -11, 15] )\n",
|
|
" x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
" print(\"Starting with A =\")\n",
|
|
" print(A)\n",
|
|
" print(\"and b =\", b)\n",
|
|
" print(\"We apply the {} algorithm to solve Ax = b.\".format(alg_name))\n",
|
|
" print()\n",
|
|
"\n",
|
|
" eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations:\".format(eps, k))\n",
|
|
" print(\"x =\\t\\t\\t\", x)\n",
|
|
" print(\"Ax =\\t\\t\\t\", np.dot(A, x))\n",
|
|
" print(\"diff(Ax, b) =\\t\\t\", diff(A @ x, b))\n",
|
|
" print(\"diff(x, x_exact) =\\t\", diff(x, x_exact))\n",
|
|
" print()\n",
|
|
" \n",
|
|
" assert diff(x, x_exact) < eps\n",
|
|
" "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"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": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Gauss-Seidel algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 0.99463393 1.99776509 -0.99803257 1.00108402]\n",
|
|
"Ax =\t\t\t [ 5.95250909 24.98206671 -10.98990699 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.04749090616931895\n",
|
|
"diff(x, x_exact) =\t 0.005366066491359844\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 0.99938302 1.99982713 -0.99978549 1.00009164]\n",
|
|
"Ax =\t\t\t [ 5.99443213 24.99877578 -10.99900762 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.005567865937722516\n",
|
|
"diff(x, x_exact) =\t 0.000616975874427883\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 6\t iterations:\n",
|
|
"x =\t\t\t [ 0.99993981 1.99998904 -0.99997989 1.00000662]\n",
|
|
"Ax =\t\t\t [ 5.99944928 24.99993935 -10.99991498 15. ]\n",
|
|
"diff(Ax, b) =\t\t 0.000550717702960668\n",
|
|
"diff(x, x_exact) =\t 6.018928065554263e-05\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 7\t iterations:\n",
|
|
"x =\t\t\t [ 0.99999488 1.99999956 -0.99999836 1.00000037]\n",
|
|
"Ax =\t\t\t [ 5.99995255 24.99999971 -10.99999375 15. ]\n",
|
|
"diff(Ax, b) =\t\t 4.744782363452771e-05\n",
|
|
"diff(x, x_exact) =\t 5.11751035947583e-06\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_GS():\n",
|
|
" \"\"\"\n",
|
|
" Check that GS returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" return test_alg(GS, \"Gauss-Seidel\")\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": 7,
|
|
"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, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Steepest Descent algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" \n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" x_cur = np.zeros(n)\n",
|
|
" x_prev = np.zeros(n)\n",
|
|
" \n",
|
|
" k = 0\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
|
|
" k += 1\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" \n",
|
|
" v = b - A @ x_prev\n",
|
|
" t = np.dot(v, v)/np.dot(v, A @ v)\n",
|
|
" x_cur = x_prev.copy() + t*v\n",
|
|
" \n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"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": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Steepest Descent algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 0.99748613 1.98300329 -0.98904751 1.01283183]\n",
|
|
"Ax =\t\t\t [ 6.01376302 24.84309304 -10.89133793 15.04071202]\n",
|
|
"diff(Ax, b) =\t\t 0.15690696195356324\n",
|
|
"diff(x, x_exact) =\t 0.016996711694861055\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 6\t iterations:\n",
|
|
"x =\t\t\t [ 0.99983175 1.99716093 -0.99850509 1.00217552]\n",
|
|
"Ax =\t\t\t [ 6.00414638 24.97397012 -10.98472385 15.00739201]\n",
|
|
"diff(Ax, b) =\t\t 0.02602988052923294\n",
|
|
"diff(x, x_exact) =\t 0.002839069878101119\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 9\t iterations:\n",
|
|
"x =\t\t\t [ 0.99991029 1.99994247 -0.99999645 1.00023877]\n",
|
|
"Ax =\t\t\t [ 5.99916754 25.00016961 -11.00032515 15.00173404]\n",
|
|
"diff(Ax, b) =\t\t 0.0017340408511579142\n",
|
|
"diff(x, x_exact) =\t 0.00023877424067331177\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 11\t iterations:\n",
|
|
"x =\t\t\t [ 0.99998551 1.99999065 -0.99999949 1.00003874]\n",
|
|
"Ax =\t\t\t [ 5.9998655 25.00002733 -11.00005329 15.00028137]\n",
|
|
"diff(Ax, b) =\t\t 0.0002813662634846281\n",
|
|
"diff(x, x_exact) =\t 3.874139053650083e-05\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_SD():\n",
|
|
" \"\"\"\n",
|
|
" Check that SD returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" return test_alg(SD, \"Steepest Descent\")\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": 9,
|
|
"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, k_max = 1000000):\n",
|
|
" \"\"\"\n",
|
|
" Return the Conjugate Gradient algorithm estimate solution x to the problem\n",
|
|
" Ax = b and the number of iterations k it took to decrease maximum\n",
|
|
" norm error below eps or to exceed iteration maximum k_max.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # Assert n by n matrix.\n",
|
|
" assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
|
|
" \n",
|
|
" n = len(A)\n",
|
|
" \n",
|
|
" x_cur = np.zeros(n)\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" r_cur = b - A @ x_cur\n",
|
|
" v = r_cur\n",
|
|
" \n",
|
|
" k = 0\n",
|
|
" while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
|
|
" k += 1\n",
|
|
" x_prev = x_cur.copy()\n",
|
|
" r_prev = r_cur\n",
|
|
" \n",
|
|
" t = np.dot(r_prev, r_prev)/np.dot(v, A @ v)\n",
|
|
" x_cur = x_prev + t*v\n",
|
|
" r_cur = r_prev - t*A @ v\n",
|
|
" s = np.dot(r_cur, r_cur)/np.dot(r_prev, r_prev)\n",
|
|
" v = r_cur + s*v\n",
|
|
" \n",
|
|
" return x_cur, k"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"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": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting with A =\n",
|
|
"[[10 -1 2 0]\n",
|
|
" [-1 11 -1 3]\n",
|
|
" [ 2 -1 10 -1]\n",
|
|
" [ 0 3 -1 8]]\n",
|
|
"and b = [ 6 25 -11 15]\n",
|
|
"We apply the Conjugate Gradient algorithm to solve Ax = b.\n",
|
|
"\n",
|
|
"For eps = 1e-01\tafter k = 4\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-02\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-03\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n",
|
|
"For eps = 1e-04\tafter k = 5\t iterations:\n",
|
|
"x =\t\t\t [ 1. 2. -1. 1.]\n",
|
|
"Ax =\t\t\t [ 6. 25. -11. 15.]\n",
|
|
"diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
|
|
"diff(x, x_exact) =\t 2.220446049250313e-16\n",
|
|
"\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"def test_CG():\n",
|
|
" \"\"\"\n",
|
|
" Check that CG returns solutions for the example system Ax = b\n",
|
|
" within the error defined by the same eps as used for the iteration.\n",
|
|
" \"\"\"\n",
|
|
" \n",
|
|
" # TODO: Do we need to plot the error over epsilon like in task 4?\n",
|
|
" \n",
|
|
" return test_alg(CG, \"Conjugate Gradient\")\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": 11,
|
|
"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": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA93UlEQVR4nO3dd3hUVf7H8fdJJyEJhCQk9NBJA0LoAhFEUBEUUURUsGFZbPuzgLq2tcuuLpZVdkVEEVRAhQUVFTGggNIJNYAgoaQBCenJzPn9ccOQCglk5s4k39fzzJOZOyWfGcL5zj3n3nOU1hohhBCiLDezAwghhHA+UhyEEEJUIsVBCCFEJVIchBBCVCLFQQghRCVSHIQQQlQixUEIIUQlUhyEEEJU4mF2gIqUUgnA34EdwAKt9arzPSc4OFi3a9fOrrmEEKK+2bhxY4bWOqSq+xxSHJRSs4FRQJrWOrrM9pHAvwB34L9a61cADeQAPkBKTV6/Xbt2bNiwoc5zCyFEfaaUOlTdfY7qVpoDjCy7QSnlDrwDXAFEAhOUUpHAaq31FcDjwHMOyieEEKIMhxQHrXUicKLC5j7APq31Aa11EbAAGKO1tpbefxLwdkQ+IYQQ5Zk55tASOFzmdgrQVyk1FhgBNAHeru7JSqkpwBSANm3a2C+lEEI0QE43IK21XgwsrsHjZimljgFXe3l59bJ/MiGEaDjMPJT1CNC6zO1WpdtqTGu9VGs9JTAwsE6DCSFEQ2dmcfgd6KSUilBKeQE3AktMzCOEEKKUQ4qDUmo+sBboopRKUUrdobUuAaYC3wG7gM+11jtq+bpXK6VmZWVl1X1oIYRwRkV5cPIQpGyA3cvh1OHzP+cCqPqwElx8fLyW8xyEEC5JayjIgtx045KTVvl62W1FOeWfP/otiLv1gn61Umqj1jq+qvucbkC6NpRSVwNXd+zY0ewoQghxltUCeZnVN/YVr1uKqngRBb7NwC8EGodAy15nr/uFnr0e1N4ub0H2HIQQoiZKCmve2Odlgu2UrTLcPKtu4Ku67tsM3O37/b3e7jkIIcRFKzwNWUcgO8X4efpYFQ1/OhRWM7bp6Xe2UW8aAa37VN/w+zQBpRz69i6USxcH6VYSQpxTcT5kH4WsFMg+Ur4InLldVaPfqKnRqDcOhbCYc3zLDwEvP8e/LweQbiUhhGuyFBsNf5WNfmkxyMus/DzfYAhsCQGtSn+2hMBWpT9bQuMw8PBy/PsxgXQrCSFci9UCOanVN/pZR4z7qfDl1jvQaOgDWxoDuBWLQEBL8PQx5S25GpcuDtKtJIQL0tr4Rl9VV8+ZbaePgbWk/PM8fc9+u+94WZlv/WUKgLe/Oe+pHpJuJSFE3SvKhYxkyNhr/Mw6fLbhzz4KJQXlH+/uBQEtKnT1VPjW36ipywzmugrpVhJC1D2tjSN6MvZWuJQWgzOUG/iHGw18eHfociUEti7f3+8bDG6yarEzkeIghDg3SwmcOgTpe842/hml1wvKHOnj6QfBnaBNfwieBCGdIbizcZKWhyzN4mpcujjImIMQdagwBzKTjca/bCE4sb/8GbyNmxuNfvQ44+eZIuDfQr791yMy5iBEQ1KuK2hPmUKQbAwKn6HcISjCaPSDO0Fwl7PXGzUxLb6oWzLmIERDYymBkwerGA+opiuo3cDSItDZKARBEdIV1MBJcRDClZ3pCkqvUAAy94O1+OzjGocZjX/0OAjpcrYQBLSUI4BElaQ4COEqigsg5Tc48DMc2WgUgewyiyfauoK6QOeRpeMBXaBZR+kKErXm0sVBBqRFvWa1wLEtRjH442f4c51xfoByh7BoaHdJaTdQ2aOCGsa0D8L+ZEBaCGehtTE4/MfPRkE4uObspHChURAxGNoPgbYDwSfA3KyiXpABaSGc1anDZ4vBH4mQc9zY3qQtRI2BiCFGUWgcam5O0eBIcRDCkXIzjCJwpiCc/MPY7hdiFIGIIcbeQdN2psYUQoqDEPZUmAOHfj1bDFK3G9u9/I0xg753GwUhtJscNSScihQHIepSSSGk/H52EPnIRmN2UXdvY4WwoU9BRAK06Gn3JSCFuBjy1ynExbBa4Pi2s8Xg0FooyTcmm2vREwY8YHQTte4Lno3MTitEjbl0cZBDWYXDaW1MNfHHz3BglXFEUcEp476QrhB369kjiuTcAuHC5FBWIc4n60j5I4pOHzW2B7aB9oONbqKIweDf3MyUQtSaHMoqRG3knYCDq892FWXuM7b7NqtwRFGEDCKLekuKgxBgHFG0Z7lREI5vBzR4NTa6h+JvLz2iKFKmpBYNhhQHIX75F3z/tLFUZas+cOkTRjFoGQfunmanE8IUUhxEw6U1rHwBVs+AqLEw5m3w8jM7lRBOQYqDaJisVvh2Gvz2vnGE0ag3wc3d7FRCOA0pDqLhsZTAkvth66fQfypc/oIMLAtRgRQH0bCUFMKiO2DXUrj0SRj8qBQGIarglIdeKKX8lFIblFKjzM4i6pGiXJh/o1EYRr4CQx6TwiBENRxSHJRSs5VSaUqppArbRyql9iil9imlppW563Hgc0dkEw1E/in4eKxxVvOYd6DfvWYnEsKpOWrPYQ4wsuwGpZQ78A5wBRAJTFBKRSqlhgM7gTQHZRP1XU46fDTKmARv3IfQ82azEwnh9Bwy5qC1TlRKtauwuQ+wT2t9AEAptQAYAzQG/DAKRr5SarnW2lrxNZVSU4ApAG3atLFjeuHSso7A3DGQlQITFkCny8xOJIRLMHNAuiVwuMztFKCv1noqgFJqMpBRVWEA0FrPAmaBMbeSfaMKl5S5H+ZeY0yMd8tiaDvA7ERCuAynPVpJaz3nfI+RWVlFtVJ3wsfXGGspTFoKLXqYnUgIl2Lm0UpHgNZlbrcq3VZjWuulWuspgYGBdRpMuLiUjTDnSmNNhdu+kcIgxAUwszj8DnRSSkUopbyAG4EltXkBpdTVSqlZWVlZdgkoXNAfq2HuaPAJhNu/hZAuZicSwiU56lDW+cBaoItSKkUpdYfWugSYCnwH7AI+11rvqM3ryp6DKGfPt/DJdRDYCm77Fpq2MzuREC7LUUcrTahm+3JguSMyiHpu+0L48m4Ii4GJi8CvmdmJhHBpTnmGdE1Jt5IAYMOHsOhOY53mW5dIYRCiDrh0cZBuJcEvM+F/D0Gn4XDzIvAJMDuREPWCSxcH2XNowM6sxfD93yDqWhg/DzwbmZ1KiHrDpYuD7Dk0UFYrfPM4JL5urMVw3Qfg4WV2KiHqFac9CU6IKllKYOkDsGWerMUghB1JcRCuo6TQGHjetQQSnpApt4WwI5cuDjJ9RgNSlAuf3Qz7VxprMciU20LYlYw5COdXkHV2LYbRb0thEMIBXHrPQTQAuRnwyVhjIr1xs40jk4QQdufSew5yKGs9l30UPrwC0vcaazFIYRDCYVy6OEi3Uj124gDMHgHZx4y1GGSRHiEcSrqVhPNJ3QkfXwuWIpi0BFrGmZ1IiAbHpfccRD10pHQtBjDWYpDCIIQppDgI5/HHavhoNHgHGGsxhHY1O5EQDZYUB+Ec9n4H88YZazHc/i0ERZidSIgGzaWLgxytVE8kLYIFN0FIV5i8HAJamJ1IiAbPpYuDHK1UD2ycAwvvMNZimLRU1mIQwkm4dHEQLu7Xt2Dpg9DxMpi4UNZiEMKJyKGswvG0hp9egsTXIPIaGPsfmXJbCCcjxUE4ltUK302H9e9Bz1vg6n+Bm7vZqYQQFUhxEI5Tdi2Gfn+BES/KlNtCOCkpDsIxSopg8Z2w82tImA5DHpfCIIQTc+niIOs5uIiiPPj8Ftj3A4x4Cfr/xexEQojzcOmjleRQVhdQkAWfXAf7foTRb0lhEMJFuPSeg3ByuZnwybWQusNYiyF6rNmJhBA1JMVB2Ef2MZg7Bk4dghvnQ+fLzU4khKgFKQ7CPpY/AlkpcPMiaHeJ2WmEELXk0mMOwkkd2QS7/weXPCSFQQgXJcVB1L2fXoRGQdD3HrOTCCEukBQHUbcOrTUOWb3kIZkrSQgX5nTFQSnVTSn1nlJqoVLqXrPziFrQGla+AH6h0Psus9MIIS6CQ4qDUmq2UipNKZVUYftIpdQepdQ+pdQ0AK31Lq31PcANwEBH5BN15I+f4dAaGPwIePmanUYIcREctecwBxhZdoNSyh14B7gCiAQmKKUiS+8bDSwDljson7hYZ/YaAlpBr8lmpxFCXCSHFAetdSJwosLmPsA+rfUBrXURsAAYU/r4JVrrK4CJjsgn6kDyCkj53dhr8PA2O40Q4iKZeZ5DS+BwmdspQF+lVAIwFvDmHHsOSqkpwBSANm3a2C2kqIEzew1N20HPm81OI4SoA053EpzWehWwqgaPmwXMAoiPj9f2TSXOaddSOL4NrnkP3D3NTiOEqANmHq10BGhd5nar0m01ppS6Wik1Kysrq06DiVqwWoxV3YI7Q+wNZqcRQtQRM4vD70AnpVSEUsoLuBFYUpsXkFlZnUDSYkjfZazRICu6CVFvOOpQ1vnAWqCLUipFKXWH1roEmAp8B+wCPtda76jl68qeg5ksJbDqJWgebawFLYSoNxwy5qC1nlDN9uVcxOGqWuulwNL4+Hg548oMW+fDiQPGrKtuTnc+pRDiIrj0/2jZczBRSSH8/Cq0iIMuV5idRghRx1y6OMiYg4k2zYWswzD0KVkLWoh6yKWLgzBJcT4kzoA2A6DDULPTCCHswKWLg3QrmeT3DyDnOAx9UvYahKinXLo4SLeSCQpzYM0b0D5BFvIRoh5z6eIgTPDb+5CXAZc+ZXYSIYQduXRxkG4lB8s/Bb/8CzqPhNa9zU4jhLAjly4O0q3kYOvehYIsuPQJs5MIIezMpYuDcKDcTFj7LkSOgfDuZqcRQtiZFAdRM7/+C4pyIEH2GoRoCFy6OMiYg4OcToX1s4xZV0O7mp1GCOEALl0cZMzBQdb8EyxFMORxs5MIIRzEpYuDcICsFNgwG3rcBM06mJ1GCOEgUhzEuSXOMJYBHfKY2UmEEA4kxUFU78QfsPlj6DUZmsg63UI0JC5dHGRA2s5+fg3cPGDQ/5mdRAjhYDUuDkqpF6rYZuq6kDIgbUfpe2HbAuh9JwSEm51GCOFgtdlzaKmUsq3oppQKBX6o+0jCKax6GTwawSUPm51ECGGC2iwTejfwnVJqP6CBDwE5trE+Op4EOxbDoEfAL9jsNEIIE5y3OCil5gKbgM3AX4BPgRLgGq31PvvGE6b46SXwDoQBU81OIoQwSU26leYACrgN+ARoB5wEblZKjbNbMmGOIxthzzIYcD80amp2GiGESc6756C1XgmsPHNbKeUBdAO6A32BhXZLJxxv5YvQKAj63WN2EiGEiWoz5gCA1roE2F56+aTOE9WCUupq4OqOHTuaGaP+OLQW9v8Iw58Hb3+z0wghTOTS5znIoax1SGtY+QL4hULvu8xOI4QwmUsXB1GH/vgZDq2BwY+Al6/ZaYQQJpPiIM7uNQS0MqbKEEI0eFIcBCSvgJTfYcij4OFtdhohhBOQ4tDQWa3GXkPTdtBjotlphBBOotZHK4l6ZvdSOL4Nrn0f3D3NTiOEcBKy59CQWS3G2dDBnSHmerPTCCGciOw5NGRJiyB9N1w/B9xMnWBXCOFknLI4KKWuAa4CAoAPtNYrzE1UD1lKjJlXm8dAtzFmpxFCOBmHdSsppWYrpdKUUkkVto9USu1RSu1TSk0D0Fp/pbW+C7gHGO+ojA3K1vlw4gBc+gS4Se+iEKI8R7YKc4CRZTeULhb0DnAFEAlMUEpFlnnIU6X3i7pUUgg/vwot4qDLFWanEUI4IYcVB611InCiwuY+wD6t9QGtdRGwABijDK8C32itNzkqY4OxaS5kHYahT4FSZqcRQjghs/sTWgKHy9xOKd12P3AZME4pVeX0oEqpKUqpDUqpDenp6fZPWl8U50PiDGgzADoMNTuNEMJJOeWAtNZ6JjDzPI+ZBcwCiI+P147IVS/8/gHkHIdxs2WvQQhRLbP3HI4ArcvcblW6rUaUUlcrpWZlZWXVebB6qTAH1vwT2l8K7QaanUYI4cTMLg6/A52UUhFKKS/gRmBJTZ8sU3bX0vr3IC/TGGsQQohzcOShrPOBtUAXpVSKUuqO0oWDpgLfAbuAz7XWO2rxmrLnUFP5p+DXmdD5CmgVb3YaIYSTU1q7fnd9fHy83rBhg9kxnNtPLxmHr969GsJjzU4jhHACSqmNWusqvy2a3a10UWTPoYZyM2HtuxB5jRQGIUSNuHRxkDGHGvr1X1CUAwnTzU4ihHARTnkoq6hDp1Nh/SyIvQFCu5qdRghRBatVk1tUQm6hhZzCEnIKS8gt89O4brFtK3v/lMHtGdQppM4zuXRxUEpdDVzdsWNHs6M4rzX/BEsRDHnc7CRC1CuFJRZyCs426LlF5Rvz06X3VdxuXC9fBPKKLDX6nW4KGnt70NjbA7/SS1GJ1S7vz6WLg9Z6KbA0Pj7+LrOzOKWsFNgwG3pOhGYdzE4jhKlKLFajUS4621DbGvbCEvKKyn87zy1t8I1GvnwRyC0sodhSs4N5vD3c8Pcpbcy9jIY9uLEXbZv5Gtu9jPvONvju5QpA2Z8+nm4oB5286tLFQZxH4uvGz8GPmZtDiAtQbLHaGuq8Iku5b9+Vv5FbbI15TqGFvDMNfJmumpp+w3ZTGI24z9lv54293Qlu7G002j5lGm0v9/INuE+Z615GQ+/h7ppDuy5dHKRb6RxOHIDNn0D87dCk9fkfL0Qd0FqTW2QhK7+YrLxiTuUXkZ1fvo88r6j8N/bcIku5LpczhaBWjbn3mW/g7rbGuZWvL429yzfevl7u5b6J+1Zo3P283Wnk6e6wb+cXIy0vje3p24kKjiLML6zOX9+li4N0K53Dz6+BmwcM+j+zkwgXo7Umv9ho4E/lFRsNfeklu/RnVdtPlf4ssZ67u+VMY16239zPy50gP9/Sbe62rpYz39p9vSp3u5zZ5siuFrMUlBSw68QutqVvMy4Z2zieexyAp/o+xfiudb/sjUsXB1GN9L2w7TPo/xfwr/tvFMI1FJQ28FU15ll5ReVvl2ncs/KLz9mfrhQE+HjSxNeTwEbGpWXTRgQ28qRJo7PbzlwCGnna+twbe3vg7VH/G/OLobXmUPYhtmVssxWD5JPJlOgSAFo2bknPkJ7ERMYQGxJL1yD7HIUoxaE+WvUyePrCwIfMTiLqQEGxhRO5ReUuJ0sb91N5Zxv0so18Vn7xebtlAnw8CCzTwIcF+hDYyKtS4162CAQ08sTf2wM3N2nc60pWYRbbM7bb9gi2p28nuygbAD9PP6KDo7kt+jZigmOICYkhuFGwQ3K5dHGQMYcqHE+CHYth0CPg55g/IlFzVqsmu6C4UmOfmVvEyTO388rfd67DHP29PQgo05B3DG1criGv2Lifufj7eOIuDbzDFVuL2XtyL9vSjSKwLWMbh7IPAeCm3OjQpAPD2w4nNiSW2OBYIgIjcHdzNyWrzK1U38y/CQ6ugYe2QqOmZqep9wpLLJzMLdPY5xVxIqewykbe+MZfjKWaPvlGnu4E+XnV6NLU14sAHw+XPRKmIdBaczz3uK17aHvGdnZm7qTQUghAcKNgYoNjiQmJITY4lqjgKPw8/Rya8VxzK7n0noOo4MhG2LMMLn1KCsMF0FpzurCEEzlnGvmqG/iyl5zCkipfSylo0sjT1phHBPvRq23T0tveBPl5Gj99vQhq7EWQrxeNvMz5hijqRl5xHjsyd9jGCbZnbCc931il0tvdm25B3RjfZTwxITF0D+5OmF+YU4+9SHGoT1a+CI2CoF+VK6s2WBarJjO3kLTsQtJOF5CabVxPPV1AWnYBaacLSc0u4ERuUbUDsV4ebjQr8829bTNfmvp60czPi6Z+XuXuC/Lzoomvl3Tb1GNWbeXAqQPGWEHpnsG+U/uwamOcp21AW/qG97V1D3Vu2hlPd0+TU9eOFIf64tCvsP9HGP538PY3O41DWK2aE3lFpGYXGI19mYY+NbuQ9NJCkJ5TWGVXTpCfF6H+3jQP8KFLc3+aNfautrH39XKNY9+FfWTmZ5YbNN6RsYOc4hwA/L38iQ2OZVibYcagcXAMTXyamBu4Drh0cZAB6VJaw8oXoHFz6H2n2WkumtWqOZlXRGrZb/el11OzC0k7XUhadgHppwurPKa+qa8nzQN8CA3woVNzf5oHGAUg1N+H0NLrIY298fKQ/npRWZGliN0ndtsKwbb0bRzJMVYvdlfudG7amavaX0VsSCwxwTG0DWiLm6p/f0suXRzkJLhSB1bBoV/gitfBy9fsNNXSWnMyr7j0m32BrZFPLfOt/0w3T1WNfhNfT5qXNvAdQ4JtjX7zAG9C/M/89MbbQ/ruzZaRn8G6Y+tYd3QdG1M3UmApMDtSjWUVZlFsLQYgzC+MmOAYbuxyI7EhsXRr1o1GHo1MTugYLl0cBGf3GgJaQa9JpsUoLLFwPKuAY1kFHMvK51hWAcezCsr07Rv9/VX16Qc28rQ19O1DmhkNfml3T2iAN6H+PoT4e+PjKY2+s8orzmND6gajIBxbR/LJZAACvALoE9bHpbpZArwCbEcRhfqGmh3HNFIcXF3yCjiyAa6eCR7edvkVhSUWUrMKbY1+2QJwLCufY6cKyMwtqvS8AB+P0m/2PvSN8CO09Fu+0cVT2r0jjb5LKrGWsCNzB+uOrmPtsbVsTd9KibUELzcvejbvyYNxD9K/RX+6Nu1q2nH64uJIcXBlVqux19A0AnrcdEEvUVRiJTW7QoN/qmwRKCAjp7DS8wJ8PAgPbER4Ex9iWjYhPNCn9GJsCwvwwc9b/rzqC601B7MP2rqKfj/+O6eLTwPQLagbt0TeQr/wfsSFxuHj4WNyWlEX5H+vK9u9FI5vg2vfhyoOkyu2nG34j57Kr9Ttc6bhr3gepL+Ph62hj2oRYGvwbY1/oDT8DUFGfgbrj623dRWdmeithV8LLm93Of1a9KNvWF+a+sg5NfWR/A93UcXFxfDDC5QEduAHPZBjifs5eqqgtAAYjX96VQ2/twdhgT6EN2lEt7AAwpv40CKwEWGBPrRo4kNYYCMaS8PfIOUV57EpbZOtq2jvyb2A0QffN7wvd8XcRf/w/rTybyWH9TYALt0K1PdDWbXWpJ0uZOfRbHYey2bXsWwOn8zneFY+A3J/5A3PvTxY9ADLF2wDwM/LnfAmxjf7LmH+tm/54U0a0SLQh7BAH/x9XOtEHGE/FqvFGDco3TPYkraFYmsxnm6exIXG8WDcg/QL70e3oG4ybtAAydxKTqLEYuVARi47jxpFYOexbHYezS430Ns6qBHtmvnRwt+DaQcmgacvW69aQngTP8Kb+BAgDb84B601f57+k7VH17Lu2Dp+O/4bp4uMcYOuQV3pF96P/uH96dm8Z4M5XLOhk7mVnMzpgmJ2Hz9tFIHSvYLdx0/bplj2cnejc1hjhnULJTI8gMgWgXQN9z/b+G+aCzsPw7ULSOgi6zWI6p0oOGEbN1h7dC3Hco8BEO4XzvC2w+kX3o++4X0J8gkyOalwNlIc7EhrzbGsgnJFYOexbA5l5tke09TXk8gWAUzq35bIFgFEhgfSPsQPz+pm2ywpNFZ5a9kLOo900DsRriK/JJ9NqZtsXUW7T+wGjCke+ob15Y7oO+jXoh9t/NvIuIE4JykOdaTYYmVfWk6lQnAqr9j2mHbNfIlqEcD1vVoR2SKAbuEBhAX41O4/6aa5kHUYRs80pv4UDZrFamHXiV22rqLNaZspthbj4eZBz9Ce3N/zfvqH9yeyWaSMG4hakeJwAbLyi9ldZlxg57FsklNzKLIY3ULeHm50DfPniugwIsONItA1PODijwIqzofEGdB2ILS/tA7eiSi2FnM85zgpOSm2KRNcwfHc46w9upb1x9fbxg26NO3CTV1vol8L43wDX0/nnUpFOD8pDuegtSblZH65AeKdx7JJOZlve0wzPy8iWwRw2yXtjPGB8AAigv3sswjL7x9AznEYN1v2GmohrziPw6cPk3I6hcOnD5e7HMs9hkVXv9KaM2vu25xhbYbRP7w/fcL7OGz5SNEwSHEoVVRiJTnt9NkuodKjhrILjMVclIKIYD+6t27ChD5tiGwRQFR4ACH+3o7pu83cD6v/AR2GQruB9v99LkRrzYmCE7YGv2IRyCzILPf4QO9AWjduTUxwDFdEXEFr/9a08m/lUkfoBHgF0Nq/tYwbCLtp0MXh573pfL3lCDuPZrMvLcc2E2gjT3e6hPkzqnuL0qOFAuga5o+vl0kf16k/Ye4Yo0KNfNWcDCYrsZZwPPd4tQUgr+TsIL9C0dyvOa39WzOk9RBb49/avzWt/VsT4BVg4jsRwjU4XXFQSrUHngQCtdbj7Pm7klNPszo5g8jwAC7tGmorBO2a+TnPKl7Zx+Cj0VCYDZOWQkhnsxPZTX5JfqVG/8ztozlHKdFnl+T0dPO0NfjxYfG2hr+VfytaNm6Jt7t9JiEUoqFwyElwSqnZwCggTWsdXWb7SOBfgDvwX631K2XuW1jT4nChJ8FZrNp5ikBVcjPgwysh+wjc8hW07m12oouitSarMKtSv/+ZIpCWn1bu8f5e/rZGv+Il1De0Xi6wUheKi4tJSUmhoMB11lAQ9uXj40OrVq3w9Cx/oqwznAQ3B3gbmFsmlDvwDjAcSAF+V0ot0VrvdFAm5y4M+Sdh7jVGl9LNC12qMBRZithzYg97T+6tVADOzOR5RmijUFr5t2JAywGVCkCgd6BJ78C1paSk4O/vT7t27WRMQqC1JjMzk5SUFCIiImr8PIcUB611olKqXYXNfYB9WusDAEqpBcAYoEbFQSk1BZgC0KZNm7oL6wwKsuGT6yBjD0yYD+0uMTtRtSxWC/uz9rMjYwdJGUkkZSax9+ReSqxGF5CHmwctG7eklX8ruod0L9f4t/JvJdM720FBQYEUBmGjlKJZs2akp6fX6nlmjjm0BA6XuZ0C9FVKNQNeBHoqpaZrrV+u6sla61nALDC6lewd1mGK8uDT8XB0C4z/GDpeZnYiG601KadTSMpMMgpBRhK7Tuwiv8Q4tLexZ2OigqOYFDmJ6OBougZ1JdwvXE6+MoEUBlHWhfw9ON2AtNY6E7inJo+td7OyFhfAgpvg8Dq47r/Q9SpT46Tnpdv2BnZk7CApM4mswiwAvN296RrUlbGdxhLVLIro4Oh6u9C6EA2RmcXhCNC6zO1WpdtqTGu9FFgaHx9/V10GM4WlGL6YDAd+gjHvQvR1Dv312UXZ7MjYwY5Mo3toe8Z20vKMAWJ35U7HJh25rM1lRAVHEd0smo5NO+LpJrPAiuqlpqby8MMPs27dOpo2bYqXlxePPfYY1157rcMyrFu3jgcffJDCwkIKCwsZP348zz77bLWP37BhA3PnzmXmzJmV7mvXrh0bNmwgOLj6kw1r8hhXYWZx+B3opJSKwCgKNwK1Wuuy3uw5WC2w+C7Y+w1cOQN6TrTrr8svyWfPiT1sz9hOUkYSOzJ3cCj7kO3+tgFtiW8eT3RwtK17yJVOEBPm01pzzTXXMGnSJD799FMADh06xJIlSxyaY9KkSXz++ed0794di8XCnj17zvn4+Ph44uOrPHinwXFIcVBKzQcSgGClVArwjNb6A6XUVOA7jENZZ2utd9TmdevFnoPVCl9PhR1fwvC/Q5+6fSvF1mL2n9rP9ozttkHjfaf22aaMCPUNJbpZNGM6jCEqOIqoZlFylJC4aCtXrsTLy4t77jnbQ9y2bVvuv/9+Dh48yC233EJubi4Ab7/9NgMGDGDVqlXMmDGD//3vfwBMnTqV+Ph4Jk+ezLRp01iyZAkeHh5cfvnlzJgxgy+++ILnnnsOd3d3AgMDSUxMrJQjLS2N8PBwANzd3YmMjAQgNzeX+++/n6SkJIqLi3n22WcZM2ZMuQyZmZlMmDCBI0eO0L9/f8oe9v/JJ58wc+ZMioqK6Nu3L++++y7u7vVrbM1RRytNqGb7cmC5IzI4Ja1h+SOw9VNIeAIGPnBRL2fVVg5lH7LtDSRlJLH7xG4KLYWAMeVCdHA0g1sNtu0VhPqG1sU7EU7quaU72Hk0u05fM7JFAM9cHXXOx+zYsYO4uLgq7wsNDeX777/Hx8eH5ORkJkyYwLnOU8rMzOTLL79k9+7dKKU4deoUAM8//zzfffcdLVu2tG2r6OGHH6ZLly4kJCQwcuRIJk2ahI+PDy+++CJDhw5l9uzZnDp1ij59+nDZZeUP/njuuee45JJLePrpp1m2bBkffPABALt27eKzzz7jl19+wdPTk/vuu4958+Zx6623nvMzcTVONyBdGy7draQ1rHgKNnwAAx+EIY/V8uma1LxU21FDSZlJ7MzYaTuPoJFHI7oFdeOGLjcQ3SyamOAYWftXmOYvf/kLa9aswcvLix9++IGpU6eyZcsW3N3d2bt37zmfGxgYiI+PD3fccQejRo1i1KhRAAwcOJDJkydzww03MHbs2Cqf+/TTTzNx4kRWrFjBp59+yvz581m1ahUrVqxgyZIlzJgxAzAO//3zzz/LPTcxMZHFixcDcNVVV9G0aVMAfvzxRzZu3Ejv3sa5R/n5+YSG1r8vWS5dHFy6W2nVy7D2begzBS57rkazrG5J28K6Y+tsBeHMhHIeyoNOTTsxMmKkbY+gfWB7PNxc+p9X1IHzfcO3l6ioKBYtWmS7/c4775CRkUF8fDxvvPEGzZs3Z+vWrVitVnx8jHNdPDw8sFqttuecOcPbw8OD3377jR9//JGFCxfy9ttvs3LlSt577z3Wr1/PsmXL6NWrFxs3buSRRx5h8+bNtGjRguXLjU6JDh06cO+993LXXXcREhJCZmYmWmsWLVpEly5dyuVOTU0973vTWjNp0iRefrnKo+zrDZc+7lApdbVSalZWVpbZUWpnzRvw86vQ82ZjIr3zFIaCkgL+vvbv3PLNLby75V3+PP0nA1oMYFqfacy7ch7rJq7j86s/5+n+TzO201g6N+0shUGYaujQoRQUFPDvf//bti0vz5gcMSsri/DwcNzc3Pj444+xWIzxr7Zt27Jz504KCws5deoUP/74IwA5OTlkZWVx5ZVX8sYbb7B161YA9u/fT9++fXn++ecJCQnh8OHDfPjhh2zZssVWGJYtW2YbK0hOTsbd3Z0mTZowYsQI3nrrLdt9mzdvrvQeBg8ebBtM/+abbzh58iQAw4YNY+HChaSlGUfznThxgkOHDlV6vqtz6RbEJfcc1r8PPzwL0ePg6pngdu76vO/kPh5NfJR9p/YxOWoyU2Kn4O/l75isQlwgpRRfffUVDz/8MK+99hohISH4+fnx6quvEhcXx3XXXcfcuXMZOXIkfn5+ALRu3ZobbriB6OhoIiIi6NmzJwCnT59mzJgxFBQUoLXmn//8JwCPPvooycnJaK0ZNmwY3bt3r5Tj448/5uGHH8bX1xcPDw/mzZuHu7s7f/vb33jooYeIjY3FarUSERFhGwg/45lnnmHChAlERUUxYMAA20wMkZGRvPDCC1x++eVYrVY8PT155513aNu2rT0/UodzyMR79nahE+853Ka5sOR+6DoKrp8D7tWfJ6C15ou9X/Da76/h5+nHS5e8xMCWso6DOL9du3bRrVs3s2MIJ1PV34UzTLwntn0BSx4wpsMYN/uchSGrMIvn1j7H94e+p394f14a9JKs8iWEcCiXLg4uc7TSrqXw5d3GBHrjPwGP6tca2JS6icdXP05GXgZ/7fVXJkVNkikphBAO59KtjtZ6qdZ6SmCgE5+0lfw9fHEbtOxlzLDqWfWZxharhfe2vsdt392Gh/Lg4ys/5rbo26QwCCFM4dJ7Dk7vj0T47GZoHgkTvwDvqgeSj+ceZ/rq6WxI3cCVEVfyt35/o7FXYweHFUKIs6Q42Muf6+HTG6FpBNz8JTRqUuXDfvrzJ/72698oshTxwsAXGN1htJyoJoQwnUsXB6cdczi6GeaNA/8wuPUr8GtW6SGFlkL+seEfzN89n25B3Xht8Gu0C2zn8KhCCFEVl+7Qdsoxh9Sd8PG14NMEJi0xCkQFB04d4KZlNzF/93xuibyFT678RAqDqFdefPFFoqKiiI2NpUePHqxfvx6AN99803YynFlOnTrFu+++W+397u7u9OjRg6ioKLp3784//vGPcmdum+Gll15y/C/VWrv8pVevXtoppCdr/VpHrWd00Tpzf6W7rVarXrR3ke79SW89aP4g/fPhn00IKeq7nTt3mvr7f/31V92vXz9dUFCgtdY6PT1dHzlyRGutddu2bXV6erqZ8fQff/yho6Kiqr3fz8/Pdj01NVUPGzZMP/30046IVq2ymS5UVX8XwAZdTbvq0nsOTuXkIZg7GrQVbl0CQe3L3Z1dlM1jiY/xzK/PEBscy8LRCxncarBJYYWwn2PHjhEcHIy3t3HIdnBwMC1atGDmzJkcPXqUSy+9lEsvvRSAFStW0L9/f+Li4rj++uvJyckBYOPGjQwZMoRevXoxYsQIjh07BkBCQgIPPvggPXr0IDo6mt9++w0wpuC+/fbb6dOnDz179uTrr78GjNlh+/TpQ48ePYiNjSU5OZlp06axf/9+evTowaOPPnrO9xIaGsqsWbN4++230VpjsVh49NFH6d27N7Gxsbz//vu29zx48GBbrtWrVwPw7bffEhcXR/fu3Rk2bNg5s86ZM4exY8cycuRIOnXqxGOPGZNxTps2jfz8fHr06MHEifZd66UsOUO6LmQfhdkjoSALJv8PwmLK3b0lbQvTVk/jeO5xpvacym1Rt8m6ysJuyp0J+800OL69bn9BWAxc8Uq1d+fk5HDJJZeQl5fHZZddxvjx4xkyZAhQfqW0jIwMxo4dyzfffGObWqOwsJDp06czZMgQvv76a0JCQvjss8/47rvvmD17NgkJCXTq1In//Oc/JCYmct9995GUlMQTTzxBZGQkN998s20K7s2bNzNt2jT69evHxIkTKSoqwmKxkJqayqhRo0hKSqoyf+PGjW1F6owmTZqwZ88evv76a9LS0njqqacoLCxk4MCBfPHFFyxevJiCggKefPJJLBYLeXl5FBQUEBcXR2JiIhEREZw4cYKgoKBqs37xxRc8//zzbN68GW9vb7p06cKaNWto3bp1lZlqq0GdIe0UA9I5afDRaMg7AZO+LlcYLFYLH+74kLc3v02YXxgfXfER3UMqz/8iRH3SuHFjNm7cyOrVq/npp58YP348r7zyCpMnTy73uHXr1rFz504GDjSmhSkqKqJ///7s2bOHpKQkhg8fDoDFYrEt2AMwYYKxPMzgwYPJzs7m1KlT1U7B3b9/f1588UVSUlIYO3YsnTp1uqj3tmLFCrZt28bChQsBYxLB5ORkevfuze23305xcTHXXHMNPXr0YNWqVQwePJiIiAgAgoKCbK9R3XThw4YN48wYamRkJIcOHaJ169YVYziESxcHbfbEe3knYO41kH0Ebl5snOhWKi0vjSdWP8H64+sZ2W4kT/d/WibME453jm/49uTu7k5CQgIJCQnExMTw0UcfVSoOWmuGDx/O/Pnzy23fvn07UVFRrF27tsrXrniot1Kq2im4u3XrRt++fVm2bBlXXnkl77//Pu3bl+/yPZ8DBw7g7u5OaGgoWmveeustRowYUelxiYmJLFu2jMmTJ/PXv/7Vtv5DRdVlXb9+va0rDozPsKSkpFZZ65KMOVyogiz4ZCxk7oMbP4W2/W13JaYkMm7JOLZlbOP5Ac/z2uDXpDCIBmPPnj0kJyfbbm/ZssU2Y6m/vz+nTxsLUvXr149ffvmFffv2AUZf/N69e+nSpQvp6em24lBcXMyOHWdXEP7ss88AWLNmDYGBgQQGBlY7BfeBAwdo3749DzzwAGPGjGHbtm3lMpxPeno699xzD1OnTkUpxYgRI/j3v/9NcXExAHv37iU3N5dDhw7RvHlz7rrrLu688042bdpEv379SExM5I8//gCMqb2BGk0XXpGnp6ftdzqKS+85mKYoF+bdYPTljp8HHYzBtSJLEW9sfINPdn1C56adeX3w67RvUrtvKUK4upycHO6//35OnTqFh4cHHTt2ZNasWQBMmTKFkSNH0qJFC3766SfmzJnDhAkTKCw0lrJ94YUX6Ny5MwsXLuSBBx4gKyuLkpISHnroIaKijIWLfHx86NmzJ8XFxcyePRug2im4P//8cz7++GM8PT0JCwvjiSeeICgoiIEDBxIdHc0VV1zB66+/Xi7/mcHf4uJiPDw8uOWWW/jrX/8KwJ133snBgweJi4tDa01ISAhfffUVq1at4vXXX8fT05PGjRszd+5cQkJCmDVrFmPHjsVqtdqWR63JdOEVTZkyhdjYWOLi4pg3b16d/ntVRwaka6u4AD69AQ6uNmZXjboWgINZB3ks8TF2ndjFTV1v4q/xf8XbvfoJ9oSwl/o8ZXdCQgIzZswgPr7KMVRxDg1qQNrhSorg81uNOZOufQ+irkVrzdf7v+al9S/h7e7NW0PfIqF1gtlJhRDiokhxqClLCSy6A5K/g1FvQPcbySnK4e/r/s7yP5bTO6w3L1/yMs39mpudVIh6a9WqVWZHaDBcujg47FBWqxW+vg92LYERL0P87WxP385jiY9xLPcYU3tM5c6YO+XcBSFEveHSRytpR8ytpDUsexi2fQZDn8La7x5mJ83m1m9uxaItfDjyQ+7ufrcUBiFEveLSew52pzV8Ox02zoFB/0dG79t44vt7WHtsLcPbDueZ/s8Q6O1Ek/4JIUQdkeJwLiv/Duv/DX3v5ZcuQ3liyXXkFufydP+nGddpnKy7IISot1y6W8muEl+H1f+gOO5W/hEczD0/3kuQTxALrlrA9Z2vl8IgxDkcP36cG2+8kQ4dOtCrVy+uvPJK9u7de0Gvdeedd7Jz5846zXe+aburk5OTw7333kuHDh2Ii4ujV69e/Oc//7moLHPmzGHq1KkAvPfee8ydO/eCXufgwYN8+umnF5WlLCkOVVn7Lqx8gT+jR3OLWzpzdn7E+C7jmX/VfDo2dbKFhYRwMlprrr32WhISEti/fz8bN27k5ZdfJjU19YJe77///S+RkZF1mvFCi8Odd95J06ZNSU5OZtOmTXz77be2M5/LutBpL+655x5uvfXWC3quFAd72zAbvpvO0s6XcH3hXg6fPsybCW/yVL+n8PHwMTudEE7vp59+wtPTk3vuuce2rXv37gwaNAitNY8++ijR0dHExMTYpsJYtWoVCQkJjBs3jq5duzJx4kTb9BIJCQmcOcm1ceOza6svXLjQNl/T/v376devHzExMTz11FO2x+Xk5DBs2DDi4uKIiYmxTY9d1bTdr7/+um0q7meeeabS+9q/fz+//fYbL7zwAm5uRtMZEhLC448/bnsPgwYNYvTo0bZids0119CrVy+ioqJsZ4kDfPjhh3Tu3Jk+ffrwyy+/2LY/++yztgn59u/fz8iRI+nVqxeDBg1i9+7dAEyePJkHHniAAQMG0L59e9skgNOmTWP16tX06NGDN954o5b/apXJmENZWxeQu+z/eKl9DEuK/yQuNI5XBr1CeOPw8z9XCCf06m+vsvvE7jp9za5BXXm8z+PV3p+UlESvXr2qvG/x4sVs2bKFrVu3kpGRQe/evRk82FjXZPPmzezYsYMWLVowcOBAfvnlFy655JIaZXrwwQd58MEHmTBhAu+9955tu4+PD19++SUBAQFkZGTQr18/Ro8ezSuvvEJSUhJbtmwBjJlSk5OT+e2339BaM3r0aBITE23ZwFgbonv37rbCUJVNmzaRlJRkm4l19uzZBAUFkZ+fT+/evbnuuusoKirimWeeYePGjQQGBnLppZfSs2fPSq81ZcoU3nvvPTp16sT69eu57777WLlyJWCsH7FmzRp2797N6NGjGTduHK+88gozZsw471QcNSXF4YwdX7Jj+QM81q4dKZzmvu73cVfsXXi4yUckRF1Zs2YNEyZMwN3dnebNmzNkyBB+//13AgIC6NOnD61atQKgR48eHDx4sMbFYe3atXz11VcA3HTTTTzyyCOA0cX1xBNPkJiYiJubG0eOHKmye2vFihWsWLHC1kjn5OSQnJxcrjhU9OKLL/LFF1+QlpbG0aNHAejTp4+tMADMnDmTL7/8EoDDhw+TnJzM8ePHSUhIICQkBIDx48dXGo/Jycnh119/5frrr7dtOzP/FBh7JG5ubkRGRl5wd935OF3Lp5TyA94FioBVWmu7zzJl3b2cj1c8yJvhzWnWKIgPBr9KfJjM3SJc37m+4dtLVFSUraujNmoyXXXZA0EKCgrO+5rz5s0jPT2djRs34unpSbt27ap8ntaa6dOnc/fdd1f7WpGRkWzduhWr1YqbmxtPPvkkTz75ZLmuLj8/P9v1VatW8cMPP7B27Vp8fX1JSEioUWYAq9VKkyZNbHs2FZX9rOw1P55DxhyUUrOVUmlKqaQK20cqpfYopfYppaaVbh4LLNRa3wWMtne2zF1L+Muqh5gRFMjgVoNYNGaxFAYhLsLQoUMpLCws18e+bds2Vq9ezaBBg/jss8+wWCykp6eTmJhInz59avzazZs3Z9euXVitVts3cjCm/160aBEACxYssG3PysoiNDQUT09PfvrpJw4dOgRQadruESNGMHv2bNtqa0eOHCEtLa3c7+7YsSPx8fE89dRTWCwWwChQ1TXOWVlZNG3aFF9fX3bv3s26desA6Nu3Lz///DOZmZkUFxfzxRdfVHpuQEAAERERtvu01mzduvWcn01tpiKvCUcNSM8BRpbdoJRyB94BrgAigQlKqUigFXC49GEWe4Zat2U249ZO5zcfb57s+RBvDntHTmoT4iIppfjyyy/54Ycf6NChA1FRUUyfPp2wsDCuvfZaYmNj6d69O0OHDuW1114jLCysRq8J8MorrzBq1CgGDBhQbnW4N998k3/+85/Exsayb98+22pqEydOZMOGDcTExDB37ly6du0KQLNmzWzTdj/66KNcfvnl3HTTTfTv35+YmBjGjRtXZUP73//+l8zMTFuhGD58OK+99lqVmUeOHElJSQndunWzLVcKEB4ezrPPPkv//v0ZOHBgtTPozps3jw8++IDu3bsTFRVlG0yvTmxsLO7u7nTv3r1OBqQdNmW3Uqod8D+tdXTp7f7As1rrEaW3p5c+NAU4qbX+n1Jqgdb6xmpebwowBaBNmza9znwjqI1FKx7m45QfeW3ELDq37Ffr5wvhjOrblN0xMTEsWbKkXF9+RXl5eTRq1AilFAsWLGD+/PnnbUwbGleasrslZ/cQwCgKfYGZwNtKqauApdU9WWs9C5gFxnoOFxJg7PB/MqowG28f2VsQwhkNHz6cmJiYcxYGgI0bNzJ16lS01jRp0sS2CJC4cE43IK21zgVuq8ljL3ZWVqWUFAYhnNj3339fo8cNGjTovH3yonbMPAnuCNC6zO1WpdtqzCGzsgrhgurDCo+i7lzI34OZxeF3oJNSKkIp5QXcCCypzQsopa5WSs3KysqyS0AhXJGPjw+ZmZlSIARgFIbMzEx8fGo3w4NDBqSVUvOBBCAYSAWe0Vp/oJS6EngTcAdma61fvJDXd+ga0kI4ueLiYlJSUmp8TL2o/3x8fGjVqhWenp7ltp9rQNphRyvZQ5kxh7uSk5PNjiOEEC7lXMXBpSfekzEHIYSwD5cuDkIIIezDpYuDDEgLIYR9uPSYwxlKqXTgEBAI1LZSBAMZZZ5b9jVqcp1abqvN/bXJe65857vt6KwVX08+W3M+2+p+r3y2Deezbau1DqnyVbXW9eYCzLqA52wo+9yyr1GT67XdZq+858p3vtuOziqfrXN8ttX9XvlsG95nW9XFpbuVqlDtdBu1eO7SKrad63ptt9Xm/po893y5a3K7Jlnks636dV35s63u98pnW/65DeGzraRedCtdDKXUBl3NoVzOyJXyulJWcK28rpQVXCuvK2UF++Wtb3sOF2LW+R/iVFwprytlBdfK60pZwbXyulJWsFPeBr/nIIQQojLZcxBCCFGJFAchhBCVSHEQQghRiRSHc1BKtVdKfaCUWmh2lppQSl2jlPqPUuozpdTlZuc5F6VUN6XUe0qphUqpe83Ocz5KKT+l1Aal1Cizs5yPUipBKbW69PNNMDvPuSil3JRSLyql3lJKTTI7z/kopQaVfq7/VUr9anaec1FKtVFKfaWUmq2Umlbb59fb4lD6gaQppZIqbB+plNqjlNp3vg9Ma31Aa32HfZPactVF3q+01ncB9wDjnTzrLq31PcANwEBnzlrqceBz+6Qsl6su8mogB/DBWH7XmbOOwVjoq9ieWUtz1cXf7erSv9v/AR85c1YgBliotb4d6FnrELU9G89VLsBgIA5IKrPNHdgPtAe8gK1AZOmH+L8Kl9Ayz1voYnn/AcQ5e1ZgNPANcJMzZwWGYyxGNRkY5ex/B4Bb6fOaA/OcPOs04G5H/D+r4/9jnwP+zpwVaAb8BKwEbqt1Bnv+Y5h9AdpV+HD7A9+VuT0dmF6D17F7caiLvIACXgUuc/asFV5rmTNnBV7EWJRqBfD1mcbXWfOWeZyXAxrci/1sbwZuKL3+mT2z1tVnC7QB/uPsWYFHgMGl12v9d+BBw9ISOFzmdgrQt7oHK6WaYTQMPZVS07XWL9s5X0W1ygvcD1wGBCqlOmqt37NnuApq+9kmAGMBb2C5PYNVoVZZtdZPAiilJgMZWmurXdNVVtvPdiwwAmgCvG3XZJXV9m92MfCWUmoQkGjPYNWobV6AO4AP7ZaoerXN+i3wrFLqJuBgbX9ZQysOtaK1zsTov3cJWuuZwEyzc9SE1noVsMrkGLWitZ5jdoaa0Fovxmh0nZ7WOg+jsXUZWutnzM5QE1rrJGDchT6/3g5IV+MI0LrM7Val25yVK+WVrPbjSnldKSu4Vl6HZm1oxeF3oJNSKkIp5YUxyLjE5Ezn4kp5Jav9uFJeV8oKrpXXsVntPahi1gWYDxzj7CFyd5RuvxLYizHq/6TZOV0xr2SVvK6W1dXyOkNWmXhPCCFEJQ2tW0kIIUQNSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIeykdGGWLaWX9Uop+f8mXIacIS2EnSilkjHm0z9mdhYhaku+yQhhP8uBbUqpN80OIkRtyXoOQtiBUmoAxsp84VrrErPzCFFbsucghH1cD+zVWpcoQ4DZgYSoDRlzEMIOlFJ9gA8ADeQD92mtN5qbSoiak+IghBCiEulWEkIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJX8PwHirA85xCEuAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"A = np.array([[ .2, .1, 1.0, 1.0, 0.0],\n",
|
|
" [ .1, 4.0, - 1.0, 1.0, - 1.0],\n",
|
|
" [ 1.0, -1.0, 60.0, .0, .0],\n",
|
|
" [ 1.0, 1.0, .0, 8.0, 4.0],\n",
|
|
" [ .0, -1.0, - 2.0, 4.0, 700.0]])\n",
|
|
"b = np.array( [ 1 , 2 , 3 , 4 , 5 ] )\n",
|
|
"x_exact = linalg.solve(A, b)\n",
|
|
"\n",
|
|
"eps_list = np.logspace(-8, -1, 8)\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"\n",
|
|
"for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\")]:\n",
|
|
" #diff_list = []\n",
|
|
" k_list = []\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" #diff_list.append(diff(x_exact, x))\n",
|
|
" k_list.append(k)\n",
|
|
" #ax.plot(eps_list, diff_list, label=alg_name)\n",
|
|
" ax.plot(eps_list, k_list, label=alg_name)\n",
|
|
" \n",
|
|
"ax.set_xscale(\"log\")\n",
|
|
"ax.invert_xaxis()\n",
|
|
"ax.set_yscale(\"log\")\n",
|
|
"ax.set_xlabel(\"$\\epsilon$\")\n",
|
|
"#ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n",
|
|
"ax.set_ylabel(\"$k$\")\n",
|
|
"ax.legend()\n",
|
|
"\n",
|
|
"fig.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"The condition number for A is K(A) = 12269.877667702964 >> 1,\n",
|
|
"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\n",
|
|
"This explains the great difference in order of required iterations k as observed above.\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"print(\"The condition number for A is K(A) =\", linalg.cond(A), \">> 1,\")\n",
|
|
"print(\"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\")\n",
|
|
"print(\"This explains the great difference in order of required iterations k as observed above.\")"
|
|
]
|
|
},
|
|
{
|
|
"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": 13,
|
|
"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": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"The number of iterations is brought down by a lot due to the conditioning.\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABDWElEQVR4nO3dd3xUVfr48c/JpEx6IQlJCL0mIQkldIGIougiKLIi2HBVVl1s+7Miu7q7Yv+u3VVcUVkFe8EVV9eCERdQUIRAgACChJ6EJKQnM+f3x50M6SSQzJ1JnvfrNa/M3HvnzpMhnOeecs9RWmuEEEKI2rzMDkAIIYT7keQghBCiAUkOQgghGpDkIIQQogFJDkIIIRqQ5CCEEKIBSQ5CCCEakOQghBCiAW+zA6hPKZUO/A3YAryptV51svdERkbqXr16tWtcQgjR0WzYsCFXax3V2D6XJAel1BJgKnBEaz241vYpwFOABfin1vphQAPFgBXIacn5e/Xqxfr169s8biGE6MiUUnub2ueqZqVXgSm1NyilLMBzwHlAIjBbKZUIfKu1Pg+4C/iLi+ITQghRi0uSg9Y6A8ivt3kksFNrvVtrXQm8CUzXWtsd+48Bfq6ITwghRF1m9jl0A/bVep0DjFJKzQDOBcKAZ5t6s1JqHjAPoEePHu0XpRBCdEJu1yGttX4feL8Fxy1WSh0ELvD19R1ef39VVRU5OTmUl5e3R5jCA1mtVuLj4/Hx8TE7FCHcnpnJYT/QvdbreMe2FtNafwx8nJaWdl39fTk5OQQHB9OrVy+UUqcXqfB4Wmvy8vLIycmhd+/eZocjhNsz8z6HH4D+SqneSilf4FJgRVudvLy8nC5dukhiEAAopejSpYvUJIVoIZckB6XUcmANMFAplaOUukZrXQ3MBz4DsoC3tdZbWnneC5RSiwsLC5vaf5qRi45E/h5Eh1BZCsf2Qs562LYSCvad/D2nwCXNSlrr2U1sXwmsPI3zNtms5C4OHz7Mbbfdxtq1awkPD8fX15c777yTiy66yGUxrF27lltuuYWKigoqKiqYNWsW999/f5PHr1+/nqVLl/L000832FdzT0lkZGST72/JMUIIB62hvBBKjhqP4iMNn9feVllc9/3TnoFhV7Z5WG7XId0aSqkLgAv69etndiiN0lpz4YUXctVVV7Fs2TIA9u7dy4oVbdZ61iJXXXUVb7/9NqmpqdhsNrZv397s8WlpaaSlpbkoOiE6ILsNSvOaLuzrP7dVNnISBQFdIDAKgqKg2/ATzwOjTzyP6NMuv4JHJwd3rzl89dVX+Pr6cv311zu39ezZk5tuuok9e/ZwxRVXUFJSAsCzzz7L2LFjWbVqFY8//jj//ve/AZg/fz5paWnMnTuXu+++mxUrVuDt7c0555zD448/zjvvvMNf/vIXLBYLoaGhZGRkNIjjyJEjxMbGAmCxWEhMTASgpKSEm266iczMTKqqqrj//vuZPn16nRjy8vKYPXs2+/fvZ8yYMdRec/z111/n6aefprKyklGjRvH8889jsVja7fsUwlTVFS0v7EvzwHnLVi1ePnUL+OjEhoV9zfOALmAxr4j26OTQUn/5eAtbDxS16TkT40K474KkZo/ZsmULw4YNa3RfdHQ0//3vf7FarWRnZzN79uxmpwDJy8vjgw8+YNu2bSilKCgoAOCvf/0rn332Gd26dXNuq++2225j4MCBpKenM2XKFK666iqsViuLFi1i0qRJLFmyhIKCAkaOHMnZZ59d571/+ctfOOOMM/jzn//MJ598wssvvwxAVlYWb731Ft999x0+Pj7ceOONvPHGG1x5ZdtXb4VoVxXHoXA/FOUYP48fbKTgPwoVjfdt4hN4olAP7w3dRzZe2AdFgTUMPKTvy6OTg7s3K9X3hz/8gdWrV+Pr68sXX3zB/Pnz2bhxIxaLhR07djT73tDQUKxWK9dccw1Tp05l6tSpAIwbN465c+dyySWXMGPGjEbf++c//5nLLruMzz//nGXLlrF8+XJWrVrF559/zooVK3j88ccBY4TXr7/+Wue9GRkZvP++cdvJb37zG8LDwwH48ssv2bBhAyNGjACgrKyM6OjoU/9yhGgPVWVQdAAKc6Bof90kUPO6sULfP9wo1IOiISa56cI+MAp8A13/e7mARyeHljYrnewKv70kJSXx3nvvOV8/99xz5ObmkpaWxhNPPEHXrl35+eefsdvtWK1WALy9vbHbT1RHa4Zeent78/333/Pll1/y7rvv8uyzz/LVV1/xwgsvsG7dOj755BOGDx/Ohg0buP322/npp5+Ii4tj5Uqjv79v377ccMMNXHfddURFRZGXl4fWmvfee4+BAwfWifvw4cMn/d201lx11VU89NBDp/09CXFKbFVGwd9ooe9IBqV5Dd8XEAmh3Yyr/F5nQEg3CI13/OwGQTHg7ev638fNeHRycHeTJk1iwYIF/OMf/+CGG24AoLS0FIDCwkLi4+Px8vLitddew2azAUafxNatW6moqKCsrIwvv/ySM844g+LiYkpLSzn//PMZN24cffoYnVC7du1i1KhRjBo1ik8//ZR9+/bxyiuv1Injk08+4fzzz0cpRXZ2NhaLhbCwMM4991yeeeYZnnnmGZRS/PTTTwwdOrTOeydMmMCyZctYuHAhn376KceOHQPgrLPOYvr06dx2221ER0eTn5/P8ePH6dmzZ7t+p6KTsNug+HDThX7hfmM/uu77/EKNgj60m9GBG9oNQhyvQxwPH6spv5Kn8ejk4O7NSkopPvzwQ2677TYeffRRoqKiCAwM5JFHHmHYsGFcfPHFLF26lClTphAYaFRNu3fvziWXXMLgwYPp3bu3s7A+fvw406dPp7y8HK01f//73wG44447yM7ORmvNWWedRWpqaoM4/vWvf3HbbbcREBCAt7c3b7zxBhaLhT/96U/ceuutpKSkYLfb6d27t7MjvMZ9993H7NmzSUpKYuzYsc55rBITE3nggQc455xzsNvt+Pj48Nxzz0lyECentXFF31hTT8224wfBXl33fT4BJ67u+519osCvnQD8gs35nTogVXv0iadKS0vT9Ttzs7KySEhIMCki4a7k78JFKksgNxtydxg/C/edKPiLDkB1vTvVLb4QElf3Kr/+Vb9/uMd05noKpdQGrXWj49Y9uuYghDCR1saIntwd9R6OZFBDeUFwrFHAx6bCwPMhtHutJBBv9AN4yarF7kSSgxCiebZqKNgLR7efKPxzHc/La4308QmEyP7QYwxEXgVRAyBygHGTlrcszeJpPDo5uHufgxAepaIY8rKNwr92IsjfVfcO3qCuRqE/eKbxsyYJBMfJ1X8H4tHJwd3vkBbC7dRpCtpeKxFkG53CNZQFInobhf6AcyByoPE8sj/4h5kWvnAdj04OQogm2Krh2J5G+gOaaArqNc74GTnASAQRvaUpqJOT5CCEJ6tpCjpaLwHk7QJ71YnjgmKMwn/wTIgaeCIRhHSTEUCiUZIc2tGiRYtYtmwZFosFLy8vXnzxRUaNGsWTTz7JvHnzCAgIMC22goICli1bxo033tjofovFQnJyMlVVVXh7e3PllVdy22234WVim/KDDz7IggULTPt801WVQ873sPsb2L/BSAJFtRZPdDYFDYQBUxz9AQOhSz9pChKt5tH3OdTqkL4uOzu7zj6zx7OvWbOGP/7xj6xatQo/Pz9yc3OprKwkLi7OLdY72LNnD1OnTiUzM7PR/UFBQRQXG/PGHzlyhDlz5jBu3Dj+8pe/uDLMJmM6VWb/XbSK3QYHNxrJ4Jdv4Ne1xv0BygIxgyFqkKMZqPaoIJn2QbRcc/c5oLX2+Mfw4cN1fVu3bm2wzZXee+89PXXq1Abbn3rqKe3j46MHDx6s09PTtdZaf/bZZ3r06NF66NCheubMmfr48eNaa63Xr1+vJ0yYoIcNG6bPOeccfeDAAa211hMnTtQ333yzTk1N1UlJSXrdunVaa62Li4v11VdfrUeMGKGHDBmiP/zwQ6211pmZmXrEiBE6NTVVJycn6x07duhZs2Zpq9WqU1NT9e23394gzsDAwDqvd+3apSMiIrTdbtfV1dX69ttv12lpaTo5OVm/8MILWmutDxw4oMePH++MKyMjQ2ut9aeffqqHDh2qU1JS9KRJk5qN9ZVXXtEXXXSRPvfcc3W/fv30HXfcobXW+q677tJeXl46NTVVz5kz55T/Xcz+u2iW3a714Syt176g9bLZWj/YXev7QozHc2O0XnmX1ttWal1WaHakooMA1usmylWPrjnUOOkd0p/eDYc2t+2HxiTDeQ83ubu4uJgzzjiD0tJSzj77bGbNmsXEiROBuiul5ebmMmPGDD799FPn1BoVFRXcc889TJw4kY8++oioqCjeeustPvvsM5YsWUJ6ejr9+/fnpZdeIiMjgxtvvJHMzEwWLFhAYmIil19+uXMK7p9++om7776b0aNHc9lll1FZWYnNZuPw4cMtrjnUCAsLY/v27Xz00UccOXKEhQsXUlFRwbhx43jnnXd4//33KS8v595778Vms1FaWkp5eTnDhg0jIyOD3r17k5+fT0RERJOxvvPOO/z1r3/lp59+ws/Pj4EDB7J69Wq6d+/eMWsOBfuMWsHub+CXDCg+ZGwP6wl9JkLvidB7gjE7qBBtTO6QNkFQUBAbNmzg22+/5euvv2bWrFk8/PDDzJ07t85xa9euZevWrYwbNw6AyspKxowZw/bt28nMzGTy5MkA2Gw254I9ALNnGyuvTpgwgaKiIgoKCpqcgnvMmDEsWrSInJwcZsyYQf/+/U/rd/v888/ZtGkT7777LmBMIpidnc2IESP43e9+R1VVFRdeeCFDhgxh1apVTJgwgd69ewMQERHhPEdT04WfddZZhIaGAsYcTnv37qV79+6nFbPbKMk1kkBNQjj2i7E9MMpIAr0nGkkhvJepYQrROZJDM1f47clisZCenk56ejrJycm89tprDZKD1prJkyezfPnyOts3b95MUlISa9asafTcqt4IE6VUk1NwJyQkMGrUKOfsrC+++KJzVteW2r17NxaLhejoaLTWPPPMM5x77rkNjsvIyOCTTz5h7ty5/PGPf3Su/1BfU7GuW7cOP78TQygtFgvV1dX13+45Koph7/9OJIPDjhqsb7AxXfSo3xsJITpBRg0JtyK3M7aT7du3U7uTfOPGjc4ZS4ODgzl+/DgAo0eP5rvvvmPnzp2AsXTnjh07GDhwIEePHnUmh6qqKrZs2eI831tvvQXA6tWrCQ0NJTQ01DkFd01T4U8//QQYBXufPn24+eabmT59Ops2baoTw8kcPXqU66+/nvnz56OU4txzz+Uf//gHVVXGUMkdO3ZQUlLC3r176dq1K9dddx3XXnstP/74I6NHjyYjI4NffjGukPPz8wGajLU5Pj4+zs90W9UVsGc1fLUIXj4HHukJy34L379kjBiatBCu+QLu2gNz3oTRN0DXREkMwu10jpqDCYqLi7npppsoKCjA29ubfv36sXjxYgDmzZvHlClTiIuL4+uvv+bVV19l9uzZVFRUAPDAAw8wYMAA3n33XW6++WYKCwuprq7m1ltvJSnJWLjIarUydOhQqqqqWLJkCUCTU3C//fbb/Otf/8LHx4eYmBgWLFhAREQE48aNY/DgwZx33nk89thjdeIvKytjyJAhzqGsV1xxBX/84x8BuPbaa9mzZw/Dhg1Da01UVBQffvghq1at4rHHHsPHx4egoCCWLl1KVFQUixcvZsaMGdjtdufyqC2ZLry+efPmkZKSwrBhw3jjjTfa9N/rlNltcGjTiRFFe9dAdZkx2VzcUBh7s9FM1H0U+PibHa0QLebRHdLuPJS1PaWnp/P444+Tltb4CDTRtNP+u9DamGril29g9yqjllBeYOyLGnSiz6DnOLm3QLi9DtshrWVuJeEKhfvrjig6fsDYHtoDEqZC73SjMzm4q5lRCtGmPDo5dFarVq0yO4SOrTQf9nx7oqkoz+gPIqBLvRFFvaWvQHRYkhyEAGNE0faVRkI4tBnQ4BtkNA+l/c4xoihRpqQWnYYkByG+ewr++2djqcr4kXDmAiMZdBsGFh+zoxPCFJIcROelNXz1AHz7OCTNgOnPgm+g2VEJ4RYkOYjOyW6H/9wN378Iw66EqU+Cl8XsqIRwG9KA2o4OHTrEpZdeSt++fRk+fDjnn38+O3bsOKVzXXvttWzdurVN4ysoKOD5559v9fuKi4u54YYb6Nu3L8OGDWP48OG89NJLpxXLq6++yvz58wF44YUXWLp06SmdZ8+ePSxbtqz5g2zV8NEfjMQwZj5c8LQkBiHqkeTQTrTWXHTRRaSnp7Nr1y42bNjAQw89xOHDh0/pfP/85z9JTExs0xhPNTlce+21hIeHk52dzY8//sh//vMf553PtZ3qtBfXX389V1555Sm996TJQWt4dy78vAzOvBfOeUBGHAnRCLdMDkqpQKXUeqXUVLNjOVVff/01Pj4+XH/99c5tqampjB8/Hq01d9xxB4MHDyY5Odk5FcaqVatIT09n5syZDBo0iMsuu8w5vUR6ejo1M88GBQU5z/nuu+8652vatWsXo0ePJjk5mYULFzqPKy4u5qyzzmLYsGEkJyfz0UcfAXD33Xeza9cuhgwZwh133AHAY489xogRI0hJSeG+++5r8Hvt2rWL77//ngceeMC58E9UVBR33XWX83cYP34806ZNcyazCy+8kOHDh5OUlOS8SxzglVdeYcCAAYwcOZLvvvvOuf3+++93Tsi3a9cupkyZwvDhwxk/fjzbtm0DYO7cudx8882MHTuWPn36OCcBvPvuu/n2228ZMmQITzzxRN3g7TYoOQpZH8OUh2HinZIYhGiCS/oclFJLgKnAEa314FrbpwBPARbgn1rrmhny7gLebqvPf+T7R9iWv62tTgfAoIhB3DXyrib3Z2ZmMnz48Eb3vf/++2zcuJGff/6Z3NxcRowYwYQJEwBjjqEtW7YQFxfHuHHj+O677zjjjDNaFNMtt9zCLbfcwuzZs3nhhRec261WKx988AEhISHk5uYyevRopk2bxsMPP0xmZiYbN24EjJlSs7Oz+f7779FaM23aNDIyMpyxAWzZsoXU1NRmV4T78ccfyczMdM7EumTJEiIiIigrK2PEiBFcfPHFVFZWct9997FhwwZCQ0M588wzGTp0aINzzZs3jxdeeIH+/fuzbt06brzxRr766isADh48yOrVq9m2bRvTpk1j5syZPPzwwzz++OMNp+KwV0PebmPuo+nPwdDLW/SdCtFZuapD+lXgWcDZkKyUsgDPAZOBHOAHpdQKoBuwFbC6KDaXW716NbNnz8ZisdC1a1cmTpzIDz/8QEhICCNHjiQ+Ph6AIUOGsGfPnhYnhzVr1vDhhx8CMGfOHG6//XbAaOJasGABGRkZeHl5sX///kabtz7//HM+//xzZyFdXFxMdnZ2neRQ36JFi3jnnXc4cuQIBw4Ydw6PHDnSmRgAnn76aT744AMA9u3bR3Z2NocOHSI9PZ2oqCgAZs2a1aA/pri4mP/973/89re/dW6rmX8KjBqJl5cXiYmJzTfX2aqMNZWryyGwC6Sc0/SxQgjARclBa52hlOpVb/NIYKfWejeAUupNYDoQBAQCiUCZUmql1tpe/5xKqXnAPIAePXo0+/nNXeG3l6SkJGdTR2u0ZLrq2tN1l5eXn/Scb7zxBkePHmXDhg34+PjQq1evRt+nteaee+7h97//fZPnSkxM5Oeff8Zut+Pl5cW9997LvffeW6epKzDwxHDQVatW8cUXX7BmzRoCAgJIT09vUcwAdrudsLAwZ82mvtrfVZNzhFVXGnc426qMZTQL9zd+nBCiDjP7HLoB+2q9zgG6aa3v1VrfCiwDXmosMQBorRdrrdO01mk1V5/uZNKkSVRUVNRpY9+0aRPffvst48eP56233sJms3H06FEyMjIYOXJki8/dtWtXsrKysNvtzityMKb/fu+99wB48803ndsLCwuJjo7Gx8eHr7/+mr179wI0mLb73HPPZcmSJc7V1vbv38+RI0fqfHa/fv1IS0tj4cKF2Gw2wEhQTRXOhYWFhIeHExAQwLZt21i7di0Ao0aN4ptvviEvL4+qqireeeedBu8NCQmhd+/ezn1aa37++edmv5s6v1N1OeRlG01KXfqCNaTZ9wohTnDLDmkArfWrWutm53BWSl2glFpcWFjoqrBaTCnFBx98wBdffEHfvn1JSkrinnvuISYmhosuuoiUlBRSU1OZNGkSjz76KDExMS06J8DDDz/M1KlTGTt2bJ3V4Z588kn+/ve/k5KSws6dO52rqV122WWsX7+e5ORkli5dyqBBgwDo0qWLc9ruO+64g3POOYc5c+YwZswYkpOTmTlzZqNrPvzzn/8kLy/PmSgmT57Mo48+2mjMU6ZMobq6moSEBOdypQCxsbHcf//9jBkzhnHjxjU5U+obb7zByy+/TGpqKklJSc7O9KakpKRgsVhITUnhiUV/Bm2HLv3AL6jZ9wkh6nLZlN2OZqV/13RIK6XGAPdrrc91vL4HQGv9UGvPfdI1pDuA5ORkVqxYUactv77S0lL8/f1RSvHmm2+yfPnykxamHVJlidHHoLyMxOBzovuqo/1dCHE63HXK7h+A/kqp3sB+4FJgTmtOUGs9h3YIz31MnjyZ5OTkZhMDwIYNG5g/fz5aa8LCwpyLAHUqFcchfzd4eRuJwdvv5O8RQjTgqqGsy4F0IFIplQPcp7V+WSk1H/gMYyjrEq31lmZO00BnWc/hv//9b4uOGz9+/Enb5Du08kLI/8VICF36GhPpCSFOiatGK81uYvtKYKUrYhAdXGk+FPxqLMUZ0RcsMm2YEKfDbTukW8KdO6SFC5XkQsFeY0bVLv0kMQjRBjw6OWitP9Zaz6sZlSM6oeLDULgP/EKMGoNMoCdEm/Do5CA1h05Mayg6YDysYRDRW1ZpE6INefT/JnevOciU3S3Xqim7tYai/UatIaALhPcyhq3Swim7hRAn5dHJwZ3JlN3tNGW31kbHc8lRCIyG0O51ZlaV5CBE25Dk0E5kyu52mrJ73RdQls/cOx7m5oUPM3bcuJZP2S2EaDGPHtbR0pvgDj34IBVZbTtlt1/CIGIWLGhyv0zZ3cZTdvftw7ovPuDG/7eArz77BHz8OXjoUMun7BZCtIpH1xzcvc+hKU1N2Q04p+z28vJyTtndUmvWrHFObz1nzombzWum7E5JSeHss89u0ZTdw4YNY9u2bWRnZzf7mYsWLWLIkCHExcU5tzU2ZXdqaiqjR492Ttm9bt0655Tdvr6+zJo1q8G5T0zZPZMhKYP5/W33cjC3EIKigVZM2S2EaDWPrjm0VHNX+O1Fpuw2tMmU3f99C6rKIbwn+Ic797doym4hxCnx6JqDOw9llSm7T3z2KU/ZHWild3wM77z/MUT0QVvDWjdltxDilHl0cnDnZiWZsttwylN226shN5s3nl3Ey+99Tuqo8a2bsjs1VTqkhTgNLpuyuz3JlN2GDjNld1WZMeW2thsT6PkGnvw9LdTR/i6EOB3uOmW3aKFONWV37bUYIvsbE+kJIVxOkoMH6DRTdstaDEK4DUkOwj3IWgxCuBWP7pA+2WiljtCf0imUHTMSg48VuvRvt8Qgfw9CtJxHJ4fmRitZrVby8vKkQHB3JblwbI9jLYb+7bYWg9aavLw8rFbryQ8WQnTcZqX4+HhycnI4evSo2aGIplQUQVkBePtDoC8cPbUZa1vKarUSHx/frp8hREfRYZODj4/PSUf3CJNoDV8/CBmPQuKFMOMl8JY+BiHcSYdNDsJN2e3w2T2w7gUYegVc8JSs3iaEG5LkIFzHVg0f3wwb34DRf4BzF9VZi0EI4T4kOQjXqK6E96+FrR9B+j0w8S5JDEK4MY9ODi1dz0GYrLIU3r4Cdn4B5z4IY/5gdkRCiJPosENZhZsoL4TXL4adX8K0ZyQxCOEhPLrmINxcSR68fhEc3gIzl8DgGWZHJIRoIUkOon0UHYSl06FgL1y6HAacY3ZEQohWkOQg2sfK26EwBy5/D3q1bA1sIYT78Og+B+Gm9v8I2/4NZ9wqiUEIDyXJQbS9rxeBfwSMut7sSIQQp0iSg2hbe9cYQ1bPuBWsIWZHI4Q4RW6XHJRSCUqpF5RS7yqlbjA7HtEKWsNXD0BgNIy4zuxohBCnwSXJQSm1RCl1RCmVWW/7FKXUdqXUTqXU3QBa6yyt9fXAJcA4V8Qn2sgv38De1TDhdvANMDsaIcRpcFXN4VVgSu0NSikL8BxwHpAIzFZKJTr2TQM+AVa6KD5xumpqDSHxMHyu2dEIIU6TS5KD1joDyK+3eSSwU2u9W2tdCbwJTHccv0JrfR5wmSviE20g+3PI+cGoNcjaz0J4PDPvc+gG7Kv1OgcYpZRKB2YAfjRTc1BKzQPmAfTo0aPdghQtUFNrCO8FQy83OxohRBtwu5vgtNargFUtOG4xsBggLS1N1gI1U9bHcGgTXPgCWHzMjkYI0QbMHK20H+he63W8Y1uLKaUuUEotLiwsbNPARCvYbcaqbpEDIOUSs6MRQrQRM5PDD0B/pVRvpZQvcCmwojUnkFlZ3UDm+3A0y1ijQVZ0E6LDcNVQ1uXAGmCgUipHKXWN1roamA98BmQBb2utt7TyvFJzMJOtGlY9CF0HG2tBCyE6DJf0OWitZzexfSWnMVxVa/0x8HFaWprccWWGn5dD/m5j1lUvt7ufUghxGjz6f7TUHExUXQHfPAJxw2DgeWZHI4RoYx6dHKTPwUQ/LoXCfTBpoawFLUQH5NHJQZikqgwyHoceY6HvJLOjEUK0A49ODtKsZJIfXobiQzDpXqk1CNFBeXRykGYlE1QUw+onoE+6LOQjRAfm0clBmOD7F6E0F85caHYkQoh25NHJQZqVXKysAL57CgZMge4jzI5GCNGOPDo5SLOSi619HsoL4cwFZkcihGhnHp0chAuV5MGa5yFxOsSmmh2NEKKdSXIQLfO/p6CyGNKl1iBEZ+DRyUH6HFzk+GFYt9iYdTV6kNnRCCFcwKOTg/Q5uMjqv4OtEibeZXYkQggX8ejkIFygMAfWL4Ehc6BLX7OjEUK4iCQH0byMx41lQCfeaXYkQggXkuQgmpb/C/z0Lxg+F8JknW4hOhOPTg7SId3OvnkUvLxh/P8zOxIhhIu1ODkopR5oZJup60JKh3Q7OroDNr0JI66FkFizoxFCuFhrag7dlFLOFd2UUtHAF20fknALqx4Cb3844zazIxFCmKA1y4T+HvhMKbUL0MArgIxt7IgOZcKW92H87RAYaXY0QggTnDQ5KKWWAj8CPwF/AJYB1cCFWuud7RueMMXXD4JfKIydb3YkQgiTtKRZ6VVAAVcDrwO9gGPA5Uqpme0WmTDH/g2w/RMYexP4h5sdjRDCJCetOWitvwK+qnmtlPIGEoBUYBTwbrtFJ1zvq0XgHwGjrzc7EiGEiVrT5wCA1roa2Ox4vN7mEbWCUuoC4IJ+/fqZGUbHsXcN7PoSJv8V/ILNjkYIYSKPvs9BhrK2Ia3hqwcgMBpGXGd2NEIIk3l0chBt6JdvYO9qmHA7+AaYHY0QwmSSHMSJWkNIvDFVhhCi05PkICD7c8j5ASbeAd5+ZkcjhHADkhw6O7vdqDWE94Ihl5kdjRDCTbR6tJLoYLZ9DIc2wUUvgsXH7GiEEG5Cag6dmd1m3A0dOQCSf2t2NEIINyI1h84s8z04ug1++yp4mTrBrhDCzbhlclBKXQj8BggBXtZaf25uRB2QrdqYebVrMiRMNzsaIYSbcVmzklJqiVLqiFIqs972KUqp7UqpnUqpuwG01h9qra8DrgdmuSrGTuXn5ZC/G85cAF7SuiiEqMuVpcKrwJTaGxyLBT0HnAckArOVUom1Dlno2C/aUnUFfPMIxA2DgeeZHY0Qwg25LDlorTOA/HqbRwI7tda7tdaVwJvAdGV4BPhUa/2jq2LsNH5cCoX7YNJCUMrsaIQQbsjs9oRuwL5ar3Mc224CzgZmKqUanR5UKTVPKbVeKbX+6NGj7R9pR1FVBhmPQ4+x0HeS2dEIIdyUW3ZIa62fBp4+yTGLgcUAaWlp2hVxdQg/vAzFh2DmEqk1CCGaZHbNYT/QvdbreMe2FlFKXaCUWlxYWNjmgXVIFcWw+u/Q50zoNc7saIQQbszs5PAD0F8p1Vsp5QtcCqxo6Ztlyu5WWvcClOYZfQ1CCNEMVw5lXQ6sAQYqpXKUUtc4Fg6aD3wGZAFva623tOKcUnNoqbIC+N/TMOA8iE8zOxohhJtTWnt+c31aWppev3692WG4t68fNIav/v5biE0xOxohhBtQSm3QWjd6tWh2s9JpkZpDC5XkwZrnIfFCSQxCiBbx6OQgfQ4t9L+noLIY0u8xOxIhhIdwy6Gsog0dPwzrFkPKJRA9yOxohBCNsNs1JZXVlFTYKK6opriimpJaP43nNue22vvnTejD+P5RbR6TRycHpdQFwAX9+vUzOxT3tfrvYKuEiXeZHYkQHUpFtY3i8hMFekll3cL8uGNf/e3G87pJoLTS1qLP9FIQ5OdNkJ83gY5HZbW9XX4/j04OWuuPgY/T0tKuMzsWt1SYA+uXwNDLoEtfs6MRwlTVNrtRKFeeKKidBXtFNaWVda/OSxwFvlHI100CJRXVVNlaNpjHz9uLYKujMPc1CvbIIF96dgkwtvsa+04U+JY6CaD2T6uPF8pFN696dHIQJ5HxmPFzwp3mxiHEKaiy2Z0FdWmlrc7Vd8MrcpuzMC+usFFaU8DXaqpp6RW2l8IoxK0nrs6D/CxEBvkZhba1VqHta6lbgFtrPfc1Cnpvi2d27Xp0cpBmpWbk74afXoe030FY95MfL0Qb0FpTUmmjsKyKwtIqCsoqKSqr20ZeWln3ir2k0lanyaUmEbSqMPeruQK3OAvn+IAAgvzqFt4BvpY6V+IB9Qr3QD8L/j4Wl12dn44jpUfYfHQzSZFJxATGtPn5PTo5SLNSM755FLy8Yfz/MzsS4WG01pRVGQV8QWmVUdA7HkWOn41tL3D8rLY339xSU5jXbjcP9LUQERjg2GZxNrXUXLUH+DZsdqnZ5sqmFrOUV5eTlZ/FpqObjEfuJg6VHAJg4aiFzBrU9sveeHRyEE04ugM2vQVj/gDBbX9FITxDuaOAb6wwLyytrPu6VuFeWFbVbHu6UhBi9SEswIdQf+PRLdyfUH8fwvxPbKt5hPj7ONvcg/y88fPu+IX56dBas7doL5tyNzmTQfaxbKp1NQDdgroxNGooyYnJpESlMCiifUYhSnLoiFY9BD4BMO5WsyMRbaC8ykZ+SWWdxzFH4V5QeqJAr13IF5ZVnbRZJsTqTWitAj4m1Eqov2+Dwr12Egjx9yHYzxsvLync20phRSGbczc7awSbj26mqLIIgECfQAZHDubqwVeTHJlMclQykf6RLonLo5OD9Dk04lAmbHkfxt8Oga75IxItZ7drisqrGhT2eSWVHKt5XVp3X3PDHIP9vAmpVZD3iw6qU5DXL9xrHsFWHyxSwLtclb2KHcd2sOmokQQ25W5ib9FeALyUF33D+jK552RSolJIiUyhd2hvLF4WU2KVuZU6muVzYM9quPVn8A83O5oOr6LaxrGSWoV9aSX5xRWNFvLGFX8Vtiba5P19LEQE+rboER7gS4jV22NHwnQGWmsOlRxyNg9tzt3M1rytVNgqAIj0jyQlMoXkqGRSIlNIikwi0CfQpTE2N7eSR9ccRD37N8D2T+DMhZIYToHWmuMV1eQX1xTyjRfwtR/FFdWNnkspCPP3cRbmvSMDGd4z3PHaj4hAH+NngC8RQb5EBPji72vOFaJoG6VVpWzJ2+LsJ9icu5mjZcYqlX4WPxIiEpg1cBbJUcmkRqYSExjj1n0vkhw6kq8WgX8EjG50ZdVOy2bX5JVUcKSogiPHyzlcZDw/fLycI0XlHDleweGicvJLKpvsiPX19qJLrSv3nl0CCA/wpUugL+GBvnX2RQT6EhbgK802HZhd29ldsNvoK3DUDHYW7MSujX6eniE9GRU7ytk8NCB8AD4WH5Ojbh1JDh3F3v/Bri9h8t/AL9jsaFzCbtfkl1ZyuKjcKOxrFfSHiyo46kgER4srGm3KiQj0JTrYj64hVgZ2DaZLkF+ThX2Ar2eMfRftI68sr06n8ZbcLRRXFQMQ7BtMSmQKZ/U4y+g0jkwmzBpmbsBtwKOTg3RIO2gNXz0AQV1hxLVmR3Pa7HbNsdJKDte+unc8P1xUwZHjFRwpKufo8YpGx9SHB/jQNcRKdIiV/l2D6RpiJIDoYCvRjudRQX74ekt7vWio0lbJtvxtzkSw6egm9hcbqxdblIUB4QP4TZ/fkBKVQnJkMj1DeuKlOt7fkkcnB7kJzmH3Ktj7HZz3GPgGmB1Nk7TWHCutclzZlzsL+cO1rvprmnkaK/TDAnzo6ijg+0VFOgv9riF+RAXX/PTDz1va7s2WW5bL2oNrWXtgLRsOb6DcVm52SC1WWFFIlb0KgJjAGJIjk7l04KWkRKWQ0CUBf29/kyN0DY9ODoITtYaQeBh+lWlhVFTbOFRYzsHCcg4WlnGwsJxDheW12vaN9v7G2vRD/X2cBX2fqC5Gge9o7okO8SM62EpUsB9WHyn03VVpVSnrD683EsLBtWQfywYgxDeEkTEjPaqZJcQ3xDmKKDog2uxwTCPJwdNlfw7718MFT4O3X7t8REW1jcOFFc5Cv3YCOFhYxsGCcvJKKhu8L8Tq7biytzKqdyDRjqt8o4nH0bwjhb5HqrZXsyVvC2sPrGXNwTX8fPRnqu3V+Hr5MrTrUG4Zdgtj4sYwKHyQaeP0xemR5ODJ7Haj1hDeG4bMOaVTVFbbOVxUr8AvqJ0EysktrmjwvhCrN7Gh/sSGWUnuFkZsqNXxMLbFhFgJ9JM/r45Ca82eoj3OpqIfDv3A8arjACREJHBF4hWMjh3NsOhhWL2tJkcr2oL87/Vk2z6GQ5vgohehkWFyVbYTBf+BgrIGzT41BX/9+yCDrd7Ogj4pLsRZ4DsL/1Ap+DuD3LJc1h1c52wqqpnoLS4wjnN6ncPouNGMihlFuFXuqemI5H+4h6qqqoIvHqA6tC9f6HEczNjFgYJyRwIwCv+jjRX8ft7EhFqJDfMnISaE2DArcaH+xIRaiQuzEhPqT5AU/J1SaVUpPx750dlUtOPYDsBogx8VO4rrkq9jTOwY4oPjZVhvJ+DRpUBHH8qqtebI8Qq2Hihi68Eisg4Wse9YGYcKyxhb8iVP+OzglsqbWfnmJgACfS3EhhlX9gNjgp1X+bFh/sSFWokJtRJs9awbcUT7sdltRr+Bo2aw8chGquxV+Hj5MCx6GLcMu4XRsaNJiEiQfoNOSOZWchPVNju7c0vYesBIAlsPFrH1QFGdjt7uEf706hJIXLA3d+++CnwC+Pk3K4gNCyQ2zEqIFPyiGVprfj3+K2sOrGHtwbV8f+h7jlca/QaDIgYxOnY0Y2LHMLTr0E4zXLOzk7mV3Mzx8iq2HTpuJAFHrWDboePOKZZ9LV4MiAnirIRoEmNDSIwLZVBs8InC/8elsHUfXPQm6QNlvQbRtPzyfGe/wZoDazhYchCA2MBYJveczOjY0YyKHUWENcLkSIW7keTQjrTWHCwsr5MEth4sYm9eqfOY8AAfEuNCuGpMTxLjQkiMDaVPVCA+Tc22WV1hrPLWbTgMmOKi30R4irLqMn48/KOzqWhb/jbAmOJhVMworhl8DaPjRtMjuIf0G4hmSXJoI1U2OzuPFDdIBAWlVc5jenUJICkuhN8OjycxLoSE2BBiQqyt+0/641Io3AfTnjam/hSdms1uIys/y9lU9NORn6iyV+Ht5c3Q6KHcNPQmxsSOIbFLovQbiFaR5HAKCsuq2FarX2DrwSKyDxdTaTOahfy8vRgUE8x5g2NIjDWSwKDYkNMfBVRVBhmPQ89x0OfMNvhNRJW9ikPFh8gpznFOmeAJDpUcYs2BNaw7tM7ZbzAwfCBzBs1hdJxxv0GAj/tOpSKap7WG6mp0ZSX2ykp0RYXxqKzEXlGJrjSe64oK/AYMwCc2ts1jkOTQDK01OcfK6nQQbz1YRM6xMucxXQJ9SYwL4eozehn9A7Eh9I4MbJ9FWH54GYoPwcwlUmtohdKqUvYd30fO8Rz2Hd9X53Gw5CA23fRKa+6sa0BXzupxFmNixzAydqTLlo/s6LTWUFV1olB2/KxfKNsrKtCVVca22q8rKpzH2StqFexVtV7XnKOy0vm8/mvszS/zWiN20SLCLp7R5t+DJAeHymo72UeOn2gScowaKio3FnNRCnpHBpLaPYzZI3uQGBdCUmwIUcF+rmm7zdsF3/4f9J0Evca1/+d5EK01+eX5zgK/fhLIK8+rc3yoXyjdg7qTHJnMeb3Po3twd+KD4z1qhE6Ibwjdg7tLv0ETtN2OvaQEW2EhtsJC7EVFjudFjtcnntuKirAVFWIvMJ7bS0pocINQaymF8vND+fqi/Hzx8q157ud87RUUhMX52hfl6+fcr3x98fLzM7bVnKPR1774dO/eNl9aPZ06OXyz4ygfbdzP1gNF7DxS7JwJ1N/HwsCYYKamxjlGC4UwKCaYAF+Tvq6CX2HpdCNDTXnEnBhMVm2v5lDJoSYTQGn1iU5+haJrYFe6B3dnYveJzsK/e3B3ugd3J8Q3xMTfRLSU1hpdVmYU3oVF2AoL6hbyRY6Cv7DIcUzhiUL++PFmr7yVry9eoSFYQkOxhITi0zUGS/8BxragIJSftUEhbBTejbx2FvrGdi9fX/Dx8fjE7XbJQSnVB7gXCNVaz2zPz8o+fJxvs3NJjA3hzEHRzkTQq0ug+6ziVXQQXpsGFUVw1ccQNcDsiNpNWXVZg0K/5vWB4gNU6xNLcvp4+TgL/LSYNGfBHx8cT7egbvhZ2mcSQtF6urLyROHtKNTthXWv3O1FhdgcV+7ObYWF6Kpm+oG8vLCEGAW8V2golrAwfHv0wBIaYrwOCTUK/9AQLCGObY6Hl1XmfzoZl9wEp5RaAkwFjmitB9faPgV4CrAA/9RaP1xr37stTQ6nehOcza7dJwk0piQXXjkfivbDFR9C9xFmR3RatNYUVhQ2aPevSQJHyo7UOT7YN9hZ6Nd/RAdEd8gFVjyBvbSU6vxj2I7lY8vPdzw3Xlfn52PLP2ZsP2b8tBcXN3s+r6AgRwF/4kreEhKCJSwUr5Ba2xz7vUJCjX2BgR5/dW42d7gJ7lXgWWBpraAswHPAZCAH+EEptUJrvdVFMbl3Yig7BksvNJqULn/XoxJDpa2S7fnb2XFsR4MEUDOTZ41o/2jig+MZ221sgwQQ6hdq0m/QeWitsR8/XquQr1/AO54fO+Z8rsubWLjHxwfv8HAsERF4R4TjHx+PJTwcS0Q4lrCwxq/kg4NR3m7XgCFwUXLQWmcopXrV2zwS2Km13g2glHoTmA60KDkopeYB8wB69OjRdsG6g/IieP1iyN0Os5dDrzPMjqhJNruNXYW72JK7hczcTDLzMtlxbAfVdqMJyNvLm25B3YgPjic1KrVO4R8fHC/TO7cxbbNhKyho+or+mGN7TcF/rACqqxs9l/L3dxb2li4R+PXrZzyPCMc7IsIo+MMdzyMi8AoKkiv5DsTMlN0N2FfrdQ4wSinVBVgEDFVK3aO1fqixN2utFwOLwWhWau9gXaayFJbNggMbYda/oN/ZZkfkpLUm53gOmXmZRiLIzSQrP4uyamNob5BPEEmRSVyVeBWDIwczKGIQsYGxcvNVK+iqKuwlJdhLS41H7ec1r0tKsR0vwnasoO7VfX4+tsLCJkfaeAUHGwV7eAQ+8fFYU5LxDo9wXulbIiKwhEfgHR5mFPb+njN6S7Q9t6vPaa3zgOtbcmyHm5W1qhzenAP71sLF/4RBvzE1nKOlR521gS25W8jMy6SwohAAP4sfgyIGMaP/DJK6JDE4cnCHXWi9Kbq6umEBXlKKvbTE8bN+AV9S95jSUnRpKbaSErTj+GY7YGvz8jKaasLD8Q4Pd1zV11zR17q6d1zhe4eFoXx92/cLER2KmclhP1B7gG68Y1uLaa0/Bj5OS0u7ri0DM4WtCt6ZC7u/hunPw+CLXfrxRZVFbMndwpY8o3loc+5mjpQaHcQWZaFfWD/O7nE2SZFJDO4ymH7h/fDx8vxZYHVVFRW7dlGxfTvVx47VLbBrF/iNXMHryoZLozZFWa14BQSceAQGYgkKxiu6q/O1V2Dd/V4BAah6r70CHMf5+6MsUiMT7cfM5PAD0F8p1RsjKVwKtGqtyw5Tc7Db4P3rYMencP7jMPSydv24suoytudvZ3PuZjJzM9mSt4W9RXud+3uG9CStaxqDIwc7m4c86QaxpthLSijfvp3yrCzKs7Ko2JpFRXZ2g6t15efXoCD3CgzEOzqq8YK7fuFdv7D395dOV+FxXDWUdTmQDkQCh4H7tNYvK6XOB57EGMq6RGu96FTO79HrOdjt8NEf4OdlMPlvMO7mNj19lb2KXQW72Jy72dlpvLNgp3PKiOiAaAZ3MZJAUmQSSV2SOsQooer8fMq3ZlGetZWKrCzKt2ZRuXevsz3eEh6ONSEBa2ICfgkJWBMS8I6MNApyH8+vEQnREqYPZdVaz25i+0pgpSticEtaw8rbjcSQvuC0E4Nd29lbtNdZG8jMzWRb/jYqbBWAMeXC4MjBTIif4KwVRAdEt8VvYhqtNVX791O+dSsV27Y5EkIW1YcPO4/xiYvDLzGBkAumYk1IxJqYgHfXrjKyRohmeHRd16OblbSGzxfC+pdh3C0w8c5Wvl1zuPSwc9RQZl4mW3O3Ou8j8Pf2JyEigUsGXsLgLoNJjkz2+LV/dXU1Fbt3O2sC5VlZlG/bhr2oyDjAywu/vn0IGDXSSAIJCVgHDcQSFmZq3EJ4Ilkm1CxfPwjfPAIj58F5j7ZoltWNRzay9uBaZ0KomVDOW3nTP7y/szYwOHIwfUL74O3lubnfXlZGxfbtlNeqDVTs2GHMVonRL+A3cKCRAGqahwYMkGkRhGgF05uV2ovH1hxWP2EkhqGXGxPpnSQxlFeX89gPj/H2jrdRKHqF9mJs3FiSIpNIjkxmYMRAj55LyFZQYNQCamoDWVlU/vKLc+I0r5AQrAkJhM+ejTXRSAa+vXtLJ68Q7UhqDq627kX49E4YPBNmLIaT3CC289hO7si4g50FO5mbNJd5KfMI9g12UbBtS2tN9cGDjgSwzfFzK9UHDjqP8Y6JqVsbGJSAT7c4j24OE8Jdddiaw+myl5RgLys7+YFtZfM78NkC6HcujF8E+ceaPFRrzce7P+bZn54jwCeAF0c+wsjYkVBUQTUVrov5NNgKC2vVBrZSkbUNW0GBsVMpfHv1ImDIUKxz5pwYMRQhC90L4Q46dc3h6PPPk/v0M+0QkahN+fjg178/fo4mIWtCItaBA/AKDDQ7NCE6tQ5bczjdPoeg8RPwDg9v26Aac2gzbHgNIvrAyOvA0vQ4+l+L9vHhzg8orizmzB5nMip2NF4e2qTiFRCA36BB+PXpI/cOCOFhOnXNwSWy/wvLZ0PcULjiffBrvL/AZrfx0uaX+MfP/yAuMI7HJj7G4MjBjR4rhBBtocPWHNzeLxnw1uXQNREue6fJxHCo5BD3fHsP6w+v5/ze5/On0X8iyDfIxcEKIcQJkhzay6/rYNmlEN4bLv8A/MMaPezrX7/mT//7E5W2Sh4Y9wDT+k6TkTlCCNN5dHJw2/scDvwEb8yE4Bi48kMI7NLgkApbBf+3/v9Yvm05CREJPDrhUXqF9nJ5qEII0RiPnnxfa/2x1npeaKgbTRR3eCv86yKwhsFVK4wEUc/ugt3M+WQOy7ct54rEK3j9/NclMQgh3IpH1xzcTu5OWDodvK1w1UcQGl9nt9aaD3Z+wMPfP4zVYuW5s55jQvwEk4IVQoimSXJoK8f2wtJpoO1w5SfGsNVaiiqL+Nuav/GfPf9hVMwoHhz/oMfPiCqE6LgkObSFogPw2gVQWQJz/w1RA+rs3nhkI3d/ezeHSg5xy7BbuDrpallXWQjh1jw6ObhFh3TxEXhtGpTmG01JMcnOXTa7jVe2vMKzPz1LTGAMr533GqlRqebFKoQQLSQd0qejNB+WXghF+437GLoNd+46UnqE3//39zz141NM7jmZdy54RxKDEMJjeHTNwVTlhfD6DMjbCXPegp5jnLsycjJYuHoh5bZy/jr2r1zY70K5d0EI4VEkOZyKyhJ44xJjzqRZb0DfM43Ntkqe2PAEr2e9zoDwATw24TH6hPU5ycmEEML9SHJorapyY66knO9h5hIYOAWAPYV7uDPjTrLys5gzaA5/TPujRy/AI4To3CQ5tEZ1Jbx9pTFn0kUvQNJFaK35aNdHPLjuQfwsfjwz6RnSu6ebHakQQpwWSQ4tZauG966B7M9g6hOQeinFlcX8be3fWPnLSkbEjOChMx6ia2BXsyMVQojT5tHJwWVDWe12+OhGyFoB5z4Eab9j89HN3JlxJwdLDjJ/yHyuTb5W7l0QQnQYMpT15B8Cn9wGm96CSQuxj76eJZlLuPLTK7FpG69MeYXfp/5eEoMQokPx6JpDu9Ma/nMPbHgVxv8/ckdczYL/Xs+ag2uY3HMy9425j1A/N5r0Twgh2ogkh+Z89TdY9w8YdQPfDZzEghUXU1JVwp/H/JmZ/WfKvQtCiA5LkkNTMh6Db/+PqmFX8nRkJK9+eQP9wvrx8jkv0y/czdaPEEKINibJoTFrnoevHuDXwdO40+soW7auYtbAWdyedjtWb6vZ0QkhRLuT5FDf+iXw2T18POAMHqjYgXeVN0+mP8lZPc8yOzIhhHAZSQ61/fwmJZ/8Px7sk8yKql8ZFj2Mh8c/TGxQrNmRCSGES0lyqLHlA7asvJk7e/Uih+PcmHoj16Vch7eXfEVCiM7H7Uo+pVQg8DxQCazSWr/R3p9p37aSf31+C0/GdqWLfwQvT3iEtJi09v5YIYRwWy65CU4ptUQpdUQplVlv+xSl1Hal1E6l1N2OzTOAd7XW1wHT2ju2vKwV/GHVrTweEcqE+PG8N/19SQxCiE7PVXdIvwpMqb1BKWUBngPOAxKB2UqpRCAe2Oc4zNaeQa3duISZa+7he6sf9w69lSfPek5uahNCCFzUrKS1zlBK9aq3eSSwU2u9G0Ap9SYwHcjBSBAbaSZ5KaXmAfMAevTocUpx7T+ymVAsvDh5MQO6jT6lcwghREdkZp9DN07UEMBICqOAp4FnlVK/AT5u6s1a68XAYoC0tDR9KgHMmPx3plYU4WeV2oIQQtTmdh3SWusS4OqWHHu6s7IqpSQxCCFEI8yclXU/0L3W63jHthZzyaysQgjRCZmZHH4A+iuleiulfIFLgRWtOYFS6gKl1OLCwsJ2CVAIITorVw1lXQ6sAQYqpXKUUtdorauB+cBnQBbwttZ6S2vOKzUHIYRoH64arTS7ie0rgZWnel6XrQQnhBCdjKwEJ4QQogGPTg5CCCHah0cnB+mQFkKI9qG0PqX7x9yKUuoosBcIBVqbKSKB3FrvrX2Oljynldtas7818TYX38leuzrW+ueT79ac77apz5XvtvN8tz211lGNnlVr3WEewOJTeM/62u+tfY6WPG/ttvaKt7n4Tvba1bHKd+se321Tnyvfbef7bht7eHSzUiOanG6jFe/9uJFtzT1v7bbW7G/Je08Wd0tetyQW+W4bP68nf7dNfa58t3Xf2xm+2wY6RLPS6VBKrddae8wc3Z4UryfFCp4VryfFCp4VryfFCu0Xb0erOZyKxWYH0EqeFK8nxQqeFa8nxQqeFa8nxQrtFG+nrzkIIYRoSGoOQgghGpDkIIQQogFJDkIIIRqQ5NAMpVQfpdTLSql3zY6lJZRSFyqlXlJKvaWUOsfseJqjlEpQSr2glHpXKXWD2fGcjFIqUCm1Xik11exYTkYpla6U+tbx/aabHU9zlFJeSqlFSqlnlFJXmR3PySilxju+138qpf5ndjzNUUr1UEp9qJRaopS6u7Xv77DJwfGFHFFKZdbbPkUptV0ptfNkX5jWerfW+pr2jdQZV1vE+6HW+jrgemCWm8eapbW+HrgEGOfOsTrcBbzdPlHWiast4tVAMWDFWH7XnWOdjrHQV1V7xuqIqy3+br91/N3+G3jNnWMFkoF3tda/A4a2OojW3o3nKQ9gAjAMyKy1zQLsAvoAvsDPQKLjS/x3vUd0rfe962Hx/h8wzN1jBaYBnwJz3DlWYDLGYlRzganu/ncAeDne1xV4w81jvRv4vSv+n7Xx/7G3gWB3jhXoAnwNfAVc3eoY2vMfw+wH0KvelzsG+KzW63uAe1pwnnZPDm0RL6CAR4Cz3T3Weuf6xJ1jBRYBTwKfAx/VFL7uGm+t43xdUOCe7nd7OXCJ4/lb7RlrW323QA/gJXePFbgdmOB43uq/A5cs9uNGugH7ar3OAUY1dbBSqgtGwTBUKXWP1vqhdo6vvlbFC9wEnA2EKqX6aa1faM/g6mntd5sOzAD8OI0Fn05Rq2LVWt8LoJSaC+Rqre3tGl1Drf1uZwDnAmHAs+0aWUOt/Zt9H3hGKTUeyGjPwJrQ2ngBrgFeabeImtbaWP8D3K+UmgPsae2Hdbbk0Cpa6zyM9nuPoLV+Gnja7DhaQmu9ClhlchitorV+1ewYWkJr/T5Goev2tNalGIWtx9Ba32d2DC2htc4EZp7q+ztsh3QT9gPda72Od2xzV54Ur8TafjwpXk+KFTwrXpfG2tmSww9Af6VUb6WUL0Yn4wqTY2qOJ8UrsbYfT4rXk2IFz4rXtbG2d6eKWQ9gOXCQE0PkrnFsPx/YgdHrf6/ZcXpivBKrxOtpsXpavO4Qq0y8J4QQooHO1qwkhBCiBSQ5CCGEaECSgxBCiAYkOQghhGhAkoMQQogGJDkIIYRoQJKDEEKIBiQ5CCGEaECSgxDtxLEwy0bHY51SSv6/CY8hd0gL0U6UUtkY8+kfNDsWIVpLrmSEaD8rgU1KqSfNDkSI1pL1HIRoB0qpsRgr88VqravNjkeI1pKagxDt47fADq11tTKEmB2QEK0hfQ5CtAOl1EjgZUADZcCNWusN5kYlRMtJchBCCNGANCsJIYRoQJKDEEKIBiQ5CCGEaECSgxBCiAYkOQghhGhAkoMQQogGJDkIIYRoQJKDEEKIBv4/V1MtDKK/vF8AAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"def CG_cond(A, b, eps, k_max = 1000000):\n",
|
|
" D = np.diag(np.diag(A))\n",
|
|
" C = np.sqrt(linalg.inv(D))\n",
|
|
" A_tilde = C @ A @ C\n",
|
|
" b_tilde = C @ b\n",
|
|
" \n",
|
|
" x_tilde, k = CG(A_tilde, b_tilde, eps, k_max)\n",
|
|
" x = C @ x_tilde\n",
|
|
" return x, k\n",
|
|
"\n",
|
|
"# Sorry for copying.\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"\n",
|
|
"for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\"), (CG_cond, \"Conjugate Gradient\")]:\n",
|
|
" k_list = []\n",
|
|
" for eps in eps_list:\n",
|
|
" x, k = alg(A, b, eps)\n",
|
|
" k_list.append(k)\n",
|
|
" ax.plot(eps_list, k_list, label=alg_name)\n",
|
|
" \n",
|
|
"ax.set_xscale(\"log\")\n",
|
|
"ax.invert_xaxis()\n",
|
|
"ax.set_yscale(\"log\")\n",
|
|
"ax.set_xlabel(\"$\\epsilon$\")\n",
|
|
"ax.set_ylabel(\"$k$\")\n",
|
|
"ax.legend()\n",
|
|
"\n",
|
|
"fig.show()\n",
|
|
"\n",
|
|
"print(\"The number of iterations is brought down by a lot due to the conditioning.\")"
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"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
|
|
}
|