diff --git a/Week 3/7 Linear Equation Systems.ipynb b/Week 3/7 Linear Equation Systems.ipynb index a318112..b878031 100644 --- a/Week 3/7 Linear Equation Systems.ipynb +++ b/Week 3/7 Linear Equation Systems.ipynb @@ -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]))" ] }, {