prb4_p.ipynb (8273B)
1 { 2 "cells": [ 3 { 4 "cell_type": "code", 5 "execution_count": 1, 6 "id": "84bcf23b", 7 "metadata": {}, 8 "outputs": [], 9 "source": [ 10 "import numpy as np\n", 11 "from numpy import testing as testing\n", 12 "from matplotlib import pyplot as plt" 13 ] 14 }, 15 { 16 "cell_type": "markdown", 17 "id": "f8656109", 18 "metadata": {}, 19 "source": [ 20 "# Sheet 4, Exercise 1" 21 ] 22 }, 23 { 24 "cell_type": "code", 25 "execution_count": 2, 26 "id": "0de86192", 27 "metadata": {}, 28 "outputs": [], 29 "source": [ 30 "def conjugate_gradient(A:np.ndarray, b:np.ndarray, steps: int) -> np.ndarray:\n", 31 " x_k = np.zeros(b.shape) \n", 32 " r_k = b - (A @ x_k)\n", 33 " p_k = r_k\n", 34 " for k in range(steps):\n", 35 " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n", 36 " x_k = x_k + (α_k * p_k) # x_k+1\n", 37 " r_k = r_k - α_k * (A @ p_k) # r_k+1\n", 38 " if not np.any(r_k): # stop if r_k+1 = 0 \n", 39 " break\n", 40 " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n", 41 " p_k = r_k - (β_k * p_k)\n", 42 " return x_k\n", 43 "\n", 44 "def poisson_mat(n:int, m : int =None) -> np.ndarray:\n", 45 " return 2 * np.eye(n, m) + (-1) * np.eye(n, m, k=1) + (-1) * np.eye(n, m, k=-1)" 46 ] 47 }, 48 { 49 "cell_type": "code", 50 "execution_count": 3, 51 "id": "60a4d70d", 52 "metadata": {}, 53 "outputs": [], 54 "source": [ 55 "for n in range(3, 21):\n", 56 " Q = poisson_mat(n)\n", 57 " b = np.ones(n)\n", 58 " x = np.linalg.inv(Q) @ b\n", 59 " x_k = conjugate_gradient(Q, b, 10)\n", 60 " testing.assert_allclose(x_k, x, rtol=1e-5, err_msg=f'CG failed at dim - {n}')" 61 ] 62 }, 63 { 64 "cell_type": "markdown", 65 "id": "67accdb6", 66 "metadata": {}, 67 "source": [ 68 "What happends to a matrix that is not positive definite? Consider the system\n", 69 "\\begin{align}\n", 70 " A = \n", 71 " \\begin{pmatrix}\n", 72 " 1 & 0 & 0\\\\\n", 73 " 0 & 1 & 0\\\\\n", 74 " 0 & 0 & -2\\\\\n", 75 " \\end{pmatrix}, \\qquad\n", 76 " b = \n", 77 " \\begin{pmatrix}\n", 78 " 1 \\\\\n", 79 " 1 \\\\\n", 80 " 1\n", 81 " \\end{pmatrix}.\n", 82 "\\end{align}\n", 83 "Because $A$ is indefinite, the coefficients $\\alpha_k$ are undefined, \n", 84 "i.e. there is an $r^k \\neq 0$ such that $(r^k)^TAr^k = 0$." 85 ] 86 }, 87 { 88 "cell_type": "code", 89 "execution_count": 4, 90 "id": "12b54df7", 91 "metadata": {}, 92 "outputs": [ 93 { 94 "name": "stderr", 95 "output_type": "stream", 96 "text": [ 97 "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: divide by zero encountered in double_scalars\n", 98 " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n", 99 "/tmp/ipykernel_219747/1775264609.py:12: RuntimeWarning: invalid value encountered in subtract\n", 100 " p_k = r_k - (β_k * p_k)\n", 101 "/tmp/ipykernel_219747/1775264609.py:6: RuntimeWarning: invalid value encountered in matmul\n", 102 " α_k = (p_k @ r_k) / (p_k @ A @ p_k)\n", 103 "/tmp/ipykernel_219747/1775264609.py:8: RuntimeWarning: invalid value encountered in matmul\n", 104 " r_k = r_k - α_k * (A @ p_k) # r_k+1\n", 105 "/tmp/ipykernel_219747/1775264609.py:11: RuntimeWarning: invalid value encountered in matmul\n", 106 " β_k = ((A @ p_k) @ r_k) / ((A @ p_k) @ p_k)\n" 107 ] 108 }, 109 { 110 "ename": "AssertionError", 111 "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.])", 112 "output_type": "error", 113 "traceback": [ 114 "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", 115 "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", 116 "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", 117 " \u001b[0;31m[... skipping hidden 2 frame]\u001b[0m\n", 118 "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", 119 "\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.])" 120 ] 121 } 122 ], 123 "source": [ 124 "A = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, -2]])\n", 125 "b = np.ones(3)\n", 126 "x_k = conjugate_gradient(A, b , 10)\n", 127 "testing.assert_allclose(x_k, b, rtol=1e-5 ,err_msg=f'CG failed for non-definite matrix')" 128 ] 129 }, 130 { 131 "cell_type": "markdown", 132 "id": "5c5e809e", 133 "metadata": {}, 134 "source": [ 135 "# Sheet 4, Exercise 4 1) check" 136 ] 137 }, 138 { 139 "cell_type": "code", 140 "execution_count": 5, 141 "id": "20f80867", 142 "metadata": {}, 143 "outputs": [ 144 { 145 "name": "stdout", 146 "output_type": "stream", 147 "text": [ 148 "[3. 2. 1.]\n" 149 ] 150 } 151 ], 152 "source": [ 153 "A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])\n", 154 "b = np.array([4, 0, 0])\n", 155 "x = conjugate_gradient(A, b, 10)\n", 156 "print(x)" 157 ] 158 } 159 ], 160 "metadata": { 161 "kernelspec": { 162 "display_name": "Python 3", 163 "language": "python", 164 "name": "python3" 165 }, 166 "language_info": { 167 "codemirror_mode": { 168 "name": "ipython", 169 "version": 3 170 }, 171 "file_extension": ".py", 172 "mimetype": "text/x-python", 173 "name": "python", 174 "nbconvert_exporter": "python", 175 "pygments_lexer": "ipython3", 176 "version": "3.10.3" 177 } 178 }, 179 "nbformat": 4, 180 "nbformat_minor": 5 181 }