diff --git a/Week 4/8 Eigenvalues and Eigenvectors.ipynb b/Week 4/8 Eigenvalues and Eigenvectors.ipynb index d810302..ef17fe3 100644 --- a/Week 4/8 Eigenvalues and Eigenvectors.ipynb +++ b/Week 4/8 Eigenvalues and Eigenvectors.ipynb @@ -39,7 +39,7 @@ "cell_type": "raw", "metadata": {}, "source": [ - "team_members = \"\"" + "team_members = \"Koen Vendrig, Kees van Kempen\"" ] }, { @@ -66,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "deletable": false, "nbgrader": { @@ -83,10 +83,8 @@ }, "outputs": [], "source": [ - "# Import packages here ...\n", - "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "import numpy as np\n", + "import numpy.linalg as linalg" ] }, { @@ -139,7 +137,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": { "deletable": false, "nbgrader": { @@ -157,13 +155,53 @@ "outputs": [], "source": [ "def inversePower(A, sigma, eps):\n", - " # YOUR CODE HERE\n", - " raise NotImplementedError()" + " \"\"\"\n", + " Estimates the eigenvectors of matrix A by the inverse power method.\n", + "\n", + " Args:\n", + " A: an n x n matrix\n", + " sigma: an initial guess for an eigenvector\n", + " eps: the desired accuracy\n", + "\n", + " Returns:\n", + " A tuple (b, k) is returned, containing:\n", + " b: the eigenvector b corresponding to the eigenvalue\n", + " closests to sigma after k iterations\n", + " k: the number of iterations done\n", + " \n", + " See also:\n", + " https://www.sciencedirect.com/topics/mathematics/inverse-power-method\n", + " \"\"\"\n", + " \n", + " # Does https://johnfoster.pge.utexas.edu/numerical-methods-book/LinearAlgebra_EigenProblem1.html help?\n", + " \n", + " # A should be n x n.\n", + " n = len(A)\n", + " assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n", + " \n", + " # Setup some initial values.\n", + " #B = linalg.inv(A - sigma*np.ones(n))\n", + " B = linalg.inv(A - sigma*np.eye(n))\n", + " #B = 1/(A - sigma*np.ones(n))\n", + " #B = 1/(A - sigma*np.eye(n))\n", + " b = np.ones(n)\n", + " k = 0\n", + " e = 0\n", + " \n", + " while e > eps or k == 0:\n", + " b_prev = b.copy()\n", + " k += 1\n", + " \n", + " b = B @ b\n", + " b /= np.sqrt(b @ b) # although b = linalg.norm(b) could be used\n", + " e = np.sqrt(np.sum( (np.abs(b_prev) - np.abs(b))**2) )\n", + " \n", + " return b, k" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": { "deletable": false, "nbgrader": { @@ -178,12 +216,48 @@ "task": false } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For sigma = 0.03\n", + "b = [ 0.45440139 -0.76618454 0.45440139]\n", + "lam = [0.62771919 0.62771831 0.62771919]\n", + "\n", + "\n", + "For sigma = 6.0\n", + "b = [0.54177432 0.64262054 0.54177432]\n", + "lam = [6.37228128 6.37228139 6.37228128]\n", + "\n", + "\n", + "For sigma = 1.99999989999999\n", + "b = [-7.07106781e-01 -4.21799612e-14 7.07106781e-01]\n", + "lam = [ 2. -0. 2.]\n", + "\n", + "\n", + "[6.37228132 2. 0.62771868]\n" + ] + } + ], "source": [ "# Use this cell for your own testing ...\n", "\n", - "# YOUR CODE HERE\n", - "raise NotImplementedError()" + "A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])\n", + "\n", + "sigma_list = [.03, 6., 1.99999989999999]\n", + "#sigma = 0.03\n", + "\n", + "for sigma in sigma_list:\n", + " print(\"For sigma =\", sigma)\n", + " b, k = inversePower(A, sigma, 1e-6)\n", + " print(\"b =\", b)\n", + " lam = np.dot(A, b) / b\n", + " print(\"lam =\", lam)\n", + " print()\n", + " print()\n", + "\n", + "print(linalg.eig(A)[0])" ] }, { @@ -206,6 +280,7 @@ "outputs": [], "source": [ "def test_inversePower():\n", + " A = np.array([[3, 2, 1], [2, 3, 2], [1, 2, 3]])\n", " # YOUR CODE HERE\n", " raise NotImplementedError()\n", " \n", @@ -506,7 +581,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" },