commit 37bfaa0982dc479d5621f54fd03623212814492b
parent fb582e38d41ae807e1907b2352f66fe2582a36c3
Author: miksa <milutin@popovic.xyz>
Date: Sun, 27 Mar 2022 21:17:41 +0200
prb4 programming and workout done need to make latex file
Diffstat:
6 files changed, 205 insertions(+), 156 deletions(-)
diff --git a/num_ana/prog/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/Untitled-checkpoint.ipynb
@@ -1,6 +0,0 @@
-{
- "cells": [],
- "metadata": {},
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/num_ana/prog/.ipynb_checkpoints/prb2_p-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/prb2_p-checkpoint.ipynb
@@ -1,142 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 1,
- "id": "517aef32",
- "metadata": {},
- "outputs": [],
- "source": [
- "import numpy as np\n",
- "import matplotlib.pyplot as plt"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "85da4812",
- "metadata": {},
- "source": [
- "# Sheet 2, Exercise 3"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "id": "e14c738f",
- "metadata": {},
- "outputs": [],
- "source": [
- "def gauss_siedel(A, b, k):\n",
- " n, m = Q.shape\n",
- " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
- " L = np.reshape([Q[i][j] if i>j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
- " U = np.reshape([Q[i][j] if i<j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
- "\n",
- " x = np.random.rand(n)\n",
- " for i in range(k):\n",
- " x = np.linalg.inv(D)@(b - (L + U)@x)\n",
- " return x\n",
- "\n",
- "def poisson_mat(n, m=None):\n",
- " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)\n",
- "\n",
- "# test\n",
- "for n in range(5, 20):\n",
- " Q = poisson_mat(n)\n",
- " b = np.ones(n)\n",
- " x = gauss_siedel(Q, b, k=1000)\n",
- " np.testing.assert_allclose(Q@x, b, rtol=1e-5, err_msg=f'GS failed at dim - {n}')"
- ]
- },
- {
- "cell_type": "markdown",
- "id": "51428a2b",
- "metadata": {},
- "source": [
- "# Sheet 2, Exercise 5"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "id": "caef181e",
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "1\t44.76606865271519\n",
- "2\t59.35975010638131\n",
- "3\t22.760834328149002\n",
- "4\t35.62184004487854\n",
- "5\t15.344146462132624\n",
- "6\t25.450588757787248\n",
- "7\t11.636387156050118\n",
- "8\t19.801556558635184\n",
- "9\t9.41247817754299\n",
- "Max. and Min. Singular value are far apart from each other for uneaven p\n"
- ]
- },
- {
- "data": {
- "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",
- "text/plain": [
- "<Figure size 504x288 with 1 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- },
- "output_type": "display_data"
- }
- ],
- "source": [
- "def neumann_polynomial_preconditioner(n, p):\n",
- " Q = poisson_mat(n)\n",
- " D = np.reshape([Q[i][j] if i==j else 0 for i in range(n) for j in range(n)], (n ,n))\n",
- " N = D-Q\n",
- " C_p = np.zeros([n, n])\n",
- " for k in range(p+1):\n",
- " C_p += np.linalg.matrix_power(N @ np.linalg.inv(D), k)\n",
- " return np.linalg.inv(D) @ C_p\n",
- " \n",
- " \n",
- "n = 20\n",
- "Q = poisson_mat(n)\n",
- "P = np.arange(1, 50)\n",
- "cond_2 = []\n",
- "for p in P:\n",
- " C_p = neumann_polynomial_preconditioner(n, p)\n",
- " cond_2.append(np.linalg.cond(C_p @ Q, p=2))\n",
- " if p in np.arange(1, 10):\n",
- " print(p, cond_2[p-1], sep='\\t')\n",
- " \n",
- "plt.figure(figsize=[7, 4])\n",
- "plt.scatter(P, cond_2)\n",
- "print(\"Max. and Min. Singular value are far apart from each other for uneaven p\")"
- ]
- }
- ],
- "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.10.2"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/num_ana/prog/prb2_p.ipynb b/num_ana/prog/prb2_p.ipynb
@@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 2,
"id": "517aef32",
"metadata": {},
"outputs": [],
@@ -21,7 +21,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 3,
"id": "e14c738f",
"metadata": {},
"outputs": [],
@@ -58,7 +58,7 @@
},
{
"cell_type": "code",
- "execution_count": 5,
+ "execution_count": 4,
"id": "caef181e",
"metadata": {},
"outputs": [
@@ -74,12 +74,21 @@
"6\t25.450588757787248\n",
"7\t11.636387156050118\n",
"8\t19.801556558635184\n",
- "9\t9.41247817754299\n",
- "Max. and Min. Singular value are far apart from each other for uneaven p\n"
+ "9\t9.41247817754299\n"
]
},
{
"data": {
+ "text/plain": [
+ "<matplotlib.collections.PathCollection at 0x7fde841e3a00>"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
"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",
"text/plain": [
"<Figure size 504x288 with 1 Axes>"
@@ -113,9 +122,16 @@
" print(p, cond_2[p-1], sep='\\t')\n",
" \n",
"plt.figure(figsize=[7, 4])\n",
- "plt.scatter(P, cond_2)\n",
- "print(\"Max. and Min. Singular value are far apart from each other for uneaven p\")"
+ "plt.scatter(P, cond_2)"
]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2469e597",
+ "metadata": {},
+ "outputs": [],
+ "source": []
}
],
"metadata": {
@@ -134,7 +150,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.10.2"
+ "version": "3.10.3"
}
},
"nbformat": 4,
diff --git a/num_ana/prog/prb4_p.ipynb b/num_ana/prog/prb4_p.ipynb
@@ -0,0 +1,181 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "84bcf23b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from numpy import testing as testing\n",
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f8656109",
+ "metadata": {},
+ "source": [
+ "# Sheet 4, Exercise 1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "0de86192",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def conjugate_gradient(A:np.ndarray, b:np.ndarray, steps: int) -> np.ndarray:\n",
+ " x_k = np.zeros(b.shape) \n",
+ " r_k = b - (A @ x_k)\n",
+ " p_k = r_k\n",
+ " for k in range(steps):\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ " x_k = x_k + (α_k * p_k) # x_k+1\n",
+ " r_k = r_k - α_k * (A @ p_k) # r_k+1\n",
+ " if not np.any(r_k): # stop if r_k+1 = 0 \n",
+ " break\n",
+ " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n",
+ " p_k = r_k - (β_k * p_k)\n",
+ " return x_k\n",
+ "\n",
+ "def poisson_mat(n:int, m : int =None) -> np.ndarray:\n",
+ " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "60a4d70d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for n in range(3, 21):\n",
+ " Q = poisson_mat(n)\n",
+ " b = np.ones(n)\n",
+ " x = np.linalg.inv(Q) @ b\n",
+ " x_k = conjugate_gradient(Q, b, 10)\n",
+ " testing.assert_allclose(x_k, x, rtol=1e-5, err_msg=f'CG failed at dim - {n}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "67accdb6",
+ "metadata": {},
+ "source": [
+ "What happends to a matrix that is not positive definite? Consider the system\n",
+ "\\begin{align}\n",
+ " A = \n",
+ " \\begin{pmatrix}\n",
+ " 1 & 0 & 0\\\\\n",
+ " 0 & 1 & 0\\\\\n",
+ " 0 & 0 & -2\\\\\n",
+ " \\end{pmatrix}, \\qquad\n",
+ " b = \n",
+ " \\begin{pmatrix}\n",
+ " 1 \\\\\n",
+ " 1 \\\\\n",
+ " 1\n",
+ " \\end{pmatrix}.\n",
+ "\\end{align}\n",
+ "Because $A$ is indefinite, the coefficients $\\alpha_k$ are undefined, \n",
+ "i.e. there is an $r^k \\neq 0$ such that $(r^k)^TAr^k = 0$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "12b54df7",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: divide by zero encountered in double_scalars\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:12: RuntimeWarning: invalid value encountered in subtract\n",
+ " p_k = r_k - (β_k * p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: invalid value encountered in matmul\n",
+ " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n",
+ "/tmp/ipykernel_219747/1775264609.py:8: RuntimeWarning: invalid value encountered in matmul\n",
+ " r_k = r_k - α_k * (A @ p_k) # r_k+1\n",
+ "/tmp/ipykernel_219747/1775264609.py:11: RuntimeWarning: invalid value encountered in matmul\n",
+ " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n"
+ ]
+ },
+ {
+ "ename": "AssertionError",
+ "evalue": "\nNot equal to tolerance rtol=1e-05, atol=0\nCG failed for non-definite matrix\nx and y nan location mismatch:\n x: array([nan, nan, nan])\n y: array([1., 1., 1.])",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
+ "Input \u001b[0;32mIn [4]\u001b[0m, in \u001b[0;36m<cell line: 4>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m b \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mones(\u001b[38;5;241m3\u001b[39m)\n\u001b[1;32m 3\u001b[0m x_k \u001b[38;5;241m=\u001b[39m conjugate_gradient(A, b , \u001b[38;5;241m10\u001b[39m)\n\u001b[0;32m----> 4\u001b[0m \u001b[43mtesting\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43massert_allclose\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx_k\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrtol\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-5\u001b[39;49m\u001b[43m \u001b[49m\u001b[43m,\u001b[49m\u001b[43merr_msg\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mCG failed for non-definite matrix\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n",
+ " \u001b[0;31m[... skipping hidden 2 frame]\u001b[0m\n",
+ "File \u001b[0;32m~/.local/lib/python3.10/site-packages/numpy/testing/_private/utils.py:745\u001b[0m, in \u001b[0;36massert_array_compare.<locals>.func_assert_same_pos\u001b[0;34m(x, y, func, hasval)\u001b[0m\n\u001b[1;32m 740\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m bool_(x_id \u001b[38;5;241m==\u001b[39m y_id)\u001b[38;5;241m.\u001b[39mall() \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m:\n\u001b[1;32m 741\u001b[0m msg \u001b[38;5;241m=\u001b[39m build_err_msg([x, y],\n\u001b[1;32m 742\u001b[0m err_msg \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mx and y \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m location mismatch:\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 743\u001b[0m \u001b[38;5;241m%\u001b[39m (hasval), verbose\u001b[38;5;241m=\u001b[39mverbose, header\u001b[38;5;241m=\u001b[39mheader,\n\u001b[1;32m 744\u001b[0m names\u001b[38;5;241m=\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mx\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124my\u001b[39m\u001b[38;5;124m'\u001b[39m), precision\u001b[38;5;241m=\u001b[39mprecision)\n\u001b[0;32m--> 745\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAssertionError\u001b[39;00m(msg)\n\u001b[1;32m 746\u001b[0m \u001b[38;5;66;03m# If there is a scalar, then here we know the array has the same\u001b[39;00m\n\u001b[1;32m 747\u001b[0m \u001b[38;5;66;03m# flag as it everywhere, so we should return the scalar flag.\u001b[39;00m\n\u001b[1;32m 748\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(x_id, \u001b[38;5;28mbool\u001b[39m) \u001b[38;5;129;01mor\u001b[39;00m x_id\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n",
+ "\u001b[0;31mAssertionError\u001b[0m: \nNot equal to tolerance rtol=1e-05, atol=0\nCG failed for non-definite matrix\nx and y nan location mismatch:\n x: array([nan, nan, nan])\n y: array([1., 1., 1.])"
+ ]
+ }
+ ],
+ "source": [
+ "A = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, -2]])\n",
+ "b = np.ones(3)\n",
+ "x_k = conjugate_gradient(A, b , 10)\n",
+ "testing.assert_allclose(x_k, b, rtol=1e-5 ,err_msg=f'CG failed for non-definite matrix')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5c5e809e",
+ "metadata": {},
+ "source": [
+ "# Sheet 4, Exercise 4 1) check"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "20f80867",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[3. 2. 1.]\n"
+ ]
+ }
+ ],
+ "source": [
+ "A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])\n",
+ "b = np.array([4, 0, 0])\n",
+ "x = conjugate_gradient(A, b, 10)\n",
+ "print(x)"
+ ]
+ }
+ ],
+ "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.10.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb4_p.pdf b/num_ana/prog/prb4_p.pdf
Binary files differ.
diff --git a/num_ana/sheets/sheet_4.pdf b/num_ana/sheets/sheet_4.pdf
Binary files differ.
