From 26c561d8ea0293fab81a77f62b63a341a21496aa Mon Sep 17 00:00:00 2001
From: Frederic Weidling <fweidli@client61.num.math.uni-goettingen.de>
Date: Mon, 7 Nov 2016 08:16:09 +0100
Subject: [PATCH] Blatt 2 Loesung
---
02_LU-Zerlegung_2dPoisson.ipynb | 336 ++++++++++++++++++++++++++++++++
1 file changed, 336 insertions(+)
create mode 100644 02_LU-Zerlegung_2dPoisson.ipynb
diff --git a/02_LU-Zerlegung_2dPoisson.ipynb b/02_LU-Zerlegung_2dPoisson.ipynb
new file mode 100644
index 0000000..6b570bc
--- /dev/null
+++ b/02_LU-Zerlegung_2dPoisson.ipynb
@@ -0,0 +1,336 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Aufgabe 8\n",
+ "## Aufgabe 8.1\n",
+ "__Aufgabe:__ Implementation der LU Zerlegung mit Spaltenpivotsuche.\n",
+ "\n",
+ "Zunächst laden der nötigen Befehle"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import array,arange,triu,tril,eye"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Es wird die Variante mit minimaler Schleifennzahl, also nur einer, gezeigt. Zusätzlich sind noch zwei Sicherungsmechanismen eingebaut, die garantieren, dass jeder Schritt im Algorithmus durchführbar ist."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "def lu_pivot(A):\n",
+ " # Input:\n",
+ " # - A: Matrix von der die LU Zerlegung bestimt werden soll\n",
+ " # Output:\n",
+ " # - A: Matrix, die die LU-Zerlegung enthält\n",
+ " # - p: Vektor der die Pivotisierung speichert\n",
+ " \n",
+ " n,m=A.shape\n",
+ " # Teste ob A quadratisch ist, wenn nicht Abbruch des Algorithmus\n",
+ " if n!=m:\n",
+ " print('Matrix nicht quadratisch')\n",
+ " return None\n",
+ " \n",
+ " p=arange(n)\n",
+ " for kkk in range(n):\n",
+ " restSpalte=abs(A[kkk:,kkk])\n",
+ " # Teste ob der näcshte Schritt durchführbar ist. Aus numerischer Sicht wäre es sinnvoller zu testen, ob\n",
+ " # restSpalte.max()<EPS, wobei EPS ein kleine vorgegebene Zahl ist.\n",
+ " if restSpalte.max()==0:\n",
+ " print('0 auf Diagonale, Algorithmus nicht durchführbar')\n",
+ " return None\n",
+ " \n",
+ " ind=restSpalte.argmax()\n",
+ " # Teste ob Pivotisiert werden muss\n",
+ " if ind!=0:\n",
+ " # Führe Pivotisierung aus, diese ist ohne Schleife möglich, indem man die beiden Zeilen komplett in \n",
+ " # einem Schritt tauscht.\n",
+ " p[[kkk,kkk+ind]]=p[[kkk+ind,kkk]] \n",
+ " A[[kkk,kkk+ind],:]=A[[kkk+ind,kkk],:] \n",
+ " \n",
+ " # Durchführung des nächsten Eliminationsschrittes, auch hier ohne Schleife. Der Trick ist den additiven Term\n",
+ " # in der Berechnung der a_ij für i,j>k durch das passende Matrixprodukt darzustellen.\n",
+ " l=A[kkk+1:,kkk]/A[kkk,kkk]; \n",
+ " A[kkk+1:,kkk]=l; \n",
+ " # Korrektur der Terme a_ij mit i,j>kkk. Dabei muss daran gedacht werden die Vektoren l und A[kkk,kkk+1:] \n",
+ " # explizit als Spalten- bzw. Zeilenvektor darzustellen. Tut man dies nicht, verwendet python punktweise\n",
+ " # Multiplikation und Broadcasting für die Rechnung (wirft also keinen Fehler aus), allerdings erhält man\n",
+ " # nicht das gewünschte Ergebnis.\n",
+ " A[kkk+1:,kkk+1:]=A[kkk+1:,kkk+1:]-l[:,None]@A[kkk,kkk+1:].reshape(1,-1); \n",
+ " \n",
+ " return A,p"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Die Ausgabe besteht nun aus der Überschriebenden Matrix A, die (L-I)+U enthält. Möchte man aus dieser die Matrix A rekonstruieren muss man zunächst L\\*U berechnen und erhält dann die Zeilen von A in der in p gespeicherten Reihenfolge auslesen. Hier ein paar Beispiele: "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[0 1 2]\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 2., -1., 0.],\n",
+ " [-1., 2., -1.],\n",
+ " [ 0., -1., 2.]])"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "A=array([[2,-1,0],[-1,2,-1],[0,-1,2]])*1.0\n",
+ "B,p=lu_pivot(A)\n",
+ "print(p)\n",
+ "(tril(B,-1)+eye(3))@triu(B)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Hier wurde wir nicht getauschen (zu erkennen an p=[0,1,2]). Wichtig ist, dass wir A als array von Kommazahlen übergeben, sonst werden Rechenfehler auftauchen."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[2 0 1]\n",
+ "[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
+ " [ 0.00000000e+00 0.00000000e+00 0.00000000e+00]\n",
+ " [ 0.00000000e+00 0.00000000e+00 -4.44089210e-16]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "A=array([[1,7,4],[1,-1,4],[8,2,3]])*1.0\n",
+ "B,p=lu_pivot(A.copy())\n",
+ "print(p)\n",
+ "C=(tril(B,-1)+eye(3))@triu(B)\n",
+ "print(C-A[p,:])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Wir sehen, dass wir bis auf (unvermeidbare) numerische Rechenfehler das erwartete Ergebnis erhalten.\n",
+ "## Aufgabe 8.1\n",
+ "__Aufgabe:__ Erstellen der 2d-Poissonmatrix bei Schrittweise 1/n auf dem Einheitsquadrat und Betrachtung der relativen Entwicklung der Anzahl der Nulleinträge von der LU-Zerlgung der Matrix.\n",
+ "\n",
+ "Wir kopieren zunächst einmal die Funktion zum Erstellen von Tridiagonalmatrizen von letzter Woche:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "from numpy import zeros, count_nonzero\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "def tridiag(a,b,c,n):\n",
+ " A=zeros([n,n]) \n",
+ " for iii in range(n):\n",
+ " if iii==0:\n",
+ " A[iii,iii]=b \n",
+ " A[iii,iii+1]=c\n",
+ " elif iii==n-1:\n",
+ " A[iii,iii]=b\n",
+ " A[iii,iii-1]=a\n",
+ " else:\n",
+ " A[iii,iii]=b\n",
+ " A[iii,iii-1]=a\n",
+ " A[iii,iii+1]=c\n",
+ " return A"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Diese werden wir nutzen um die Matrizen zu erstellen, die auf der Diagonale der 2d-Poissonmatrix sitzen."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "def poisson2d_matrix(n):\n",
+ " A=zeros([(n-1)**2,(n-1)**2]) # Hier speichern wie die Poissonmatrix rein\n",
+ " B=tridiag(-1,4,-1,n-1) # Die Matrix entlang der Diagonalen\n",
+ " C=-eye(n-1) # Die Matrix auf den Matrixnebendiagonalen\n",
+ " # Jetzt bauen wir Matrix nach dem selben Prinzip zusammen, wie in tridiag:\n",
+ " for iii in range(n-1): # für \"blockweise\" durchlaufen der Matrix\n",
+ " ind=iii*(n-1) # Umrechnen von Blockindex zu Matrixindex\n",
+ " A[ind:ind+(n-1),ind:ind+n-1]=B #Setzen des Diagonalblocks\n",
+ " # Setzen der Nebendiagonalblöcke\n",
+ " if iii==0:\n",
+ " A[ind:ind+(n-1),ind+(n-1):ind+2*(n-1)]=C\n",
+ " elif iii==n-2:\n",
+ " A[ind:ind+(n-1),ind-(n-1):ind]=C\n",
+ " else:\n",
+ " A[ind:ind+(n-1),ind+(n-1):ind+2*(n-1)]=C\n",
+ " A[ind:ind+(n-1),ind-(n-1):ind]=C\n",
+ " return n**2*A"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Test das die richtige Matrix bei rauskommt mit n=4 (sollte also eine $(4-1)^2\\cdot(4-1)^2$ Matrix sein)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 64., -16., 0., -16., -0., -0., 0., 0., 0.],\n",
+ " [-16., 64., -16., -0., -16., -0., 0., 0., 0.],\n",
+ " [ 0., -16., 64., -0., -0., -16., 0., 0., 0.],\n",
+ " [-16., -0., -0., 64., -16., 0., -16., -0., -0.],\n",
+ " [ -0., -16., -0., -16., 64., -16., -0., -16., -0.],\n",
+ " [ -0., -0., -16., 0., -16., 64., -0., -0., -16.],\n",
+ " [ 0., 0., 0., -16., -0., -0., 64., -16., 0.],\n",
+ " [ 0., 0., 0., -0., -16., -0., -16., 64., -16.],\n",
+ " [ 0., 0., 0., -0., -0., -16., 0., -16., 64.]])"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "poisson2d_matrix(4)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Jetzt: Verhältnis der nicht Nulleinträge bestimmen, da der jupyter Server nicht so schnell ist nur bis n=51"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAFkCAYAAACw3EhvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmYXFWd//H3lx2CBJUZcEFHRBGECXQzaIQQNtlkdW9k\n34lhaWQVfwyKI8oiIDuiso3tgLJqCCCEHWToZksCCiIEEghCoEMWBJLz++NUTzptJ53uruqqW/V+\nPU8eqVu3636Pp6C/+dxz742UEpIkSUtVuwBJklQbbAokSRJgUyBJkkpsCiRJEmBTIEmSSmwKJEkS\nYFMgSZJKbAokSRJgUyBJkkpsCiRJEtDPpiAiToyIhyNiZkRMj4jrI+LT3d5fJiJ+EhFPRMSsiJga\nEVdExIfKX7okSSqn/iYFo4DzgM8B2wDLArdFxIql91cCNgS+D2wE7A6sA9xYlmolSVLFxGAeiBQR\nqwGvApunlO5bxD4bA38CPp5SemnAB5MkSRU12DUFqwIJmLEE+7w5yGNJkqQKGnBSEBEB3Ay8L6U0\nehH7LA/cD0xOKe29iH0+CGwHPA+8PaBiJElqTCsA/wbcmlJ6fbAftswgfvZCYD1g097ejIhlgGvJ\nKcGYxXzOdsB/D6IOSZIa3beAXw/2QwbUFETE+cCOwKiU0su9vN/VEKwJbJVSmrWYj3se4Oqrr2bd\nddcdSDk1p7W1lbPPPrvaZZRNPY2nnsYCjqeW1dNYwPHUqqeeeoo999wTSr9LB6vfTUGpIdgVGJ1S\nmtLL+10NwVrAlimlN/r4yLcB1l13XZqamvpbTk0aPnx43YwF6ms89TQWcDy1rJ7GAo6nAMpy+r1f\nTUFEXAi0ALsAsyNi9dJbnSmlt0sNwe/IlyXuBCzbbZ8ZKaV3y1G0JEkqv/4mBYeS1wjc1WP7fsCV\nwEfIzQDAY6X/jdLPbAncM6AqJUlSxfWrKUgpLfYSxpTSC8DSg6pIkiRVhc8+qICWlpZql1BW9TSe\nehoLOJ5aVk9jAcfTKAZ1R8OyFBDRBLS3t7fX26IPSZIqqqOjg+bmZoDmlFLHYD/PpECSJAE2BZIk\nqcSmQJIkATYFkiSpxKZAkiQBNgWSJKnEpkCSJAE2BZIkqcSmQJIkATYFkiSpxKZAkiQBNgWSJKnE\npkCSJAE2BZIkqcSmQJIkATYFkiSpxKZAkiQBNgWSJKnEpkCSJAE2BZIkqcSmQJIkAf1sCiLixIh4\nOCJmRsT0iLg+Ij7dy34/iIhpETEnIm6PiLXLV7IkSaqE/iYFo4DzgM8B2wDLArdFxIpdO0TE8cBY\n4GBgE2A2cGtELFeWiiVJUkUs05+dU0o7dn8dEfsCrwLNwH2lzUcCp6aUfl/aZ29gOrAbcM0g65Uk\naUjMnQuzZ8Nqq1W7kqEz2DUFqwIJmAEQEZ8A1gDu6NohpTQT+BMwcpDHkiRpSDz0EGy0ERx8cLUr\nGVoDbgoiIoBzgPtSSpNLm9cgNwnTe+w+vfSeJEk1a+5cOPZY2HRTGD4cfvjDalc0tPp1+qCHC4H1\ngE3LVIskSVXz0EOw777wt7/Bj34E3/kOLDOY35IFNKDhRsT5wI7AqJTSy93eegUIYHUWTgtWBx5d\n3Ge2trYyfPjwhba1tLTQ0tIykBIlSVoic+fCySfDT38KG28Mjz4K661X7ar+WVtbG21tbQtt6+zs\nLOsxIqXUvx/IDcGuwOiU0nO9vD8NOCOldHbp9SrkBmHvlNK1vezfBLS3t7fT1NQ0gCFIkjQwDz4I\n++2X04Ef/KB46UBHRwfNzc0AzSmljsF+Xr+GHhEXAi3ALsDsiFi99FZnSunt0j+fA3wvIp4FngdO\nBV4CbhxssZIklUNR0oGh1t9+6FDyQsK7emzfD7gSIKV0ekSsBFxCvjrhXmCHlNI7gytVkqTB60oH\nnn8eTjsNjj66WOlAJfX3PgVLdLVCSukU4JQB1CNJUkX0lg6su261q6ot9kaSpLpnOrBkfCCSJKlu\nzZ0LxxyT7zuw6qo5HTjuOBuCRfH/FklSXXrggZwOvPAC/OQnOR1YeulqV1XbTAokSXWlKx3YbDN4\n//tzOnDssTYES8KkQJJUN0wHBsekQJJUeHPn5hsPmQ4MjkmBJKnQTAfKx6RAklRI3dOBD3wAHnvM\ndGCwTAokSYVz//2w//6mA+VmUiBJKow5c3I6MGqU6UAlmBRIkgrh/vvz2oEpU+D006G11Wag3EwK\nJEk1rXs68MEP5nTgmGNsCCrBpECSVLNMB4aWSYEkqebMmZMXD44aBautZjowVEwKJEk1pSsdePFF\nOOMMOOoom4GhYlIgSaoJvaUD3/mODcFQMimQJFXdfffl+w6YDlSXSYEkqWrmzMmLBzff3HSgFpgU\nSJKqwnSg9pgUSJKGlOlA7TIpkCQNGdOB2mZSIEmqONOBYjApkCRV1H335fsOvPQSnHkmHHmkzUCt\nMimQJFVE93TgX/81pwM+4ri29bspiIhREXFTREyNiPkRsUuP94dFxPkR8WJEzImISRFxSPlKliTV\nuvvugxEj4OKLczpwzz2wzjrVrkp9GUhSMAx4DBgDpF7ePxvYFtgD+Ezp9fkRsdNAi5QkFcOcOXnx\noOlAMfV7TUFKaTwwHiAiopddRgJXpJTuLb2+LCIOBTYBfj/QQiVJte3ee/PagalT4ayz4IgjbAaK\nphJrCh4AdomIDwNExJbAp4BbK3AsSVKVdaUDo0fD6qvD44/7iOOiqsTVB4cDlwIvRcR7wDzgoJTS\n/RU4liSpikwH6kslmoIjgM8BOwFTgM2BCyNiWkrpzkX9UGtrK8OHD19oW0tLCy0tLRUoUZI0GLNn\nw0knwc9+BiNHwrhx8OlPV7uq+tbW1kZbW9tC2zo7O8t6jEipt7WCS/jDEfOB3VJKN5VerwB0lrbd\n0m2/nwMfSSnt2MtnNAHt7e3tNDU1DbgWSdLQ6J4O/OhHpgPV1NHRQXNzM0BzSqljsJ9X7jUFy5b+\nzOuxfV4FjiVJGkKzZ+cbD40eDWusAU884dqBetPv0wcRMQxYG+i68mCtiBgBzEgpvRgRdwNnRsTh\nwAvAFsDewFHlKVmSNNRcO9AYBvK3942BR4F28n0KzgI6gO+X3v8G8L/A1cAk4DjgxJTSpYOuVpI0\npEwHGstA7lNwN4tpJlJKrwIHDKYoSVL13XNPfqLhtGnw05/C4YfbDNQ7z/NLkhbSMx14/HEfcdwo\nfEqiJOn/dE8Hzj7bdKDRmBRIkpg9Oy8eNB1obCYFktTg7r47pwMvvwznnANjx9oMNCqTAklqUF3p\nwBZbwIc/nNOBI4+0IWhkJgWS1IB6pgOHHw5L+dfEhudXQJIayKLSARsCgUmBJDUM0wH1xa+DJNW5\n2bNzA2A6oL6YFEhSHTMdUH/41ZCkOtQ9HfjIR/IzC0wH1BeTAkmqM3fdBQccYDqg/vNrIkl1Ytas\n3ABsuaXpgAbGpECS6sBdd+W1A6+8Aueem+9KaDOg/vIrI0kF1j0d+OhHczpwxBE2BBoYkwJJKijT\nAZWbXx9JKphZs3IDsOWWsOaapgMqH5MCSSqQrnRg+nT42c/g29+2GVD5+FWSpALoLR3wUkOVm0mB\nJNW4CRPyfQdMB1Rpfq0kqUbNmpUbgK22Mh3Q0DApkKQaZDqgavArJkk1xHRA1dTvr1lEjIqImyJi\nakTMj4hdetln3Yi4MSLejIhZEfGniPhoeUqWpPo0YQJssAFcfnlOByZMgE9+stpVqZEMpPccBjwG\njAFSzzcj4pPAvcBkYHNgA+BU4O2BlylJ9at7OvDxj5sOqHr6vaYgpTQeGA8QEdHLLj8E/pBSOrHb\ntr8NrDxJqm8TJuT7Drz6Kpx3HowZYzOg6inrV6/UJHwJeCYixkfE9Ih4KCJ2LedxJKnoeksHvE2x\nqq3cX79/BVYGjgfGAV8Ergeui4hRZT6WJBXSnXcuWDtw3nn5tWsHVAvKfUliV5NxQ0rpZ6V/fiIi\nvgAcSl5rIEkNadYsOO44uOgiGD0a7rgD1lqr2lVJC5S7KXgNeA94qsf2p4BNF/eDra2tDB8+fKFt\nLS0ttLS0lLVASaqGO+/M9x149VU4/3w47DBPFah/2traaGtrW2hbZ2dnWY8RKf3TBQRL/sMR84Hd\nUko3ddt2P/BsSmmfbtuuA+aklPbs5TOagPb29naampoGXIsk1aK33srpwMUX53Tgl780HVD5dHR0\n0NzcDNCcUuoY7Of1OymIiGHA2kDXlQdrRcQIYEZK6UXgDOA3EXEvMAHYAdgJGD3YYiWpSO68M19Z\n8NprpgMqhoF8PTcGHgXayfcpOAvoAL4PkFK6gbx+4DjgCWB/4MsppQfLUbAk1bq33soNwNZbwyc+\nka8s8DbFKoKB3KfgbvpoJlJKlwOXD6wkSSquO+7Iawdeew0uuAAOPdRmQMXhV1WSyqArHdhmm5wO\nPPmkNyJS8fiUREkaJNMB1Qu/tpI0QG+9lRuAbbbJVxSYDqjoTAokaQD++MecDrz+Olx4IRxyiM2A\nis+vsCT1w8yZuQH44hfzrYmffNJLDVU/TAokaQmZDqje+XWWpD50TwfWXhsmTjQdUH0yKZCkxbj9\ndjjwQJgxIz/I6OCDbQZUv/xqS1IvutKBbbfN6cCTT3qpoeqfSYEk9dAzHTjkEIjo++ekorPnlaSS\nmTPz6YFtt4VPfWpBOmBDoEZhUiBJwG235XTgjTfyY44PPthmQI3HpEBSQ5s5Ew46CLbbDj796ZwO\neLpAjcqkQFLDMh2QFmZSIKnhdHYuSAfWWSffd8B0QDIpkNRgbr01pwNvvgmXXJKbA5sBKTMpkNQQ\nutKB7beHz3wmpwOeLpAWZlIgqe6ZDkhLxqRAUt3q7MzNwPbbw7rrmg5IfTEpkFSXxo/PiUBnJ1x6\naW4ObAakxTMpkFRXOjvz44132GFBOuDpAmnJmBRIqhumA9LgmBRIKrw33zQdkMrBpEBSod1yS24A\nZs6En/88Nwc2A9LA9DspiIhREXFTREyNiPkRscti9r24tM8RgytTkhbWlQ7suCN89rM5HfB0gTQ4\nAzl9MAx4DBgDpEXtFBG7A58Dpg6sNEnq3S23wPrrw7XX5nRg/Hj42MeqXZVUfP1uClJK41NKJ6eU\nbgR67ckj4iPAucAewHuDK1GSsjffhP33z+nA+uubDkjlVvY1BRERwJXA6Smlp8J/WyWVwbhx+cZD\nb70Fl12WmwP/8yKVVyWuPjgBeCeldH4FPltSg3nzTdhvP/jSlxakAy4mlCqjrElBRDQDRwAb9fdn\nW1tbGT58+ELbWlpaaGlpKVN1korGdEBaoK2tjba2toW2dXZ2lvUYkdIi1wr2/cMR84HdUko3lV4f\nCZzFwgsQlwbmA1NSSmv18hlNQHt7eztNTU0DrkVS/XjzTWhthcsvh+22y4sJ11yz2lVJtaejo4Pm\n5maA5pRSx2A/r9xrCq4Ebu+x7bbS9l+V+ViS6tC4cfm+A7NmwS9+kU8dmA5IQ6PfTUFEDAPWZsGV\nB2tFxAhgRkrpReCNHvu/C7ySUnpmsMVKql9vvJHTgSuuyE81vPRS0wFpqA0kKdgYmEA+RZDIpwsA\nrgD272X/gZ+fkNQQ/vCHvHbAdECqrn43BSmlu+nHVQu9rSOQJPjndODnP4ePfrTaVUmNy2cfSKqK\n7unAL38J++5rOiBVm09JlDSk3ngD9tkHdtoJRoyASZM8XSDVCpMCSUOmKx2YPdt0QKpFJgWSKq5n\nOjBxoumAVItMCiRV1O9/n9OBOXPgV7/KzYHNgFSbTAokVcQbb8Dee8POO8NGG+W1A54ukGqbSYGk\nsrv5ZjjkENMBqWhMCiSVTVc6sMsupgNSEZkUSCqL7unA5Zfn5sBmQCoWkwJJgzJjBuy1V04Hmppy\nOuDpAqmYTAokDdjNN+crC+bONR2Q6oFJgaR+654ONDebDkj1wqRAUr/cdFNeO/D22/lBRnvtZTMg\n1QuTAklLZMYM2HNP2HXXBemApwuk+mJSIKlPN96Y04F//AOuvDI3BzYDUv0xKZC0SF3pwG67wX/8\nR04HPF0g1S+TAkm9Mh2QGo9JgaSFvP46fOtbpgNSIzIpkPR/TAekxmZSIGmhdGCTTWDyZNMBqRGZ\nFEgN7oYb4NBD4Z134KqrcnNgMyA1JpMCqUG9/jrssQfsvntOByZN8nSB1OhMCqQGdP31OR149124\n+urcHNgMSOp3UhARoyLipoiYGhHzI2KXbu8tExE/iYgnImJWaZ8rIuJD5S1b0kB0pQNf/jJ8/vM5\nHfB0gaQuAzl9MAx4DBgDpB7vrQRsCHwf2AjYHVgHuHEQNUoqg+uvh/XWg/Hjczpwww3wIdt1Sd30\n+/RBSmk8MB4gYuG/X6SUZgLbdd8WEWOBP0XER1NKLw2iVkkD8PrrMHYs/OY3+amGF19sMyCpd0Ox\npmBVcqLw5hAcS1I3110Hhx3m2gFJS6aiVx9ExPLAj4Ffp5RmVfJYkhZ47TX45jfhK1+BkSPzfQdc\nOyCpLxVLCiJiGeBackowplLHkbSw3/0upwPz5sGvf52bA5sBSUuiIk1Bt4ZgTWCrJUkJWltbGT58\n+ELbWlpaaGlpqUSJUt35+9/z2oFrrsl3JrzoIlhjjWpXJalc2traaGtrW2hbZ2dnWY8RKfW8gKAf\nPxwxH9gtpXRTt21dDcFawJYppRl9fEYT0N7e3k5TU9OAa5Ea2W9/C2PG5HTg/PNNB6RG0dHRQXNz\nM0BzSqljsJ/X76QgIoYBawNd/8lZKyJGADOAl4HfkS9L3AlYNiJWL+03I6X07mALlrTA3/8O3/42\nXHttvjPhRRfB6qv3/XOS1JuBnD7YGJhAXiuQgLNK268g359g59L2x0rbo/R6S+CewRQraYFrr83p\nQEr5csOvf910QNLgDOQ+BXez+KsWfJ6CVEGvvprTgd/+Nt+Z8MILTQcklYfPPpAK5JprckMA8D//\nA1/7mumApPLxb/VSAbz6am4AvvENGD06P7PA0wWSys2kQKphKS1IByLyP3/ta9WuSlK9MimQatT0\n6fDVr+bLC7faKqcDNgSSKsmkQKoxKeX1AmPHwlJLmQ5IGjomBVINmT49P6+gpQW23tp0QNLQMimQ\nakDXvQbGjoWll873IPjqV6tdlaRGY1IgVdkrr+T7DeyxB3zxizkdsCGQVA0mBVKVpARtbXD44bDM\nMvlmRF/5SrWrktTITAqkKnjllfysgm99C7bdNqcDNgSSqs2kQBpCKcGvf53TgWWXhd/9Lp86kKRa\nYFIgDZGXX4bddoM994Ttt4fJk20IJNUWkwKpwlKC//5vOOIIWG45uO66fOpAkmqNSYFUQdOmwa67\nwl57wQ475LUDNgSSapVJgVQBKcHVV+d0YPnl4frr86kDSaplJgVSmU2bBrvsAnvvDTvtlNcO2BBI\nKgKTAqlMUoKrroIjj4QVVoAbb8zNgSQVhUmBVAZTp8LOO8M+++R0YNIkGwJJxWNSIA1CSnDFFXDU\nUbDiiqYDkorNpEAaoKlTcyqw3375CgPTAUlFZ1Ig9VP3dGClleDmm3NzIElFZ1Ig9cNLL8GXvrRw\nOmBDIKlemBRISyAl+NWvoLUVVl7ZdEBSfTIpkPrw4ouw445wwAH5WQUTJ9oQSKpP/W4KImJURNwU\nEVMjYn5E/NPSqoj4QURMi4g5EXF7RKxdnnKloZMS/OIXsP768MQT8Ic/5LTg/e+vdmWSVBkDSQqG\nAY8BY4DU882IOB4YCxwMbALMBm6NiOUGUac0pF58MT+r4MAD4StfyWsHdtyx2lVJUmX1e01BSmk8\nMB4gIqKXXY4ETk0p/b60z97AdGA34JqBlypVXlc6cPTRsMoqMG5cbg4kqRGUdU1BRHwCWAO4o2tb\nSmkm8CdgZDmPJZXblCmw/fZw0EHwta/ltQM2BJIaSbkXGq5BPqUwvcf26aX3pJqTEvz853ntwOTJ\ncMstOS1YddVqVyZJQ6tmLklsbW1l+PDhC21raWmhpaWlShWpEUyZktcN3H57vrrgrLOgx9dQkmpC\nW1sbbW1tC23r7Ows6zEipX9aK7jkPxwxH9gtpXRT6fUngL8CG6aUnui2313Aoyml1l4+owlob29v\np6mpacC1SP3RlQ4cc0xuAi67DLbbrtpVSVL/dHR00NzcDNCcUuoY7OeV9fRBSulvwCvA1l3bImIV\n4HPAA+U8ljRQL7wA224LhxwC3/hGXjtgQyBJAzh9EBHDgLWBrisP1oqIEcCMlNKLwDnA9yLiWeB5\n4FTgJeDGslQsDVBKcOmlOR14//vh1ltzcyBJygaypmBjYAJ5QWECziptvwLYP6V0ekSsBFwCrArc\nC+yQUnqnDPVKA/L883ntwB135KsLzjjDtQOS1NNA7lNwN32cdkgpnQKcMrCSpPJJCS65BI491nRA\nkvrisw9Ut55/HrbZBg47DPbYI68dsCGQpEWzKVDdmT8fLroo33fg2Wfz5YaXXJLvUChJWjSbAtWV\nv/0tpwNjxsCee+Z0YJttql2VJBWDTYHqwvz5cOGFsMEG8NxzOR24+GJ43/uqXZkkFYdNgQqvKx34\n9rdhr73gySdNByRpIGwKVFjz58MFFyxIB/74x7yWwHRAkgbGpkCF9NxzsNVWMHYs7LNPTge23rrv\nn5MkLZpNgQpl/nw4//ycDrzwQr4Z0QUXmA5IUjnYFKgw/vrXnA4cfjjsu29OB7baqtpVSVL9sClQ\nzZs/H847D/793/Ojju+8M6cDK69c7cokqb7YFKimPfssbLklHHEE7L8/PPFEfi1JKj+bAtWk+fPh\nZz/L6cBLL8GECTktMB2QpMqxKVDNefZZ2GILOPJIOOCAnA5ssUW1q5Kk+mdToJoxfz6ce25OB6ZO\nhbvuyunAsGHVrkySGoNNgWrCM8/A6NFw1FFw0EE5HRg9utpVSVJjsSlQVc2bB+ecAyNGwMsvw913\n57TAdECShp5NgarmL3/JacDRR8PBB8Pjj8Pmm1e7KklqXDYFGnLz5sFPf5rTgenTczpwzjmmA5JU\nbTYFGlJ/+UtOA445Bg49NKcDo0ZVuypJEtgUaIh0TwdefRXuuQfOPhtWWqnalUmSutgUqOL+/Oec\nBhxzDBx2WE4HNtus2lVJknqyKVDFzJsHZ54JG24Ir70G996b0wLTAUmqTTYFqoinn85pwHHHwZgx\n8NhjsOmm1a5KkrQ4ZW8KImKpiDg1Ip6LiDkR8WxEfK/cx1FtmjcPzjgjpwMzZsB998FZZ5kOSFIR\nLFOBzzwBOATYG5gMbAxcHhFvppTOr8DxVCOefhr22w/+9Kd874FTT4UVV6x2VZKkJVWJpmAkcGNK\naXzp9ZSI2APYpALHUg147728VuDkk+HjH8/pwBe+UO2qJEn9VYk1BQ8AW0fEpwAiYgSwKTCuAsdS\nlU2enNcKnHgijB2b1w7YEEhSMVUiKfgxsArwdETMIzceJ6WUflOBY6lK3nsvrx045RRYay24/374\n/OerXZUkaTAq0RR8A9gD+CZ5TcGGwLkRMS2ldFUFjqch9uSTee3Ao4/CscfmxmCFFapdlSRpsCrR\nFJwOnJZSurb0elJE/BtwIrDIpqC1tZXhw4cvtK2lpYWWlpYKlKiBePdd+MlP4Ac/gE99Ch58EDZx\npYgkDYm2tjba2toW2tbZ2VnWY0RKqbwfGPEa8N2U0qXdtp0I7JNS+kwv+zcB7e3t7TQ1NZW1FpXP\n44/ndOCJJ+D44/OiwuWXr3ZVktTYOjo6aG5uBmhOKXUM9vMqkRTcDHwvIl4CJgFNQCtwWQWOpQp7\n5x047TT44Q/hM5/Jlxvm758kqd5UoikYC5wKXAD8KzANuKi0TQXy6KOw7775CoPvfhdOOgmWW67a\nVUmSKqXsTUFKaTZwdOmPCugf/8jJwGmnwfrrw8MPw0YbVbsqSVKlVSIpUIE98kheO/D003ndwAkn\nmA5IUqPwgUgCcjrw3e/mew0su2xuDk4+2YZAkhqJSYF4+OGcDjzzDHz/+/nJhssuW+2qJElDzaSg\ngb39dr68cOTI/OCijo68mNCGQJIak0lBg3rwQdh/f3juOfiv/4JjjoFl/DZIUkMzKWgwc+fmBmDT\nTWGVVfJlhyecYEMgSTIpaCj33ZfTgSlT8u2KW1ttBiRJC5gUNIA5c3IDsPnmsNpq+fHGxx5rQyBJ\nWpi/FurcPffkdGDqVDjzTDjySFh66WpXJUmqRSYFdWrWLDj8cBg9Gj70ofxAo6OPtiGQJC2aSUEd\nmjABDjgAXnkFzj0Xxo6FpWz/JEl98FdFHXnrLRgzBrbaCtZcMz/m+IgjbAgkSUvGpKBO3HFHTgde\new3OOy83BzYDkqT+8NdGwc2cCYccAttsA2utldMBTxdIkgbCpKDAbrsNDjwQ3ngDLroIDj7YZkCS\nNHD+Cimgzs7cDGy3HayzDkycCIceakMgSRock4KCGTcuJwIzZ8Kll+bmIKLaVUmS6oF/tyyIN97I\njzf+0pdg/fVzOnDQQTYEkqTyMSkogJtvzosJZ8+GX/wiNwc2A5KkcjMpqGEzZsDee8Muu8BGG8Gk\nSfmWxTYEkqRKMCmoUTfcAIcdBm+/DZdfnpsDmwFJUiWZFNSY116DPfaA3XeHjTfO6cA++9gQSJIq\nz6Sghlx3XU4H3n0XrroKvvUtmwFJ0tAxKagBf/87fOMb8JWvwBe+AJMnw5572hBIkoZWRZqCiPhw\nRFwVEa9FxJyIeDwimipxrKK79lpYb7387IK2tpwWrLFGtauSJDWisjcFEbEqcD/wD2A7YF3gO8Ab\n5T5WkU2fDl/9Knz96zB6dF478M1vmg5IkqqnEmsKTgCmpJQO7LbthQocp5BSgt/8Bg4/PDcA11wD\nX/tatauSJKkypw92Bh6JiGsiYnpEdETEgX3+VAN45RX48pfz1QVbb53XDtgQSJJqRSWagrWAw4A/\nA9sCFwE/i4i9KnCsQkgJrr46rx144AH47W/hf/4H/uVfql2ZJEkLVOL0wVLAwyml/1d6/XhErA8c\nCly1qB9qbW1l+PDhC21raWmhpaWlAiUOnWnT8hMMb745JwTnngurrVbtqiRJRdPW1kZbW9tC2zo7\nO8t6jEgciQq0AAAJ60lEQVQplfcDI54HbkspHdxt26HASSmlNXvZvwlob29vp6mpfi5QSAmuvBKO\nOgqWXx4uvhh2263aVUmS6klHRwfNzc0AzSmljsF+XiVOH9wPrNNj2zo00GLDl17KTzPcd1/Yeee8\ndsCGQJJU6yrRFJwNfD4iToyIT0bEHsCBwPkVOFZNSSk/xfCzn4XHHoObbsppwQc+UO3KJEnqW9mb\ngpTSI8DuQAvwJHAScGRK6TflPlYtmTIFtt8eDjwwX2EwaVJOCSRJKoqKPPsgpTQOGFeJz641KcFl\nl8F3vgOrrALjxsEOO1S7KkmS+s9nHwzC88/DttvCwQfnOxNOmmRDIEkqLpuCAZg/Hy66CDbYAP78\nZxg/PqcFPa6olCSpUGwK+um552CbbWDMmHzfgYkTYbvtql2VJEmDZ1OwhObPh/PPz+nAc8/B7bfD\nJZfkdQSSJNUDm4Il8Ne/wpZb5ocY7bsvPPlkTgskSaonNgWLMX9+vi3xBhvAiy/CnXfCBRfA+95X\n7cokSSo/m4JF+MtfYPPN822KDzgAnngipwWSJNUrm4Ie5s2Ds86CESPyo47vugvOOw9WXrnalUmS\nVFk2Bd08/TRsthkceywccgg8/jiMHl3tqiRJGho2BcB778Hpp8OGG8Lrr8M998A558CwYdWuTJKk\nodPwTcHkybDppnDCCTB2bH6Q0WabVbsqSZKGXsM2Be+9B6edBhttBJ2dcP/9cOaZsNJK1a5MkqTq\naMimYOJEGDkSvve9fHXBo4/m15IkNbKGagrefRd++ENoaoLZs+GBB+AnP4EVV6x2ZZIkVV9FHp1c\nix5/HPbbL99v4Ljj4OSTYYUVql2VJEm1o+6Tgnfege9/HzbeOCcFDz0EP/qRDYEkST3VdVLw6KM5\nHZg4EU48Ma8hWH75alclSVJtqsuk4J138umBTTaBlODhh+HUU20IJElanLpLCh55JKcDTz8NJ50E\n3/0uLLdctauSJKn21U1S8I9/5Abg85+HZZaB//1fOOUUGwJJkpZUXSQFDz+c04FnnsmNwPHHw7LL\nVrsqSZKKpdBJwdtv5wZg5Mh8NUF7e15MaEMgSVL/FbYpePDBfIvic87Jiwgfegg22KDaVWVtbW3V\nLqGs6mk89TQWcDy1rJ7GAo6nUVS8KYiIEyJifkT8tByfN3cuHHNMfojR+94HHR15LUEtpQP19mWr\np/HU01jA8dSyehoLOJ5GUdE1BRHxH8DBwOPl+Lz774f994cXXoAf/xiOPjovKpQkSYNXsaQgIlYG\nrgYOBN4c7OdNmQJbbAEf+EC+KdFxx9kQSJJUTpU8fXABcHNK6c5yfNjHPga33w733QfrrluOT5Qk\nSd1V5O/aEfFNYENg4yXYfQWAp556qs8dV1klP9io1nV2dtLR0VHtMsqmnsZTT2MBx1PL6mks4Hhq\nVbffnWV5ok+klMrxOQs+MOKjwCPANimliaVtE4BHU0pH97L/HsB/l7UISZIay7dSSr8e7IdUoinY\nFbgOmAdEafPSQCptWz51O2hEfBDYDngeeLusxUiSVN9WAP4NuDWl9PpgP6wSTcEw4OM9Nl8OPAX8\nOKXU93kCSZI05Mq+piClNBuY3H1bRMwGXrchkCSpdg3VHQ3LG0dIkqSyK/vpA0mSVEyFffaBJEkq\nL5sCSZIEDGFTEBGjIuKmiJhaekDSLr3s84OImBYRcyLi9ohYe6jq66++xhMRvypt7/5nXLXqXZyI\nODEiHo6ImRExPSKuj4hP97Jfzc/PkoylYHNzaEQ8HhGdpT8PRMT2Pfap+Xnp0td4ijQ3PS3q4W9F\nmp/uehtPkeYnIv6zl1p7LoIvxNz0NZZyzstQJgXDgMeAMfSy8DAijgfGkh+gtAkwG7g1IpYbwhr7\nY7HjKbkFWB1Yo/SnZWhK67dRwHnA54BtgGWB2yJixa4dCjQ/fY6lpChz8yJwPNAENAN3AjdGxLpQ\nqHnpstjxlBRlbv7Poh7+VsD5Afp8mF2R5mciC9e6WdcbBZybRY6lpDzzklIa8j/AfGCXHtumAa3d\nXq8CzAW+Xo0ayzCeXwHXVbu2AY5ntdKYNiv6/CxiLIWdm1L9rwP7FXleFjOews0NsDLwZ2ArYALw\n027vFW5++hhPYeYH+E+gYzHvF2ZulmAsZZuXmlhTEBGfIHc2d3RtSynNBP4EjKxWXWWwRSnCfjoi\nLoyID1S7oCW0Kjn9mAGFn5+FxtJN4eYmIpaK/FyRlYAHCj4v/zSebm8VbW56ffhbgeenr4fZFWl+\nPhX5FO9fI+LqiFgTCjs3vY6lm7LMS608fHgN8n+4p/fYPr30XhHdAvwO+BvwSeA0YFxEjEyl1q4W\nRUQA5wD3pZS6zlkVcn4WMRYo2NxExPrAg+Tbmb4F7J5S+nNEjKSY89LreEpvF21uFvfwt8L9e9PH\neKBY8/MQsC859fgQcApwT+n7V7S56W0s90bEZ1O+YWDZ5qVWmoK6k1K6ptvLSRHxJPBXYAtyJFer\nLgTWAzatdiFl0OtYCjg3TwMjgOHAV4ErI2Lz6pY0KL2OJ6X0dJHmJvLD384hP/zt3WrXM1hLMp4i\nzU9K6dZuLydGxMPAC8DXyd/BwuhjLL8q57zUxOkD4BXyw5NW77F99dJ7hZdS+hvwGlCTq1sBIuJ8\nYEdgi5TSy93eKtz8LGYs/6TW5yal9F5K6bmU0qMppZPIi7+OpIDzAosdT2/71vLcNAP/AnRExLsR\n8S4wGjgyIt4h/62zSPOz2PGUkreF1Pj8LCSl1An8hVxrIf/d6dJjLL29P+B5qYmmoDSAV4Ctu7ZF\nxCrkFeQPLOrniqTUhX8QWOwvqGop/RLdFdgypTSl+3tFm5/FjWUR+9f03PRiKfLTRgs1L4uxFLB8\nb2/U+Nz8EdiAHLePKP15BLgaGJFSeo5izU9f4+ntqrFanp+FRMTK5F+S04r+7063sfT6//ug5mUI\nV08OI3/JNiSvBj+q9HrN0vvHkVch70z+Yt4APAMsN1Q1lms8pfdOJ3/BPk7+4j1CflLkstWuvZex\nXAi8Qb6cb/Vuf1botk8h5qevsRRwbn5UGsvHgfXJ5wrfA7Yq0rwsyXiKNjeLGF/P1fqFmp/Fjado\n8wOcAWxeqvULwO3k9OaDRZubxY2l3PMylIMaTf7lOa/Hn1922+cU8mUic4BbgbWrPRkDGQ95AdV4\ncif6NvAccBHwL9WuexFj6W0c84C9e+xX8/PT11gKODeXlWqcW6r5NkoNQZHmZUnGU7S5WcT47qRb\nU1C0+VnceIo2P0Ab8FLpuzYF+DXwiSLOzeLGUu558YFIkiQJqJE1BZIkqfpsCiRJEmBTIEmSSmwK\nJEkSYFMgSZJKbAokSRJgUyBJkkpsCiRJEmBTIEmSSmwKJEkSYFMgSZJK/j+hizYYMD1YVgAAAABJ\nRU5ErkJggg==\n",
+ "text/plain": [
+ "<matplotlib.figure.Figure at 0x7feca8aee0b8>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "k=5\n",
+ "N=arange(k)*10+11\n",
+ "x=zeros(k)\n",
+ "\n",
+ "for iii in range(k):\n",
+ " A=poisson2d_matrix(N[iii])\n",
+ " B,p=lu_pivot(A.copy())\n",
+ " x[iii]=(1.0*count_nonzero(B))/count_nonzero(A)\n",
+ " \n",
+ "plt.plot(N,x)\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "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