{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 2\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 3 0 -1 0]\n", " [ 0 1 0 1]\n", " [-1 0 8 1]\n", " [ 0 1 1 1]]\n" ] } ], "source": [ "A = np.array([[3, 0, -1, 0],[0, 1, 0, 1], [-1, 0, 8, 1],[0, 1, 1, 1]]);\n", "print(A)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# This powermethod function outputs a pair l, w where l is the the eigenvalue and w is the eigenvector.\n", "\n", "def powermethod(A, v, n): \n", " l = np.zeros(n)\n", " w=(1/np.linalg.norm(v))*v\n", " for k in range(n):\n", " v = A @ w\n", " newl = np.dot(v,w)\n", " l[k] = newl\n", " w=(1/np.linalg.norm(v))*v\n", " return l,w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the powermethod to find $e_1$ and $\\lambda_1$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lambda1 = 8.32680562975727 \n", " v_1 = [-0.18277562 0.01848086 0.97361021 0.13540568]\n" ] } ], "source": [ "n,n = np.shape(A) # should be square\n", "v0 = np.random.rand(n)\n", "elist1, v1 = powermethod(A,v0,20)\n", "l1 = elist1[-1]\n", "\n", "print(\"Lambda1 =\",l1, \"\\n v_1 = \", v1)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ True, True, True, True])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isclose(A @ v1, l1*v1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The First Eigenvector/Eigenvalue pair \n", "So we have $e_1 \\approx \\begin{pmatrix}-0.182 \\\\ 0.018 \\\\ 0.973 \\\\ 0.135\\end{pmatrix}$ and $\\lambda_1 \\approx 8.33$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# this does the powermethod on A-sI\n", "def powermethodshift(A, v, s, n):\n", " m,m = A.shape\n", " B = A - s*np.eye(m,m)\n", " l = np.zeros(n)\n", " w=(1/np.linalg.norm(v))*v\n", " for k in range(n):\n", " v = B @ w\n", " newl = np.dot(v,w)\n", " l[k] = newl\n", " w=(1/np.linalg.norm(v))*v\n", " return l+ s,w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we set $B=A-\\lambda_1 I$ then $B$ has the same eigenvectors as $A$ but the eigenvalues are shifted back by $\\lambda_1$. This means that the principal eigenvector of $B$ is gauranteed to be something other than $e_1$ ($e_1$ is an eigenvector of $B$ with eigenvalue $0$). So, we'll find a second eigenvector of $A$ this way." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lambda2 = -0.0666785904540319 \n", " e_2 = [ 0.03058325 0.68060099 0.093789 -0.72598249]\n" ] } ], "source": [ "elist2, v2 = powermethodshift(A,v0,l1,80)\n", "l2 = elist2[-1] \n", "\n", "print(\"Lambda2 =\",l2, \"\\n e_2 = \", v2)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1.74469706e-09 1.08596754e-08 -1.90512418e-09 9.86121242e-09]\n" ] } ], "source": [ "print(A @ v2 - l2*v2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Second Eigenvector/Eigenvalue pair \n", "So we have $e_2 \\approx \\begin{pmatrix}0.035 \\\\ 0.681 \\\\ 0.94 \\\\ -0.726\\end{pmatrix}$ and $\\lambda_1 \\approx -0.067$" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# this does the powermethod on (A-sI)^{-1}\n", "def invpowermethodshift(A, v, s, n):\n", " m,m = A.shape\n", " B = np.linalg.inv(A - s*np.eye(m,m))\n", " l = np.zeros(n)\n", " w=(1/np.linalg.norm(v))*v\n", " for k in range(n):\n", " v = B @ w\n", " newl = 1/np.dot(v,w)\n", " l[k] = newl\n", " w=(1/np.linalg.norm(v))*v\n", " return l+s,w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The operator $(A-sI)^{-1}$ has the same eigenvectors as $A$, but the corresponding eigenvalues are $\\frac{1}{\\lambda-s}$, which has the greatest magnitude when $\\lambda \\approx s$. So, the power method on this operator will find the eigenvector of $A$ whose eigenvalue is closest to $s$. We can look for an eigenvector whose eigenvalue is imbetween $\\lambda_1$ and $\\lambda_2$ by using this shifted inverse power method function with shift equal to the midpoint of $\\lambda_1$ and $\\lambda_2$. " ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "elist3, v3 = invpowermethodshift(A,v0,(l1+l2)/2,80)\n", "l3 = elist3[-1] " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.8318153216290556 \n", " [0.97567239 0.06966243 0.16409315 0.1276087 ]\n" ] } ], "source": [ "print(l3, \"\\n\", v3)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-4.44089210e-16 2.77555756e-17 0.00000000e+00 0.00000000e+00]\n" ] } ], "source": [ "print(A @ v3 - l3*v3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Third Eigenvector/Eigenvalue pair \n", "So we have $e_3 \\approx \\begin{pmatrix}0.976 \\\\ 0.070 \\\\ 0.164 \\\\ 0.128\\end{pmatrix}$ and $\\lambda_1 \\approx 2.832$\n", "\n", "Let's take stock" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8.32680562975727 -0.0666785904540319 2.8318153216290556\n" ] } ], "source": [ "print(l1, l2, l3)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGdCAYAAADuR1K7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfA0lEQVR4nO3de3TT9f3H8VdSaMuUhFGgoVIEPZ4VBKm2UIrb9Nj+rJfN9YgTOSiIPXLmryBYxmwRYZuXejkosiIdnm3OAz0iOhky1h1WdkBH5dKCk6uesw0QfklhSAJllNrk9wczM5qWgqRp3zwf5+Rw+Obzbd5fE5PnSb8JjlAoFBIAAIARzngPAAAAcCERNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADClW7wHiIdgMKhDhw6pZ8+ecjgc8R4HAAC0QygU0vHjx5WWlians/X3Zy7KuDl06JDS09PjPQYAADgPBw4c0IABA1q9/qKMm549e0o68x/H5XLFeRoAANAegUBA6enp4dfx1lyUcfP5r6JcLhdxAwBAF3O2U0o4oRgAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYApxAwAATCFuAACAKcQNAAAwhbgBAACmEDcAAMAU4gYAAJhC3AAAAFOIGwAAYEqHxM2iRYs0aNAgJScnKycnR5s3b25z/YoVK5SRkaHk5GQNHz5ca9asaXXtj370IzkcDi1YsOACTw0AALqimMfN8uXLVVJSonnz5qm+vl4jRoxQQUGBGhoaoq7fuHGjxo8fr6KiIm3btk2FhYUqLCzUjh07vrL27bff1vvvv6+0tLRYHwYAAOgiYh43L7zwgh588EFNnjxZQ4cOVWVlpb7xjW/o17/+ddT1L730km655RbNmjVLQ4YM0RNPPKHrrrtOFRUVEesOHjyoadOmadmyZerevXusDwMAAHQRMY2b06dPq66uTvn5+f+9QadT+fn5qq2tjbpPbW1txHpJKigoiFgfDAZ13333adasWbr66qvPOkdTU5MCgUDEBQAA2BTTuDly5IhaWlqUmpoasT01NVVerzfqPl6v96zrn332WXXr1k0PP/xwu+YoLy+X2+0OX9LT08/xSAAAQFfR5T4tVVdXp5deekmvvvqqHA5Hu/YpKyuT3+8PXw4cOBDjKQEAQLzENG769OmjhIQE+Xy+iO0+n08ejyfqPh6Pp8317777rhoaGjRw4EB169ZN3bp10759+zRz5kwNGjQo6s9MSkqSy+WKuAAAAJtiGjeJiYnKyspSTU1NeFswGFRNTY1yc3Oj7pObmxuxXpLWrl0bXn/ffffpb3/7m7Zv3x6+pKWladasWfrTn/4Uu4MBAABdQrdY30BJSYkmTZqk7OxsjRo1SgsWLFBjY6MmT54sSZo4caIuu+wylZeXS5KmT5+uG264QfPnz9ftt9+u119/XVu3btWSJUskSSkpKUpJSYm4je7du8vj8ehb3/pWrA8HAAB0cjGPm3Hjxunw4cOaO3euvF6vMjMzVV1dHT5peP/+/XI6//sG0pgxY1RVVaU5c+Zo9uzZuuqqq7Ry5UoNGzYs1qMCAAADHKFQKBTvITpaIBCQ2+2W3+/n/BsAALqI9r5+d7lPSwEAALSFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApHRI3ixYt0qBBg5ScnKycnBxt3ry5zfUrVqxQRkaGkpOTNXz4cK1ZsyZ8XXNzsx599FENHz5cl1xyidLS0jRx4kQdOnQo1ocBAAC6gJjHzfLly1VSUqJ58+apvr5eI0aMUEFBgRoaGqKu37hxo8aPH6+ioiJt27ZNhYWFKiws1I4dOyRJJ0+eVH19vR5//HHV19frd7/7nfbu3as77rgj1ocCAAC6AEcoFArF8gZycnI0cuRIVVRUSJKCwaDS09M1bdo0lZaWfmX9uHHj1NjYqNWrV4e3jR49WpmZmaqsrIx6G1u2bNGoUaO0b98+DRw48KwzBQIBud1u+f1+uVyu8zwyAADQkdr7+h3Td25Onz6turo65efn//cGnU7l5+ertrY26j61tbUR6yWpoKCg1fWS5Pf75XA41KtXr6jXNzU1KRAIRFwAAIBNMY2bI0eOqKWlRampqRHbU1NT5fV6o+7j9XrPaf2pU6f06KOPavz48a1WXHl5udxud/iSnp5+HkcDAAC6gi79aanm5mbdfffdCoVCWrx4cavrysrK5Pf7w5cDBw504JQAAKAjdYvlD+/Tp48SEhLk8/kitvt8Pnk8nqj7eDyedq3/PGz27dundevWtfm7t6SkJCUlJZ3nUQAAgK4kpu/cJCYmKisrSzU1NeFtwWBQNTU1ys3NjbpPbm5uxHpJWrt2bcT6z8Pm448/1p///GelpKTE5gAAAECXE9N3biSppKREkyZNUnZ2tkaNGqUFCxaosbFRkydPliRNnDhRl112mcrLyyVJ06dP1w033KD58+fr9ttv1+uvv66tW7dqyZIlks6EzV133aX6+nqtXr1aLS0t4fNxevfurcTExFgfEgAA6MRiHjfjxo3T4cOHNXfuXHm9XmVmZqq6ujp80vD+/fvldP73DaQxY8aoqqpKc+bM0ezZs3XVVVdp5cqVGjZsmCTp4MGDWrVqlSQpMzMz4rb+8pe/6MYbb4z1IQEAgE4s5t9z0xnxPTcAAHQ9neJ7bgAAADoacQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACY0i3eA5gSbJH2bZRO+KRLU6XLx0jOhHhPBYt4rAHojDrJc1OHvHOzaNEiDRo0SMnJycrJydHmzZvbXL9ixQplZGQoOTlZw4cP15o1ayKuD4VCmjt3rvr3768ePXooPz9fH3/8cSwP4ex2rZIWDJN++z3praIzfy4YdmY7cCHxWAPQGXWi56aYx83y5ctVUlKiefPmqb6+XiNGjFBBQYEaGhqirt+4caPGjx+voqIibdu2TYWFhSosLNSOHTvCa5577jktXLhQlZWV2rRpky655BIVFBTo1KlTsT6c6Hatkt6YKAUORW4P/N+Z7bzo4ELhsQagM+pkz02OUCgUiuUN5OTkaOTIkaqoqJAkBYNBpaena9q0aSotLf3K+nHjxqmxsVGrV68Obxs9erQyMzNVWVmpUCiktLQ0zZw5Uz/+8Y8lSX6/X6mpqXr11Vd1zz33nHWmQCAgt9stv98vl8v19Q4w2HKmTL98h/5HSA4Fe/bXpw/W8WsDfD3BFn1zSZacJw7JEeXqkBwKudJ06n+38VgD0HGCLUpelCnH8ejPTZJDcqVJMz782s9N7X39juk5N6dPn1ZdXZ3KysrC25xOp/Lz81VbWxt1n9raWpWUlERsKygo0MqVKyVJ//jHP+T1epWfnx++3u12KycnR7W1tVHjpqmpSU1NTeG/BwKBr3NYkfZtbDVsJMmhkBKOH9LUZxbp/eDQC3e7uOiMdu7S64ltP9YcgYN64OcLeawB6DBne26SQlLg4JnXy8Hf6ZCZYvprqSNHjqilpUWpqakR21NTU+X1eqPu4/V621z/+Z/n8jPLy8vldrvDl/T09PM6nqhO+Nq1rJ+OXbjbxEWpvY8hHmsAOlK7n3Pa+Xp5IVwUn5YqKyuLeDcoEAhcuMC5NPXsayQ9cW+e5g68/sLcJi5K3ff3kN6oOOu65yb/j565/NsdMBEASM59l0jLzv7c1N7XywshpnHTp08fJSQkyOeLrDWfzyePxxN1H4/H0+b6z//0+Xzq379/xJrMzMyoPzMpKUlJSUnnexhtu3zMmd8lBv5PUrTTl878rtGdcQPnQeDrybihXY+15Cu/w2MNQMe58jvtem7S5WM6bKSY/loqMTFRWVlZqqmpCW8LBoOqqalRbm5u1H1yc3Mj1kvS2rVrw+sHDx4sj8cTsSYQCGjTpk2t/syYciZItzz7n798+VSq//z9lmd4scHXx2MNQGfUCZ+bYv5R8JKSEr3yyiv67W9/q927d+uhhx5SY2OjJk+eLEmaOHFixAnH06dPV3V1tebPn689e/bopz/9qbZu3aqpU6dKkhwOh2bMmKEnn3xSq1at0ocffqiJEycqLS1NhYWFsT6c6IbeId39muTqH7ndlXZm+9A74jMX7OGxBqAz6mTPTTE/52bcuHE6fPiw5s6dK6/Xq8zMTFVXV4dPCN6/f7+czv821pgxY1RVVaU5c+Zo9uzZuuqqq7Ry5UoNGzYsvOYnP/mJGhsbNWXKFB07dkzf/va3VV1dreTk5FgfTuuG3iFl3N4pvpkRxvFYA9AZdaLnpph/z01ndEG/5wYAAHSI9r5+8w9nAgAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmBKzuDl69KgmTJggl8ulXr16qaioSCdOnGhzn1OnTqm4uFgpKSm69NJLNXbsWPl8vvD1H3zwgcaPH6/09HT16NFDQ4YM0UsvvRSrQwAAAF1QzOJmwoQJ2rlzp9auXavVq1drw4YNmjJlSpv7PPLII3rnnXe0YsUKrV+/XocOHdKdd94Zvr6urk79+vXT0qVLtXPnTj322GMqKytTRUVFrA4DAAB0MY5QKBS60D909+7dGjp0qLZs2aLs7GxJUnV1tW677TZ98sknSktL+8o+fr9fffv2VVVVle666y5J0p49ezRkyBDV1tZq9OjRUW+ruLhYu3fv1rp169o9XyAQkNvtlt/vl8vlOo8jBAAAHa29r98xeeemtrZWvXr1CoeNJOXn58vpdGrTpk1R96mrq1Nzc7Py8/PD2zIyMjRw4EDV1ta2elt+v1+9e/e+cMMDAIAurVssfqjX61W/fv0ib6hbN/Xu3Vter7fVfRITE9WrV6+I7ampqa3us3HjRi1fvlx/+MMf2pynqalJTU1N4b8HAoF2HAUAAOiKzumdm9LSUjkcjjYve/bsidWsEXbs2KEf/OAHmjdvnm6++eY215aXl8vtdocv6enpHTIjAADoeOf0zs3MmTN1//33t7nmiiuukMfjUUNDQ8T2zz77TEePHpXH44m6n8fj0enTp3Xs2LGId298Pt9X9tm1a5fy8vI0ZcoUzZkz56xzl5WVqaSkJPz3QCBA4AAAYNQ5xU3fvn3Vt2/fs67Lzc3VsWPHVFdXp6ysLEnSunXrFAwGlZOTE3WfrKwsde/eXTU1NRo7dqwkae/evdq/f79yc3PD63bu3KmbbrpJkyZN0lNPPdWuuZOSkpSUlNSutQAAoGuLyaelJOnWW2+Vz+dTZWWlmpubNXnyZGVnZ6uqqkqSdPDgQeXl5em1117TqFGjJEkPPfSQ1qxZo1dffVUul0vTpk2TdObcGunMr6JuuukmFRQU6Pnnnw/fVkJCQrui63N8WgoAgK6nva/fMTmhWJKWLVumqVOnKi8vT06nU2PHjtXChQvD1zc3N2vv3r06efJkeNuLL74YXtvU1KSCggK9/PLL4evffPNNHT58WEuXLtXSpUvD2y+//HL985//jNWhAACALiRm79x0ZrxzAwBA1xPX77kBAACIF+IGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAlJjFzdGjRzVhwgS5XC716tVLRUVFOnHiRJv7nDp1SsXFxUpJSdGll16qsWPHyufzRV37r3/9SwMGDJDD4dCxY8dicAQAAKArilncTJgwQTt37tTatWu1evVqbdiwQVOmTGlzn0ceeUTvvPOOVqxYofXr1+vQoUO68847o64tKirSNddcE4vRAQBAF+YIhUKhC/1Dd+/eraFDh2rLli3Kzs6WJFVXV+u2227TJ598orS0tK/s4/f71bdvX1VVVemuu+6SJO3Zs0dDhgxRbW2tRo8eHV67ePFiLV++XHPnzlVeXp4+/fRT9erVq93zBQIBud1u+f1+uVyur3ewAACgQ7T39Tsm79zU1taqV69e4bCRpPz8fDmdTm3atCnqPnV1dWpublZ+fn54W0ZGhgYOHKja2trwtl27dunnP/+5XnvtNTmd7Ru/qalJgUAg4gIAAGyKSdx4vV7169cvYlu3bt3Uu3dveb3eVvdJTEz8yjswqamp4X2ampo0fvx4Pf/88xo4cGC75ykvL5fb7Q5f0tPTz+2AAABAl3FOcVNaWiqHw9HmZc+ePbGaVWVlZRoyZIjuvffec97P7/eHLwcOHIjRhAAAIN66ncvimTNn6v77729zzRVXXCGPx6OGhoaI7Z999pmOHj0qj8cTdT+Px6PTp0/r2LFjEe/e+Hy+8D7r1q3Thx9+qDfffFOS9PnpQn369NFjjz2mn/3sZ1F/dlJSkpKSktpziAAAoIs7p7jp27ev+vbte9Z1ubm5OnbsmOrq6pSVlSXpTJgEg0Hl5ORE3ScrK0vdu3dXTU2Nxo4dK0nau3ev9u/fr9zcXEnSW2+9pX//+9/hfbZs2aIHHnhA7777rq688spzORQAAGDUOcVNew0ZMkS33HKLHnzwQVVWVqq5uVlTp07VPffcE/6k1MGDB5WXl6fXXntNo0aNktvtVlFRkUpKStS7d2+5XC5NmzZNubm54U9KfTlgjhw5Er69c/m0FAAAsCsmcSNJy5Yt09SpU5WXlyen06mxY8dq4cKF4eubm5u1d+9enTx5MrztxRdfDK9tampSQUGBXn755ViNCAAADIrJ99x0dnzPDQAAXU9cv+cGAAAgXogbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgCnEDAABMIW4AAIApxA0AADCFuAEAAKYQNwAAwBTiBgAAmELcAAAAU4gbAABgSrd4DxAPoVBIkhQIBOI8CQAAaK/PX7c/fx1vzUUZN8ePH5ckpaenx3kSAABwro4fPy63293q9Y7Q2fLHoGAwqEOHDqlnz55yOBwxu51AIKD09HQdOHBALpcrZreDC4P7q+vgvuo6uK+6ls5+f4VCIR0/flxpaWlyOls/s+aifOfG6XRqwIABHXZ7LperUz5IEB33V9fBfdV1cF91LZ35/mrrHZvPcUIxAAAwhbgBAACmEDcxlJSUpHnz5ikpKSneo6AduL+6Du6rroP7qmuxcn9dlCcUAwAAu3jnBgAAmELcAAAAU4gbAABgCnEDAABMIW5iaNGiRRo0aJCSk5OVk5OjzZs3x3skfEl5eblGjhypnj17ql+/fiosLNTevXvjPRba4ZlnnpHD4dCMGTPiPQpacfDgQd17771KSUlRjx49NHz4cG3dujXeY+FLWlpa9Pjjj2vw4MHq0aOHrrzySj3xxBNn/febOjPiJkaWL1+ukpISzZs3T/X19RoxYoQKCgrU0NAQ79HwBevXr1dxcbHef/99rV27Vs3Nzbr55pvV2NgY79HQhi1btuiXv/ylrrnmmniPglZ8+umnuv7669W9e3f98Y9/1K5duzR//nx985vfjPdo+JJnn31WixcvVkVFhXbv3q1nn31Wzz33nH7xi1/Ee7TzxkfBYyQnJ0cjR45URUWFpDP/nlV6erqmTZum0tLSOE+H1hw+fFj9+vXT+vXr9d3vfjfe4yCKEydO6LrrrtPLL7+sJ598UpmZmVqwYEG8x8KXlJaW6q9//avefffdeI+Cs/je976n1NRU/epXvwpvGzt2rHr06KGlS5fGcbLzxzs3MXD69GnV1dUpPz8/vM3pdCo/P1+1tbVxnAxn4/f7JUm9e/eO8yRoTXFxsW6//faI/7/Q+axatUrZ2dn64Q9/qH79+unaa6/VK6+8Eu+xEMWYMWNUU1Ojjz76SJL0wQcf6L333tOtt94a58nO30X5D2fG2pEjR9TS0qLU1NSI7ampqdqzZ0+cpsLZBINBzZgxQ9dff72GDRsW73EQxeuvv676+npt2bIl3qPgLP7+979r8eLFKikp0ezZs7VlyxY9/PDDSkxM1KRJk+I9Hr6gtLRUgUBAGRkZSkhIUEtLi5566ilNmDAh3qOdN+IG+I/i4mLt2LFD7733XrxHQRQHDhzQ9OnTtXbtWiUnJ8d7HJxFMBhUdna2nn76aUnStddeqx07dqiyspK46WTeeOMNLVu2TFVVVbr66qu1fft2zZgxQ2lpaV32viJuYqBPnz5KSEiQz+eL2O7z+eTxeOI0FdoydepUrV69Whs2bNCAAQPiPQ6iqKurU0NDg6677rrwtpaWFm3YsEEVFRVqampSQkJCHCfEF/Xv319Dhw6N2DZkyBC99dZbcZoIrZk1a5ZKS0t1zz33SJKGDx+uffv2qby8vMvGDefcxEBiYqKysrJUU1MT3hYMBlVTU6Pc3Nw4ToYvC4VCmjp1qt5++22tW7dOgwcPjvdIaEVeXp4+/PBDbd++PXzJzs7WhAkTtH37dsKmk7n++uu/8rUKH330kS6//PI4TYTWnDx5Uk5nZA4kJCQoGAzGaaKvj3duYqSkpESTJk1Sdna2Ro0apQULFqixsVGTJ0+O92j4guLiYlVVVen3v/+9evbsKa/XK0lyu93q0aNHnKfDF/Xs2fMr50JdcsklSklJ4RypTuiRRx7RmDFj9PTTT+vuu+/W5s2btWTJEi1ZsiTeo+FLvv/97+upp57SwIEDdfXVV2vbtm164YUX9MADD8R7tPPGR8FjqKKiQs8//7y8Xq8yMzO1cOFC5eTkxHssfIHD4Yi6/Te/+Y3uv//+jh0G5+zGG2/ko+Cd2OrVq1VWVqaPP/5YgwcPVklJiR588MF4j4UvOX78uB5//HG9/fbbamhoUFpamsaPH6+5c+cqMTEx3uOdF+IGAACYwjk3AADAFOIGAACYQtwAAABTiBsAAGAKcQMAAEwhbgAAgCnEDQAAMIW4AQAAphA3AADAFOIGAACYQtwAAABTiBsAAGDK/wPwsXGy+z0jzQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "xpoints = np.array([l1, l2, l3])\n", "ypoints = np.array([0,0,0])\n", "\n", "\n", "plt.plot(xpoints, ypoints)\n", "plt.plot(xpoints, ypoints, 'o')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looking for an eigenvector with an eigenvalue between $\\lambda_2$ and $\\lambda_3$, we find the last eigenvector, eigenvalue pair" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "elist4, v4 = invpowermethodshift(A, v0, (l2 + l3)/2, 30)\n", "l4 = elist4[-1]" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ True, True, True, True])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.isclose(A @ v4, l4*v4)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.19" } }, "nbformat": 4, "nbformat_minor": 2 }