kees-study/cds-numerical-methods
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

07: Explicitly copy array, fix typos: GS works!

Browse Source 
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.
This commit is contained in:
Kees van Kempen
2022-02-18 15:48:59 +01:00
parent cc133e99ef
commit a4a38f2c37
1 changed files with 48 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

90
Week 3/7 Linear Equation Systems.ipynb
View File

@ -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"
},