prb2_p.ipynb (18835B)
1 { 2 "cells": [ 3 { 4 "cell_type": "code", 5 "execution_count": 2, 6 "id": "517aef32", 7 "metadata": {}, 8 "outputs": [], 9 "source": [ 10 "import numpy as np\n", 11 "import matplotlib.pyplot as plt" 12 ] 13 }, 14 { 15 "cell_type": "markdown", 16 "id": "85da4812", 17 "metadata": {}, 18 "source": [ 19 "# Sheet 2, Exercise 3" 20 ] 21 }, 22 { 23 "cell_type": "code", 24 <<<<<<< HEAD 25 "execution_count": 2, 26 ======= 27 "execution_count": 3, 28 >>>>>>> origin/master 29 "id": "e14c738f", 30 "metadata": {}, 31 "outputs": [], 32 "source": [ 33 "def gauss_siedel(A, b, k):\n", 34 " n, m = Q.shape\n", 35 " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n", 36 " L = np.reshape([Q[i][j] if i>j else 0 for i in range(n) for j in range(n)], (n ,n))\n", 37 " U = np.reshape([Q[i][j] if i<j else 0 for i in range(n) for j in range(n)], (n ,n))\n", 38 "\n", 39 " x = np.random.rand(n)\n", 40 " for i in range(k):\n", 41 " x = np.linalg.inv(D)@(b - (L + U)@x)\n", 42 " return x\n", 43 "\n", 44 "def poisson_mat(n, m=None):\n", 45 " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)\n", 46 "\n", 47 "# test\n", 48 "for n in range(5, 20):\n", 49 " Q = poisson_mat(n)\n", 50 " b = np.ones(n)\n", 51 " x = gauss_siedel(Q, b, k=1000)\n", 52 " np.testing.assert_allclose(Q@x, b, rtol=1e-5, err_msg=f'GS failed at dim - {n}')" 53 ] 54 }, 55 { 56 "cell_type": "markdown", 57 "id": "51428a2b", 58 "metadata": {}, 59 "source": [ 60 "# Sheet 2, Exercise 5" 61 ] 62 }, 63 { 64 "cell_type": "code", 65 <<<<<<< HEAD 66 "execution_count": 6, 67 ======= 68 "execution_count": 4, 69 >>>>>>> origin/master 70 "id": "caef181e", 71 "metadata": {}, 72 "outputs": [ 73 { 74 "name": "stdout", 75 "output_type": "stream", 76 "text": [ 77 <<<<<<< HEAD 78 "1\t44.76606865271526\n", 79 "2\t59.35975010638207\n", 80 "3\t22.760834328149066\n", 81 "4\t35.62184004487846\n", 82 "5\t15.344146462132612\n", 83 "6\t25.450588757787237\n", 84 "7\t11.636387156050118\n", 85 "8\t19.801556558635152\n", 86 "9\t9.412478177542988\n", 87 "10\t16.208077941288785\n", 88 "11\t7.930495393514105\n", 89 "12\t13.721435287552058\n", 90 "13\t6.872470144914724\n", 91 "14\t11.898893608722659\n", 92 "15\t6.079418049762978\n", 93 "16\t10.506063704904786\n", 94 "17\t5.463014409761165\n", 95 "18\t9.407246486889383\n", 96 "19\t4.970264368579988\n", 97 "20\t8.518437685556936\n", 98 "21\t4.567443892886474\n", 99 "22\t7.784851843354031\n", 100 "23\t4.2320702628688585\n", 101 "24\t7.1692347898822\n", 102 "25\t3.948578490221603\n", 103 "26\t6.645370572800227\n", 104 "27\t3.705850699860457\n", 105 "28\t6.194275175956779\n", 106 "29\t3.495733758960119\n", 107 "Max. and Min. Singular value are far apart from each other for uneaven p\n" 108 ======= 109 "1\t44.76606865271519\n", 110 "2\t59.35975010638131\n", 111 "3\t22.760834328149002\n", 112 "4\t35.62184004487854\n", 113 "5\t15.344146462132624\n", 114 "6\t25.450588757787248\n", 115 "7\t11.636387156050118\n", 116 "8\t19.801556558635184\n", 117 "9\t9.41247817754299\n" 118 >>>>>>> origin/master 119 ] 120 }, 121 { 122 "data": { 123 <<<<<<< HEAD 124 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAD4CAYAAABWiRm9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASiElEQVR4nO3df2xdd3nH8fczNxVWYXO7elGSsqUblRGC0SCrAhWhrqWEAVo9hCrQNgWpUvYHm4qYsib8M5i2EZaNH9IkpoyyBQloqxKSCiRC1R9ik6aCgwuBhqylakWdtDGjFnSyIA3P/vAxczI7vjf2Offrc94vKfK933vt+xwd5X7uec73fG9kJpIklehXBl2AJEnLMaQkScUypCRJxTKkJEnFMqQkScW6pMkXu/LKK3Pr1q1NvqQkqWBHjx79UWaOLvd4oyG1detWJicnm3xJSVLBIuLpCz1uu0+SVCxDSpJULENKklSsnkIqIkYi4t6I+H5EHI+IN0TEFRFxf0Q8Xv28vO5iJUnd0uuR1CeBr2bmK4HXAseB3cADmXkN8EB1X5KkNbPi7L6I+DXgTcB7ATLz58DPI+IW4IbqaQeAh4E76ijyYh2ammbfkROcnJ1j88gwu7aPMbFty6DLkiT1qJcjqauBGeBfI2IqIj4dEZcBGzPzVPWcZ4GNS/1yROyMiMmImJyZmVmbqntwaGqaPQePMT07RwLTs3PsOXiMQ1PTjdUgSVqdXkLqEuB1wKcycxvwP5zX2sv57/tY8js/MnN/Zo5n5vjo6LLXa625fUdOMHfm7Dljc2fOsu/IicZqkCStTi8h9QzwTGY+Ut2/l/nQei4iNgFUP0/XU+LFOTk719e4JKk8K4ZUZj4L/DAixqqhm4DHgPuAHdXYDuBwLRVepM0jw32NS5LK0+vsvj8HPhcR3wGuBf4O2AvcHBGPA2+u7hdj1/YxhjcMnTM2vGGIXdvHlvkNSVJpelq7LzMfBcaXeOimNa1mDS3M4nN2nyStX40uMNu0iW1bDCVJWsdcFkmSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVCxDSpJULENKklQsQ0qSVKxLenlSRDwF/BQ4C7yYmeMRcQVwN7AVeAq4NTOfr6dMSVIX9XMk9XuZeW1mjlf3dwMPZOY1wAPVfUmS1sxq2n23AAeq2weAiVVXI0nSIr2GVAJfi4ijEbGzGtuYmaeq288CG5f6xYjYGRGTETE5MzOzynIlSV3S0zkp4I2ZOR0RvwHcHxHfX/xgZmZE5FK/mJn7gf0A4+PjSz5HkqSl9HQklZnT1c/TwJeA64DnImITQPXzdF1FSpK6acWQiojLIuJlC7eBtwDfBe4DdlRP2wEcrqtISVI39dLu2wh8KSIWnv/5zPxqRHwTuCcibgOeBm6tr0xJUhetGFKZ+STw2iXG/xu4qY6iJEkCV5yQJBXMkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBWr16/qKMKhqWn2HTnBydk5No8Ms2v7GBPbtgy6LElSTdZNSB2ammbPwWPMnTkLwPTsHHsOHgMwqCSppdZNu2/fkRO/DKgFc2fOsu/IiQFVJEmq27oJqZOzc32NS5LWv3UTUptHhvsalyStf+smpHZtH2N4w9A5Y8Mbhti1fWxAFUmS6rZuJk4sTI5wdp8kdce6CSmYDypDSZK6Y920+yRJ3WNISZKKZUhJkoplSEmSimVISZKKZUhJkoplSEmSimVISZKKZUhJkorVc0hFxFBETEXEl6v7V0fEIxHxRETcHRGX1lemJKmL+jmSuh04vuj+R4GPZ+YrgOeB29ayMEmSegqpiLgKeDvw6ep+ADcC91ZPOQBM1FCfJKnDej2S+gTwl8Avqvu/Dsxm5ovV/WeAJVd+jYidETEZEZMzMzOrqVWS1DErhlREvAM4nZlHL+YFMnN/Zo5n5vjo6OjF/AlJUkf18lUd1wN/EBFvA14C/CrwSWAkIi6pjqauAqbrK1OS1EUrHkll5p7MvCoztwLvBh7MzD8CHgLeVT1tB3C4tiolSZ20muuk7gA+EBFPMH+O6s61KUmSpHl9fTNvZj4MPFzdfhK4bu1LkiRpnitOSJKKZUhJkoplSEmSimVISZKKZUhJkoplSEmSitXXFHQt7dDUNPuOnODk7BybR4bZtX2MiW1LLmUoSeqDIbVKh6am2XPwGHNnzgIwPTvHnoPHAAwqSVol232rtO/IiV8G1IK5M2fZd+TEgCqSpPYwpFbp5OxcX+OSpN4ZUqu0eWS4r3FJUu8MqVXatX2M4Q1D54wNbxhi1/axAVUkSe3hxIlVWpgc4ew+SVp7htQamNi2xVCSpBrY7pMkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBXLkJIkFcuQkiQVy5CSJBVrxZCKiJdExDci4tsR8b2I+HA1fnVEPBIRT0TE3RFxaf3lSpK6pJcjqZ8BN2bma4FrgbdGxOuBjwIfz8xXAM8Dt9VWpSSpk1YMqZz3QnV3Q/UvgRuBe6vxA8BEHQVKkrqrp3NSETEUEY8Cp4H7gR8As5n5YvWUZ4Alv1ApInZGxGRETM7MzKxByZKkrugppDLzbGZeC1wFXAe8stcXyMz9mTmemeOjo6MXV6UkqZP6mt2XmbPAQ8AbgJGIWPhm36uA6bUtTZLUdb3M7huNiJHq9jBwM3Cc+bB6V/W0HcDhmmqUJHXUJSs/hU3AgYgYYj7U7snML0fEY8BdEfE3wBRwZ411SpI6aMWQyszvANuWGH+S+fNTkiTVwhUnJEnFMqQkScUypCRJxTKkJEnFMqQkScUypCRJxTKkJEnF6uViXhXg0NQ0+46c4OTsHJtHhtm1fYyJbUuu6StJrWFIrQOHpqbZc/AYc2fOAjA9O8eeg8cADCpJrWa7bx3Yd+TELwNqwdyZs+w7cmJAFUlSMwypdeDk7Fxf45LUFobUOrB5ZLivcUlqC0NqHdi1fYzhDUPnjA1vGGLX9rEBVSRJzXDixDqwMDnC2X2SusaQWicmtm0xlCR1ju0+SVKxDClJUrFs953HlR0kqRyG1CKu7CBJZbHdt4grO0hSWQypRVzZQZLKYkgt4soOklQWQ2oRV3aQpLI4cWIRV3aQpLIYUudxZQdJKoftPklSsVYMqYh4eUQ8FBGPRcT3IuL2avyKiLg/Ih6vfl5ef7mSpC7p5UjqReAvMvNVwOuB90XEq4DdwAOZeQ3wQHVf69ihqWmu3/sgV+/+CtfvfZBDU9ODLklSx60YUpl5KjO/Vd3+KXAc2ALcAhyonnYAmKipRjVgYbWN6dk5kv9bbcOgkjRIfZ2TioitwDbgEWBjZp6qHnoW2Li2palJrrYhqUQ9h1REvBT4IvD+zPzJ4scyM4Fc5vd2RsRkREzOzMysqljVx9U2JJWop5CKiA3MB9TnMvNgNfxcRGyqHt8EnF7qdzNzf2aOZ+b46OjoWtSsGrjahqQS9TK7L4A7geOZ+bFFD90H7Khu7wAOr315aoqrbUgqUS8X814P/AlwLCIercY+COwF7omI24CngVtrqVCNcLUNSSWK+dNJzRgfH8/JycnGXk+SVLaIOJqZ48s97ooTkqRiGVKSpGIZUpKkYhlSkqRiGVKSpGL5fVJq3KGpaae6S+qJIaVGLSxku7BO4MJCtoBBJen/sd2nRrmQraR+eCQ1AF1ud7mQraR+eCTVsK5/b5ML2UrqhyHVsK63u1zIVlI/bPc1rOvtLheyldQPQ6phm0eGmV4ikLrU7prYtsVQktQT230Ns90lSb3zSKphtrua0eUZlFKbGFIDYLurXl4wLLWH7T61TtdnUEptYkipdbo+g1JqE0NKreMFw1J7GFJqHWdQSu3hxAm1TlMzKJ1BKNXPkGqxLr+J1j2D0hmEUjNs97VU1xeyrZszCKVmGFIt5ZtovZxBKDXDkGop30Tr5QxCqRmGVEv5JlqvJmcQHpqa5vq9D3L17q9w/d4HbdmqUwyplnIadr0mtm3hI+98DVtGhglgy8gwH3nna2qZQei5RXXZirP7IuIzwDuA05n56mrsCuBuYCvwFHBrZj5fX5nqlwvZ1q+JNRgvdG7Rfaku6GUK+r8B/wR8dtHYbuCBzNwbEbur+3esfXlaDReyXf88t6iuWzGkMvPrEbH1vOFbgBuq2weAhzGkOqnL12I1oakvyXQ/qlQXe05qY2aeqm4/C2xc7okRsTMiJiNicmZm5iJfTiXyfEn9mji36H5UyVY9cSIzE8gLPL4/M8czc3x0dHS1L6eCeC1W/ZqYoOF+VMkudlmk5yJiU2aeiohNwOm1LErrg+dLmlH3ucUm96NtRfXrYo+k7gN2VLd3AIfXphytJ16L1Q5N7UfbiroYK4ZURHwB+E9gLCKeiYjbgL3AzRHxOPDm6r46pqlrsbyYtV5N7UfbiroYvczue88yD920xrVonWniWixXG69fU9fUNdVWtKXYLn5Vh1al7vMlXszajCauqWtiOr0fatrHZZFUNCdntEcTbcWmWoq2oJvjkZSK1tTFrGCbqG5NtBWb+FDj0VqzDCkVbdf2sXPeEKC+yRm+8dSv7rZiEx9qmmxB+8HJdp8K19Rq4848a4cmWopNTgBpYsp+6a1Lj6RUvCZO6jvzrB2aaCk21YJu4ohtPXQQDCkJZ561Sd0fappqQTfxwWk9zJ613SfhzDP1rqkWdBMrgayH2bMeSUk480z9aaIF3cQRW5OzZy+WISVVnHnWH8+v1auJD05NtS5Xw5CSGtLEG0LTM8/qPmLrehDW/cGpqSWxVsOQkhrizLP+NNm67HIYNtG6XA1DSmqQM89611Tr0qPCshlSUos01b5p4oitqdZlm44K2xiEhpTUMs48609bjgrb2h71OilJfWviWqGmvoyxLdcjNXkdXpPfsOyRlKSL0paZZ205KmxTe3QxQ0pSsZpoXbbleqQ2tUcXM6QkdV4bjgqbmtnZ9CoVhpQkNaANQQjNr1JhSElSS7SlPbqYISVJ6kuTq1Q4BV2SVCxDSpJULENKklQsQ0qSVCxDSpJUrMjM5l4sYgZ4+rzhK4EfNVZEWbq87dDt7e/ytkO3t99tP9dvZebocr/QaEgtWUDEZGaOD7SIAenytkO3t7/L2w7d3n63vb9tt90nSSqWISVJKlYJIbV/0AUMUJe3Hbq9/V3eduj29rvtfRj4OSlJkpZTwpGUJElLMqQkScUaWEhFxFsj4kREPBERuwdVx6BExFMRcSwiHo2IyUHXU7eI+ExEnI6I7y4auyIi7o+Ix6uflw+yxross+0fiojpav8/GhFvG2SNdYmIl0fEQxHxWER8LyJur8Zbv+8vsO1d2fcviYhvRMS3q+3/cDV+dUQ8Ur333x0Rl17w7wzinFREDAH/BdwMPAN8E3hPZj7WeDEDEhFPAeOZ2YmL+iLiTcALwGcz89XV2N8DP87MvdUHlcsz845B1lmHZbb9Q8ALmfkPg6ytbhGxCdiUmd+KiJcBR4EJ4L20fN9fYNtvpRv7PoDLMvOFiNgA/AdwO/AB4GBm3hUR/wx8OzM/tdzfGdSR1HXAE5n5ZGb+HLgLuGVAtagBmfl14MfnDd8CHKhuH2D+P3DrLLPtnZCZpzLzW9XtnwLHgS10YN9fYNs7Iee9UN3dUP1L4Ebg3mp8xX0/qJDaAvxw0f1n6NDOqyTwtYg4GhE7B13MgGzMzFPV7WeBjYMsZgD+LCK+U7UDW9fuOl9EbAW2AY/QsX1/3rZDR/Z9RAxFxKPAaeB+4AfAbGa+WD1lxfd+J04Mzhsz83XA7wPvq1pCnZXzfecuXQ/xKeB3gGuBU8A/DrSamkXES4EvAu/PzJ8sfqzt+36Jbe/Mvs/Ms5l5LXAV8x20V/b7NwYVUtPAyxfdv6oa64zMnK5+nga+xPwO7Jrnqr79Qv/+9IDraUxmPlf9B/4F8C+0eP9X5yO+CHwuMw9Ww53Y90tte5f2/YLMnAUeAt4AjETEJdVDK773DyqkvglcU83yuBR4N3DfgGppXERcVp1IJSIuA94CfPfCv9VK9wE7qts7gMMDrKVRC2/QlT+kpfu/Onl+J3A8Mz+26KHW7/vltr1D+340Ikaq28PMT5Q7znxYvat62or7fmArTlTTLj8BDAGfycy/HUghAxARv8380RPAJcDn2779EfEF4Abml+p/Dvgr4BBwD/CbzH+Fy62Z2boJBsts+w3Mt3sSeAr400XnaFojIt4I/DtwDPhFNfxB5s/NtHrfX2Db30M39v3vMj8xYoj5A6J7MvOvq/e/u4ArgCngjzPzZ8v+HZdFkiSVyokTkqRiGVKSpGIZUpKkYhlSkqRiGVKSpGIZUpKkYhlSkqRi/S9bhJKASc+kxQAAAABJRU5ErkJggg==\n", 125 ======= 126 "text/plain": [ 127 "<matplotlib.collections.PathCollection at 0x7fde841e3a00>" 128 ] 129 }, 130 "execution_count": 4, 131 "metadata": {}, 132 "output_type": "execute_result" 133 }, 134 { 135 "data": { 136 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD4CAYAAAC5S3KDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUAUlEQVR4nO3db4xcV3nH8e9TxxEroN2ELJa9JnUQkVHUNHG1ioISoeAQnAIiVoQiKJVcKZLf0ApUarB5U1GVJtQSkBdVJYvQGok/iYLjRBRhIicRNKoCazbggHET0kRk7cRLyQqoVuCYpy/mrrN29s/s3J2dM3O/H8mauWfu3Tk+szO/PeeeeyYyE0mSSvMHva6AJEnzMaAkSUUyoCRJRTKgJElFMqAkSUW6YDWf7JJLLslNmzat5lNKkgp25MiRX2TmyHyPrWpAbdq0ifHx8dV8SklSwSLiuYUec4hPklQkA0qSVCQDSpJUpLYCKiKGI+K+iPhpRByLiLdFxMUR8VBEPFXdXtTtykqSmqPdHtRdwLcy863AVcAxYDdwODMvBw5X25IkrYglZ/FFxB8Bbwf+CiAzfwf8LiJuAW6odtsPPAp8ohuVbMfBiUn2HjrOiekZNgwPsWvbZrZvGe1VdSRJNbXTg7oMmAL+LSImIuILEfFaYF1mnqz2eQFYN9/BEbEzIsYjYnxqamplan2egxOT7DlwlMnpGRKYnJ5hz4GjHJyY7MrzSZK6r52AugD4M+BfM3ML8H+cN5yXre/smPd7OzJzX2aOZebYyMi812LVtvfQcWZOnzmnbOb0GfYeOt6V55MkdV87AfU88HxmPl5t30crsF6MiPUA1e2p7lRxaSemZ5ZVLkkq35IBlZkvAD+PiM1V0Y3AT4AHgR1V2Q7gga7UsA0bhoeWVS5JKl+7s/j+BvhyRPwIuBr4J+BO4KaIeAp4Z7XdE7u2bWZo7ZpzyobWrmHXts0LHCFJKl1ba/Fl5hPA2DwP3biitenQ7Gw9Z/FJ0uBY1cViu2n7llEDSZIGiEsdSZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKdEE7O0XEs8CvgTPAy5k5FhEXA/cAm4Bngdsy86XuVFOS1DTL6UG9IzOvzsyxans3cDgzLwcOV9uSJK2IOkN8twD7q/v7ge21ayNJUqXdgErg2xFxJCJ2VmXrMvNkdf8FYN18B0bEzogYj4jxqampmtWVJDVFW+eggOszczIi3gg8FBE/nftgZmZE5HwHZuY+YB/A2NjYvPtIknS+tnpQmTlZ3Z4C7geuAV6MiPUA1e2pblVSktQ8SwZURLw2Il4/ex94F/Ak8CCwo9ptB/BAtyopSWqedob41gH3R8Ts/l/JzG9FxPeBeyPiduA54LbuVVOS1DRLBlRmPgNcNU/5/wI3dqNSkiS5koQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlI7a5mXoSDE5PsPXScE9MzbBgeYte2zWzfMtrrakmSuqBvAurgxCR7Dhxl5vQZACanZ9hz4CiAISVJA6hvhvj2Hjp+NpxmzZw+w95Dx3tUI0lSN/VNQJ2YnllWuSSpv/VNQG0YHlpWuSSpv/VNQO3atpmhtWvOKRtau4Zd2zb3qEaSpG7qm0kSsxMhnMUnSc3QNwEFrZAykCSpGfpmiE+S1CwGlCSpSAaUJKlIBpQkqUgGlCSpSAaUJKlIBpQkqUgGlCSpSG0HVESsiYiJiPhGtX1ZRDweEU9HxD0RcWH3qilJaprl9KA+Ahybs/0Z4HOZ+RbgJeD2layYJKnZ2gqoiNgIvAf4QrUdwFbgvmqX/cD2LtRPktRQ7fagPg98HPh9tf0GYDozX662nwfmXSQvInZGxHhEjE9NTdWpqySpQZYMqIh4L3AqM4908gSZuS8zxzJzbGRkpJMfIUlqoHZWM78OeF9EvBt4DfCHwF3AcERcUPWiNgKT3aumJKlpluxBZeaezNyYmZuADwAPZ+aHgEeA91e77QAe6FotJUmNU+c6qE8AfxsRT9M6J3X3ylRJkqRlfmFhZj4KPFrdfwa4ZuWrJEmSK0lIkgplQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkoq0rOugBtXBiUn2HjrOiekZNgwPsWvbZrZvmXftW0nSKml8QB2cmGTPgaPMnD4DwOT0DHsOHAUwpCSphxo/xLf30PGz4TRr5vQZ9h463qMaSZLAgOLE9MyyyiVJq6PxAbVheGhZ5ZKk1dH4gNq1bTNDa9ecUza0dg27tm3uUY0kSeAkibMTIZzFJ0llaXxAQSukDCRJKkvjh/gkSWUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRTKgJElFMqAkSUUyoCRJRVoyoCLiNRHxvYj4YUT8OCI+VZVfFhGPR8TTEXFPRFzY/epKkpqinR7Ub4GtmXkVcDVwc0RcC3wG+FxmvgV4Cbi9a7WUJDXOkgGVLb+pNtdW/xLYCtxXle8HtnejgpKkZmrrHFRErImIJ4BTwEPAz4DpzHy52uV5YN7vq4iInRExHhHjU1NTK1BlSVITtBVQmXkmM68GNgLXAG9t9wkyc19mjmXm2MjISGe1lCQ1zrJm8WXmNPAI8DZgOCJmv/BwIzC5slWTJDVZO7P4RiJiuLo/BNwEHKMVVO+vdtsBPNClOkqSGqidr3xfD+yPiDW0Au3ezPxGRPwE+FpE/CMwAdzdxXpKkhpmyYDKzB8BW+Ypf4bW+ShJklacK0lIkopkQEmSimRASZKKZEBJkopkQEmSimRASZKKZEBJkorUzoW6WsTBiUn2HjrOiekZNgwPsWvbZrZvmXfdXEnSMhhQNRycmGTPgaPMnD4DwOT0DHsOHAUwpCSpJof4ath76PjZcJo1c/oMew8d71GNJGlwGFA1nJieWVa5JKl9BlQNG4aHllUuSWqfAVXDrm2bGVq75pyyobVr2LVtc49qJEmDw0kSNcxOhHAWnyStPAOqpu1bRg0kSeoCh/gkSUVqRA/Ki2klqf8MfEB5Ma0k9aeBH+LzYlpJ6k8DH1BeTCtJ/WngA8qLaSWpPw18QHkxrST1p4GfJOHFtJLUnwY+oMCLaSWpHw38EJ8kqT8ZUJKkIi05xBcRbwK+BKwDEtiXmXdFxMXAPcAm4Fngtsx8qXtVHSyubiFJi2unB/Uy8LHMvAK4FvhwRFwB7AYOZ+blwOFqW22YXd1icnqG5JXVLQ5OTPa6apJUjCUDKjNPZuYPqvu/Bo4Bo8AtwP5qt/3A9i7VceC4uoUkLW1Z56AiYhOwBXgcWJeZJ6uHXqA1BDjfMTsjYjwixqempurUdWC4uoUkLa3tgIqI1wFfBz6amb+a+1hmJq3zU6+Smfsycywzx0ZGRmpVdlC4uoUkLa2tgIqItbTC6cuZeaAqfjEi1lePrwdOdaeKg8fVLSRpaUsGVEQEcDdwLDM/O+ehB4Ed1f0dwAMrX73BtH3LKHfceiWjw0MEMDo8xB23XuksPkmaI1qjc4vsEHE98F3gKPD7qviTtM5D3QtcCjxHa5r5Lxf7WWNjYzk+Pl63zpKkARERRzJzbL7HlrwOKjP/E4gFHr6xTsUkSVqIK0lIkopkQEmSimRASZKKZEBJkorUiO+DGjQuNCupCQyoPjO70OzsWn6zC80ChpSkgeIQX59xoVlJTWEPagmlDae50KykprAHtYgSv7fJhWYlNYUBtYgSh9NcaFZSUzjEt4gSh9NmhxdLGnaUpG4woBaxYXiIyXnCqNfDadu3jBpIkgaeQ3yLGMThtIMTk1x358Nctvs/uO7Oh3t6Pk2SFmMPahGDNpzmNVSS+okBtYRBGk5bbNLHoPwfJQ0Oh/gapMRJH5K0EAOqQbyGSlI/MaAaZBAnfUgaXJ6DapBBm/QhabAZUA1TZ9JHaesSShpsBlQXDdIHulPUJa02z0F1SYkLzdZR4rqEkgabAdUlg/aB7hR1SavNgOqSQftAd4q6pNVmQHXJoH2g152i7hqAkpZryYCKiC9GxKmIeHJO2cUR8VBEPFXdXtTdavafQbvmaPuWUe649UpGh4cIYHR4iDtuvbKtCRKDdj5O0uqIzFx8h4i3A78BvpSZf1KV/TPwy8y8MyJ2Axdl5ieWerKxsbEcHx9fgWr3h0GaxVfHdXc+PO/XlowOD/HY7q09qJGkUkTEkcwcm++xJaeZZ+Z3ImLTecW3ADdU9/cDjwJLBlTTeM1RS53zcYPUDpKWp9ProNZl5snq/gvAuoV2jIidwE6ASy+9tMOna5ZBu+ao0y9+HLR2kLQ8tSdJZGuMcMFxwszcl5ljmTk2MjJS9+kaYdCmqHd6Pm7Q2kHS8nTag3oxItZn5smIWA+cWslKNd2gTVHvdA3Auu3g8KDU3zoNqAeBHcCd1e0DK1YjdTwkVrJOzsfVaQeHB6X+1840868C/wVsjojnI+J2WsF0U0Q8Bbyz2tYK8Zqjljrt4PCg1P/amcX3wQUeunGF66JKna/FGKSeQ512cHhQ6n+uZl6oTqeoL9Zz6McP2E7bweFBqf+51NGAGbQJFp3q5fDgoAyxSr1mD2rA1O05DMqwVq+GB+19SSvHgBowu7ZtPucDEtrrOQziB2svhgfrDrEO0h8JUl0O8Q2YThd1ddbbK+oMD65E78tFdaUWe1ADqJOeg7PeXlFneLBXva9Ban9plgElwFlv5+t0eLDTIVbo/I+Euu1vuKlUDvEJcNbbSqnzvVmdfsllnfavO6w4SK+dymMPSoCz3lbSave+6rR/3WFFe27qJgNKZznrrbc6/SOhTvsbbiqZAaXaenHeBQbzQ66TPxLqtH/Twq3E11wL8xyUauvFeRfo3bmX0s671Gn/Ouce67x23Qq3xfTyfFtpvzP9wh6UVkQ/zXqDzv+CL7XX1mn71zn32G89t1722kr8nekHBpR6qlfXHPXbh9zs8YZb5697r4Yk+3E4s5RQNKDUc73offXbh5zh9opOX/denW/rt55+STNrDSj1rSZ9yPVjuC113GqHW6+GJPutp1/SV/YYUOprTfmQ67dw63YoLva6L3RsO6/5QsfWed37radf0lf2GFBqrE7Crd/Ou0Bvwq3UCQlLBdtSz7vQ675YoNY5dqnfmYWO7VWPb6UZUNIy9dN5F+hNuPXjhISljl3odW8n2Do9drHfmcWObef3pU5vcbUmURhQ0ipqSrj144SEUs/3LPQ7s9ixj+3eenafhXptdXqLqzWJwoCS+kQ/hVs/Tkjot/M9Sx272O9Lp73Fdo5dSa4kITXA9i2jPLZ7K/9z53t4bPfWtj9IOl2lolerW/Ti2DoravTq2H6ZRGEPStKi6vTcVrvH14tjezXDrx97qcsVmbniP3QhY2NjOT4+vmrPJ0mroVerNtS5Rm2+cGunl1vn2PlExJHMHJv3MQNKkpqnlKWQDChJUpEWC6hakyQi4uaIOB4RT0fE7jo/S5KkuToOqIhYA/wL8OfAFcAHI+KKlaqYJKnZ6vSgrgGezsxnMvN3wNeAW1amWpKkpqsTUKPAz+dsP1+VnSMidkbEeESMT01N1Xg6SVKTdP1C3czcl5ljmTk2MjLS7aeTJA2IOhfqTgJvmrO9sSpb0JEjR34REc+18bMvAX5Ro25NYTu1x3Zqj+3UHtupPe220x8v9EDH08wj4gLgv4EbaQXT94G/yMwfd/QDz/3Z4wtNO9QrbKf22E7tsZ3aYzu1ZyXaqeMeVGa+HBF/DRwC1gBfXIlwkiQJaq7Fl5nfBL65QnWRJOmsUlcz39frCvQJ26k9tlN7bKf22E7tqd1Oq7rUkSRJ7Sq1ByVJajgDSpJUpOICygVo5xcRX4yIUxHx5JyyiyPioYh4qrq9qJd1LEFEvCkiHomIn0TEjyPiI1W5bTVHRLwmIr4XET+s2ulTVfllEfF49f67JyIu7HVdey0i1kTERER8o9q2jeYREc9GxNGIeCIixquyWu+7ogLKBWgX9e/AzeeV7QYOZ+blwOFqu+leBj6WmVcA1wIfrn6HbKtz/RbYmplXAVcDN0fEtcBngM9l5luAl4Dbe1fFYnwEODZn2zZa2Dsy8+o51z/Vet8VFVC4AO2CMvM7wC/PK74F2F/d3w9sX806lSgzT2bmD6r7v6b1wTKKbXWObPlNtbm2+pfAVuC+qrzx7RQRG4H3AF+otgPbaDlqve9KC6i2FqDVWesy82R1/wVgXS8rU5qI2ARsAR7HtnqVaujqCeAU8BDwM2A6M1+udvH9B58HPg78vtp+A7bRQhL4dkQciYidVVmt912tC3VVjszMiPCagUpEvA74OvDRzPxV6w/fFtuqJTPPAFdHxDBwP/DW3taoLBHxXuBUZh6JiBt6XJ1+cH1mTkbEG4GHIuKncx/s5H1XWg9q2QvQNtyLEbEeoLo91eP6FCEi1tIKpy9n5oGq2LZaQGZOA48AbwOGq3U2wfffdcD7IuJZWqcbtgJ3YRvNKzMnq9tTtP7guYaa77vSAur7wOXVLJkLgQ8AD/a4TiV7ENhR3d8BPNDDuhShOkdwN3AsMz875yHbao6IGKl6TkTEEHATrfN1jwDvr3ZrdDtl5p7M3JiZm2h9Fj2cmR/CNnqViHhtRLx+9j7wLuBJar7viltJIiLeTWvcd3YB2k/3tkZliIivAjfQWsL+ReDvgYPAvcClwHPAbZl5/kSKRomI64HvAkd55bzBJ2mdh7KtKhHxp7ROWq+h9YfqvZn5DxHxZlq9hYuBCeAvM/O3vatpGaohvr/LzPfaRq9Wtcn91eYFwFcy89MR8QZqvO+KCyhJkqC8IT5JkgADSpJUKANKklQkA0qSVCQDSpJUJANKklQkA0qSVKT/B002BwKIALUUAAAAAElFTkSuQmCC\n", 137 >>>>>>> origin/master 138 "text/plain": [ 139 "<Figure size 504x288 with 1 Axes>" 140 ] 141 }, 142 "metadata": { 143 "needs_background": "light" 144 }, 145 "output_type": "display_data" 146 } 147 ], 148 "source": [ 149 "def neumann_polynomial_preconditioner(n, p):\n", 150 " Q = poisson_mat(n)\n", 151 " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n", 152 <<<<<<< HEAD 153 " N = D - Q\n", 154 ======= 155 " N = D-Q\n", 156 >>>>>>> origin/master 157 " C_p = np.zeros([n, n])\n", 158 " for k in range(p+1):\n", 159 " C_p += np.linalg.matrix_power(N @ np.linalg.inv(D), k)\n", 160 " return np.linalg.inv(D) @ C_p\n", 161 " \n", 162 " \n", 163 "n = 20\n", 164 "Q = poisson_mat(n)\n", 165 <<<<<<< HEAD 166 "P = np.arange(1, 30)\n", 167 ======= 168 "P = np.arange(1, 50)\n", 169 >>>>>>> origin/master 170 "cond_2 = []\n", 171 "for p in P:\n", 172 " C_p = neumann_polynomial_preconditioner(n, p)\n", 173 " cond_2.append(np.linalg.cond(C_p @ Q, p=2))\n", 174 " if p in np.arange(1, 10):\n", 175 " print(p, cond_2[p-1], sep='\\t')\n", 176 " \n", 177 "plt.figure(figsize=[7, 4])\n", 178 "plt.scatter(P, cond_2)" 179 ] 180 }, 181 { 182 "cell_type": "code", 183 "execution_count": null, 184 "id": "2469e597", 185 "metadata": {}, 186 "outputs": [], 187 "source": [] 188 } 189 ], 190 "metadata": { 191 "kernelspec": { 192 "display_name": "Python 3 (ipykernel)", 193 "language": "python", 194 "name": "python3" 195 }, 196 "language_info": { 197 "codemirror_mode": { 198 "name": "ipython", 199 "version": 3 200 }, 201 "file_extension": ".py", 202 "mimetype": "text/x-python", 203 "name": "python", 204 "nbconvert_exporter": "python", 205 "pygments_lexer": "ipython3", 206 "version": "3.10.3" 207 } 208 }, 209 "nbformat": 4, 210 "nbformat_minor": 5 211 }