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

08: Task 1 somewhat working. What to do with zeroes?

Browse Source
This commit is contained in:
Kees van Kempen
2022-02-25 11:44:51 +01:00
parent b6f0b0d7d9
commit ab959984a6
1 changed files with 89 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

103
Week 4/8 Eigenvalues and Eigenvectors.ipynb
View File

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