From a4a38f2c3770384a90e279e26bdca7c31827eea6 Mon Sep 17 00:00:00 2001 From: Kees van Kempen Date: Fri, 18 Feb 2022 15:48:59 +0100 Subject: [PATCH] 07: Explicitly copy array, fix typos: GS works! Under the hood, Numpy arrays are C arrays, and therefore explicitly copying is needed to create copies and not have them point to the same memory. Ugh ugh ugh this cost me 1.5 hour. --- Week 3/7 Linear Equation Systems.ipynb | 90 ++++++++++++++------------ 1 file changed, 48 insertions(+), 42 deletions(-) diff --git a/Week 3/7 Linear Equation Systems.ipynb b/Week 3/7 Linear Equation Systems.ipynb index b878031..016726c 100644 --- a/Week 3/7 Linear Equation Systems.ipynb +++ b/Week 3/7 Linear Equation Systems.ipynb @@ -116,7 +116,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 2, "metadata": { "deletable": false, "nbgrader": { @@ -134,7 +134,7 @@ "outputs": [], "source": [ "def diff(a, b):\n", - " return np.max(np.abs(a - b)" + " return np.max(np.abs(a - b))" ] }, { @@ -189,11 +189,13 @@ } }, "source": [ + "---\n", + "\n", "We start from our linear equations:\n", "\n", "$$Ax = b$$\n", "\n", - "We seperate A into different components (diagonal, strictly lower triangular and strictly upper triangular):\n", + "We separate A into different components (diagonal, strictly lower triangular and strictly upper triangular):\n", "\n", "$$A = D - L - U$$\n", "\n", @@ -211,7 +213,7 @@ "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 rewrite our previous equation to:\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].$$" @@ -249,7 +251,7 @@ }, { "cell_type": "code", - "execution_count": 56, + "execution_count": 21, "metadata": { "deletable": false, "nbgrader": { @@ -285,12 +287,11 @@ " k = 1\n", " while diff(x_cur, x_prev) > eps:\n", " k += 1\n", - " x_prev = x_cur\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] * ( - sum([A[i, j] * x_cur[j] for j in range(i)]) - sum([A[i, j] * x_cur[j] for j in range(i + 1, n)]) + b[i])\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", - " \n", - " print(x_cur, x_prev)\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" ] }, @@ -340,7 +341,7 @@ }, { "cell_type": "code", - "execution_count": 57, + "execution_count": 30, "metadata": { "deletable": false, "nbgrader": { @@ -368,59 +369,64 @@ "and b = [ 6 25 -11 15]\n", "We apply the Gauss-Seidel algorithm to solve Ax = b\n", "\n", - "[ 1.04727273 2.27272727 -1.1 1.875 ] [ 1.04727273 2.27272727 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -1.1 1.875 ] [ 1.04727273 1.75657025 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.875 ] [ 1.04727273 1.75657025 -0.94629752 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.09799897] [ 1.04727273 1.75657025 -0.94629752 1.09799897]\n", - "For eps = 1e-01\tafter k = 2\t iterations, x = [ 1.04727273 1.75657025 -0.94629752 1.09799897] \twhich gives Ax = [ 6.82356198 22.51529442 -10.22299897 15. ]\n", - "[ 1.04727273 2.27272727 -1.1 1.875 ] [ 1.04727273 2.27272727 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -1.1 1.875 ] [ 1.04727273 1.75657025 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.875 ] [ 1.04727273 1.75657025 -0.94629752 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.09799897] [ 1.04727273 1.75657025 -0.94629752 1.09799897]\n", - "For eps = 1e-02\tafter k = 2\t iterations, x = [ 1.04727273 1.75657025 -0.94629752 1.09799897] \twhich gives Ax = [ 6.82356198 22.51529442 -10.22299897 15. ]\n", - "[ 1.04727273 2.27272727 -1.1 1.875 ] [ 1.04727273 2.27272727 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -1.1 1.875 ] [ 1.04727273 1.75657025 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.875 ] [ 1.04727273 1.75657025 -0.94629752 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.09799897] [ 1.04727273 1.75657025 -0.94629752 1.09799897]\n", - "For eps = 1e-03\tafter k = 2\t iterations, x = [ 1.04727273 1.75657025 -0.94629752 1.09799897] \twhich gives Ax = [ 6.82356198 22.51529442 -10.22299897 15. ]\n", - "[ 1.04727273 2.27272727 -1.1 1.875 ] [ 1.04727273 2.27272727 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -1.1 1.875 ] [ 1.04727273 1.75657025 -1.1 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.875 ] [ 1.04727273 1.75657025 -0.94629752 1.875 ]\n", - "[ 1.04727273 1.75657025 -0.94629752 1.09799897] [ 1.04727273 1.75657025 -0.94629752 1.09799897]\n", - "For eps = 1e-04\tafter k = 2\t iterations, x = [ 1.04727273 1.75657025 -0.94629752 1.09799897] \twhich gives Ax = [ 6.82356198 22.51529442 -10.22299897 15. ]\n" + "For eps = 1e-01\tafter k = 4\t iterations:\n", + "x =\t\t [ 0.99463393 1.99776509 -0.99803257 1.00108402]\n", + "Ax =\t\t [ 5.95250909 24.98206671 -10.98990699 15. ]\n", + "diff(Ax, b) =\t 0.04749090616931895\n", + "\n", + "For eps = 1e-02\tafter k = 5\t iterations:\n", + "x =\t\t [ 0.99938302 1.99982713 -0.99978549 1.00009164]\n", + "Ax =\t\t [ 5.99443213 24.99877578 -10.99900762 15. ]\n", + "diff(Ax, b) =\t 0.005567865937722516\n", + "\n", + "For eps = 1e-03\tafter k = 6\t iterations:\n", + "x =\t\t [ 0.99993981 1.99998904 -0.99997989 1.00000662]\n", + "Ax =\t\t [ 5.99944928 24.99993935 -10.99991498 15. ]\n", + "diff(Ax, b) =\t 0.000550717702960668\n", + "\n", + "For eps = 1e-04\tafter k = 7\t iterations:\n", + "x =\t\t [ 0.99999488 1.99999956 -0.99999836 1.00000037]\n", + "Ax =\t\t [ 5.99995255 24.99999971 -10.99999375 15. ]\n", + "diff(Ax, b) =\t 4.744782363452771e-05\n", + "\n" ] } ], "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(\"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", - " 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, k = GS(A, b, eps)\n", - " print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations, x =\".format(eps, k), x, \"\\twhich gives Ax =\", np.dot(A, x))" + " 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": 53, + "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "0.7565702479338843" + "1.9462975206611572" ] }, - "execution_count": 53, + "execution_count": 63, "metadata": {}, "output_type": "execute_result" } @@ -764,7 +770,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" },