From 1d29c7262ed6bdef3269bc935467a4b5bd1fb472 Mon Sep 17 00:00:00 2001
From: Frederic Weidling <fweidli@client61.num.math.uni-goettingen.de>
Date: Mon, 23 Jan 2017 12:29:53 +0100
Subject: [PATCH] blatt11

---
 11_Vektoriterationen.ipynb | 252 +++++++++++++++++++++++++++++++++++++
 1 file changed, 252 insertions(+)
 create mode 100644 11_Vektoriterationen.ipynb

diff --git a/11_Vektoriterationen.ipynb b/11_Vektoriterationen.ipynb
new file mode 100644
index 0000000..848fa24
--- /dev/null
+++ b/11_Vektoriterationen.ipynb
@@ -0,0 +1,252 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Aufgabe 43\n",
+    "## Teil 1: Vektoriteration\n",
+    "Zunächst die Implementation der Vektoriteration. Da wir kein explizites Stoppkriterium in der Vorlesung besprochen haben (wobei das heuristische Stoppkriterium vom Fixpunktsatz von Banach anwendbar ist!) verwenden wir eine maximale Iterationszahl als Endbedingung. Als Rückgabe verwenden wir eine Liste, die die Folge der Eigenwerte enthält.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "from numpy.linalg import norm\n",
+    "\n",
+    "def vektoriteration(A,x0,max_It=10):\n",
+    "    eigenvalue_list=[]\n",
+    "    eigenvalue=norm(x0)\n",
+    "    for j in range(max_It):\n",
+    "        x0=(A@x0)/eigenvalue\n",
+    "        eigenvalue=norm(x0)\n",
+    "        eigenvalue_list.append(eigenvalue)\n",
+    "    return eigenvalue_list,x0"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Nun für die Poissonmatrix. Wir starten mit dem Vektor $(0,0,1)$. Als Erinnerung: der maximale Eigenwert der Matrix ist $2+\\sqrt{2}$ mit Eigenvektor $(1/2,-1/\\sqrt{2},1/2)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Folge der Eigenwerte: [2.2360679774997898, 2.8982753492378879, 3.1922525261132124, 3.3292810259172643, 3.3838222398784663, 3.4036279699578267, 3.4105622571268381, 3.4129583893366058, 3.4137825902599368, 3.4140656449140097]\n",
+      "Bestimmter Eigenvektor: [ 1.69075164 -2.41405427  1.72323681]\n"
+     ]
+    }
+   ],
+   "source": [
+    "from numpy import array\n",
+    "\n",
+    "A=array([[2.0,-1,0],[-1,2,-1],[0,-1,2]])\n",
+    "x0=array([0,0,1.0])\n",
+    "\n",
+    "a,b=vektoriteration(A,x0)\n",
+    "\n",
+    "print('Folge der Eigenwerte: {0}\\nBestimmter Eigenvektor: {1}'.format(a,b))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Da wir lineare Konvergenz des Verfahrens bewiesen haben (Satz 6.11) erwarten wir, dass die Paare $(n,\\log(|\\lambda_n-\\lambda_{1}|))$ etwa auf einer Geraden liegen. Der Wert den wir für die Steigung erwarten ist $\\log(\\lambda_2/\\lambda_1)=\\log(2/(2+\\sqrt{2}))\\approx-0.53$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Die berechnete Steigung ist -1.01109209482217\n"
+     ]
+    }
+   ],
+   "source": [
+    "from numpy import log,arange,sqrt\n",
+    "from scipy.stats import linregress\n",
+    "\n",
+    "a,b=vektoriteration(A,x0)\n",
+    "s1=arange(10)+1\n",
+    "t1=log(abs(array(a)-(2+sqrt(2))))\n",
+    "\n",
+    "reg=linregress(s1,t1)\n",
+    "\n",
+    "print('Die berechnete Steigung ist {}'.format(reg[0]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Wir sehen, dass in diesem Beispiel eine schnellere Konvergenz erhalten. Dies liegt unter anderem daran, dass der Startvektor auch einen Anteil in Richtung des dritten Eigenvektors hat."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Teil 2: Inverse Vektoriteration"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "from numpy.linalg import solve\n",
+    "from numpy import eye\n",
+    "\n",
+    "def inverse_vektoriteration(A,x0,ew,max_It=10):\n",
+    "    I=eye(len(x0))\n",
+    "    for j in arange(max_It):\n",
+    "        x0=solve(A-ew*I,x0)\n",
+    "        x0=x0/norm(x0)\n",
+    "        ew=x0@(A@x0)\n",
+    "    return ew,x0\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Wir erhalten den Eigenwert 3.414213562373095 mit Eigenvektor [-0.5         0.70710678 -0.5       ]\n"
+     ]
+    }
+   ],
+   "source": [
+    "ew,ev=inverse_vektoriteration(A,x0,4)\n",
+    "print('Wir erhalten den Eigenwert {0} mit Eigenvektor {1}'.format(ew,ev))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "ename": "LinAlgError",
+     "evalue": "Singular matrix",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-6-bd12161260ef>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mew\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minverse_vektoriteration\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Wir erhalten den Eigenwert {0} mit Eigenvektor {1}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mew\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mev\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m<ipython-input-4-70e2dd1da485>\u001b[0m in \u001b[0;36minverse_vektoriteration\u001b[0;34m(A, x0, ew, max_It)\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mI\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax_It\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m         \u001b[0mx0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mew\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mI\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m         \u001b[0mx0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mew\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m/opt/anaconda/lib/python3.5/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m    382\u001b[0m     \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'DD->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'dd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    383\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 384\u001b[0;31m     \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    385\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    386\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m/opt/anaconda/lib/python3.5/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m     88\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     91\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     92\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
+     ]
+    }
+   ],
+   "source": [
+    "ew,ev=inverse_vektoriteration(A,x0,2)\n",
+    "print('Wir erhalten den Eigenwert {0} mit Eigenvektor {1}'.format(ew,ev))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Da wir hier bereits mit einem Eigenwert starten kann das Gleichungssystem nicht gelöst werden, da die Matrix die invertiert werden müsste singulär ist."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Wir erhalten den Eigenwert 2.0000000000000004 mit Eigenvektor [-0.70710678  0.          0.70710678]\n"
+     ]
+    }
+   ],
+   "source": [
+    "ew,ev=inverse_vektoriteration(A,x0,1.6)\n",
+    "print('Wir erhalten den Eigenwert {0} mit Eigenvektor {1}'.format(ew,ev))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Wir erhalten den Eigenwert 0.5857864376269049 mit Eigenvektor [-0.5        -0.70710678 -0.5       ]\n"
+     ]
+    }
+   ],
+   "source": [
+    "ew,ev=inverse_vektoriteration(A,x0,1)\n",
+    "print('Wir erhalten den Eigenwert {0} mit Eigenvektor {1}'.format(ew,ev))"
+   ]
+  }
+ ],
+ "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.5.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
-- 
GitLab