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

07: Draft task 3

Browse Source
This commit is contained in:
Kees van Kempen
2022-02-18 12:33:50 +01:00
parent 73466c7506
commit cc133e99ef
1 changed files with 100 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -87,7 +87,8 @@
},
"outputs": [],
"source": [
"import numpy as np"
"import numpy as np\n",
"import numpy.linalg as linalg"
]
},
{
@ -133,7 +134,7 @@
"outputs": [],
"source": [
"def diff(a, b):\n",
" return np.max(a - b)"
" return np.max(np.abs(a - b)"
]
},
{
@ -248,7 +249,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 56,
"metadata": {
"deletable": false,
"nbgrader": {
@ -266,8 +267,31 @@
"outputs": [],
"source": [
"def GS(A, b, eps):\n",
" # YOUR CODE HERE\n",
" raise NotImplementedError()"
" # 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:\n",
" k += 1\n",
" x_prev = x_cur\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",
" return x_cur, k"
]
},
{
@ -316,7 +340,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 57,
"metadata": {
"deletable": false,
"nbgrader": {
@ -331,18 +355,78 @@
"task": false
}
},
"outputs": [],
"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",
"[ 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"
]
}
],
"source": [
"# Here, Koen, I made the matrix.\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",
"# Do your plotting here ...\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",
"# YOUR CODE HERE\n",
"raise NotImplementedError()"
"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))"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.7565702479338843"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(x,np.array([1,1,1,1]))"
]
},
{