diff --git a/Week 3/7 Linear Equation Systems.ipynb b/Week 3/7 Linear Equation Systems.ipynb index 016726c..eb19eb6 100644 --- a/Week 3/7 Linear Equation Systems.ipynb +++ b/Week 3/7 Linear Equation Systems.ipynb @@ -70,7 +70,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 39, "metadata": { "deletable": false, "nbgrader": { @@ -88,7 +88,8 @@ "outputs": [], "source": [ "import numpy as np\n", - "import numpy.linalg as linalg" + "import numpy.linalg as linalg\n", + "from matplotlib import pyplot as plt" ] }, { @@ -268,7 +269,13 @@ }, "outputs": [], "source": [ - "def GS(A, b, eps):\n", + "def GS(A, b, eps, k_max = 10000):\n", + " \"\"\"\n", + " Return the estimate solution x to the problem Ax = b and the number\n", + " of iterations k it took to decrease maximum norm error below eps\n", + " or to exceed mi\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", @@ -341,7 +348,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 61, "metadata": { "deletable": false, "nbgrader": { @@ -356,6 +363,62 @@ "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": [ + "
" + ] + }, + "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": 62, + "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", @@ -367,7 +430,7 @@ " [ 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", + "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 [ 0.99463393 1.99776509 -0.99803257 1.00108402]\n", @@ -392,71 +455,36 @@ ] } ], - "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", - "\n", - "print(\"Starting with A =\")\n", - "print(A)\n", - "print(\"and b =\", b)\n", - "print(\"We apply the Gauss-Seidel algorithm to solve Ax = b.\")\n", - "print()\n", - "\n", - "eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n", - "for eps in eps_list:\n", - " x, k = GS(A, b, eps)\n", - " print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations:\".format(eps, k))\n", - " print(\"x =\\t\\t\", x)\n", - " print(\"Ax =\\t\\t\", np.dot(A, x))\n", - " print(\"diff(Ax, b) =\\t\", diff(A @ x, b))\n", - " print()" - ] - }, - { - "cell_type": "code", - "execution_count": 63, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "1.9462975206611572" - ] - }, - "execution_count": 63, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "diff(x,np.array([1,1,1,1]))" - ] - }, - { - "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", + " 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", + " 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 Gauss-Seidel algorithm to solve Ax = b.\")\n", + " print()\n", + "\n", + " eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n", + " for eps in eps_list:\n", + " x, k = GS(A, b, eps)\n", + " print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations:\".format(eps, k))\n", + " print(\"x =\\t\\t\", x)\n", + " print(\"Ax =\\t\\t\", np.dot(A, x))\n", + " print(\"diff(Ax, b) =\\t\", diff(A @ x, b))\n", + " print()\n", + " \n", + " assert diff(x, x_exact) < eps\n", " \n", "test_GS()" ]