{ "cells": [ { "cell_type": "markdown", "id": "5d32decb", "metadata": {}, "source": [ "# Métodos iterativos \n", "#### https://meet.noysi.com/metodosnumericos1" ] }, { "cell_type": "markdown", "id": "9d51f4b8", "metadata": {}, "source": [ "Vamos a aplicar el método de Jacobi para resolver el sistema\n", "\\begin{cases}\n", "3 x - y + z = 3,\\\\\n", "x + 5 y - z = 2,\\\\\n", "x + y + 4 z = 1.\n", "\\end{cases}\n" ] }, { "cell_type": "markdown", "id": "da65535b", "metadata": {}, "source": [ "Definimos las matrices del método:" ] }, { "cell_type": "code", "execution_count": 20, "id": "3f720df4", "metadata": {}, "outputs": [], "source": [ "A = matrix(RDF,[[3,-1,1],[1,5,-1],[1,1,4]])\n", "M = matrix(RDF,[[3,0,0],[0,5,0],[0,0,4]])\n", "N = M - A\n", "b = vector(RDF,[3,2,1])" ] }, { "cell_type": "markdown", "id": "5fd0f394", "metadata": {}, "source": [ "En este caso está justificado calcular la inversa, por ser inmediata y no añadir mucho error." ] }, { "cell_type": "code", "execution_count": 21, "id": "0f5fe495", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n", "0.3333333333333333 & 0.0 & 0.0 \\\\\n", "0.0 & 0.2 & 0.0 \\\\\n", "0.0 & 0.0 & 0.25\n", "\\end{array}\\right)$$" ], "text/plain": [ "[0.3333333333333333 0.0 0.0]\n", "[ 0.0 0.2 0.0]\n", "[ 0.0 0.0 0.25]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Se podría calcular directamente\n", "Mi = matrix.diagonal([k^-1 for k in M.diagonal()], sparse=true)\n", "show(Mi)" ] }, { "cell_type": "markdown", "id": "12d43b88", "metadata": {}, "source": [ "Tomamos como aproximación inicial el vector nulo:" ] }, { "cell_type": "code", "execution_count": 22, "id": "01c8ca84", "metadata": {}, "outputs": [], "source": [ "x0 = vector([0,0,0])" ] }, { "cell_type": "markdown", "id": "c4dea0ed", "metadata": {}, "source": [ "Y aplicamos el método" ] }, { "cell_type": "code", "execution_count": 23, "id": "a27970e0", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.0, 0.4, 0.25)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x1 = Mi*(N*x0+b)\n", "x1" ] }, { "cell_type": "code", "execution_count": 24, "id": "693dc08e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.0, 0.4, 0.25)" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x1 = M\\(N*x0+b)\n", "x1" ] }, { "cell_type": "code", "execution_count": 5, "id": "b1274220", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.0499999999999998, 0.25, -0.09999999999999998)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x2 = Mi*(N*x1+b)\n", "x2" ] }, { "cell_type": "code", "execution_count": 6, "id": "e3495f57", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.1166666666666667, 0.17000000000000004, -0.07499999999999996)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x3 = Mi*(N*x2+b)\n", "x3" ] }, { "cell_type": "markdown", "id": "55540f05", "metadata": {}, "source": [ "Comparamos con la solución exacta." ] }, { "cell_type": "code", "execution_count": 7, "id": "a9fadeab", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "((1.078125, 0.17187500000000003, -0.062499999999999986),\n", " (1.1166666666666667, 0.17000000000000004, -0.07499999999999996),\n", " 0.0405613818113294)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A\\b, x3, (A\\b - x3).norm().n()" ] }, { "cell_type": "code", "execution_count": 66, "id": "71ecf691", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.4659658624115353, 0.9746794344808964)" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MJ = matrix(RDF,[[1,0,0],[0,1,0],[0,0,1]])\n", "MG = matrix(RDF,[[1,0,0],[0,1,0],[1,0,1]])\n", "NJ = matrix([[0,1.9,1.6],[0,0,0.5],[-1,0,0.8]])\n", "NG = matrix([[0,1.9,1.6],[0,0,0.5],[0,0,0.8]])\n", "max([k.abs() for k in (MJ.inverse()*NJ).eigenvalues()]),max([k.abs() for k in (MG.inverse()*NG).eigenvalues()])" ] }, { "cell_type": "code", "execution_count": 67, "id": "da35aab2", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n", "1.0 & -1.9 & -1.6 \\\\\n", "0.0 & 1.0 & -0.5 \\\\\n", "1.0 & 0.0 & 0.19999999999999996\n", "\\end{array}\\right) \\left(\\begin{array}{rrr}\n", "1.0 & -1.9 & -1.6 \\\\\n", "0.0 & 1.0 & -0.5 \\\\\n", "1.0 & 0.0 & 0.19999999999999996\n", "\\end{array}\\right)$$" ], "text/plain": [ "[ 1.0 -1.9 -1.6]\n", "[ 0.0 1.0 -0.5]\n", "[ 1.0 0.0 0.19999999999999996] [ 1.0 -1.9 -1.6]\n", "[ 0.0 1.0 -0.5]\n", "[ 1.0 0.0 0.19999999999999996]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "show(MJ - NJ,MG - NG)" ] }, { "cell_type": "code", "execution_count": null, "id": "1d6ec3d3", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 38, "id": "042b07c5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[4.1248854197645715, -0.7615571818318911, 0.636671762067317]" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matrix(RDF,A(x=1)).eigenvalues()" ] }, { "cell_type": "markdown", "id": "bdefdaff", "metadata": {}, "source": [ "
\n", " Ejercicio 1. \n", " a) Aplica 5 pasos del método de Jacobi para aproximar una solución del sistema y compara los valores obtenidos con la solución exacta $$\\left(\\begin{array}{rrr} 15.0 & -1.0 & 1.0 \\\\ 2.0 & -5.0 & 1.0 \\\\ 1.0 & 1.0 & -3.0 \\end{array}\\right)x = \\left(\\begin{array}{r} 1.0 \\\\ 2.0 \\\\ 3.0 \\end{array}\\right).$$\n", " \n", " b) Aplica 5 pasos del método de Gauss-Seidel para aproximar una solución del sistema anterior y compara los valores obtenidos con la solución exacta y con el método de Jacobi.\n", " \n", " c) Encuentra una matriz tal que el método de Jacobi converja y el de Gauss-Seidel no. \n", " \n", " d) Encuentra una matriz tal que el método de Gauss-Seidel converja y el de Jacobi no. \n", " \n", " e) En los apartados a) y b), aplica los métodos hasta que el error estimado con el criterio de parada sea menor que $10^{-3}$. Compara el error estimado con el error real.\n", " \n", " f) (Para valientes) Elige una matriz estrictamente diagonal dominante y un término independiente. Representa en una gráfica la estimación del error del criterio de parada para los métodos de Jacobi y Gauss-Seidel y el error exacto, para un número de pasos entre 3 y 20.\n", "
" ] }, { "cell_type": "code", "execution_count": 16, "id": "26879fa8", "metadata": {}, "outputs": [], "source": [ "x0 = vector([0,0,0])" ] }, { "cell_type": "code", "execution_count": 17, "id": "204b262a", "metadata": {}, "outputs": [], "source": [ "x1 = copy(x0)" ] }, { "cell_type": "code", "execution_count": 18, "id": "f08e3bc4", "metadata": {}, "outputs": [], "source": [ "x1[0] = 1" ] }, { "cell_type": "code", "execution_count": 19, "id": "0382edff", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(0, 0, 0)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x0" ] }, { "cell_type": "code", "execution_count": null, "id": "5d45cf67", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a9c92afe", "metadata": {}, "source": [ "# Representación de matrices\n", "Un método para visualizar matrices cuando la dimensión es muy grande es mediante su representación como imagen. \n", "\n", "Para ello, la matriz $A\\in\\mathcal{M}_n$, consideramos una cuadrícula de $nxn$ píxeles. Numeramos los píxeles comenzando desde la esquina superior izquierda (ojo - esto depende fuertemente del lenguaje y librerías que estemos utilizando). Así, la posición $(1,1)$ será la esquina superior izquierda y la $(1,n)$ la superior derecha.\n", "\n", "En cada posición dibujaremos un color según el valor de dicha posición en la matriz. \n", "\n", "Sage utiliza el color negro para el mayor valor, el blanco para el más pequeño y tonos de grises para los intermedios. (Esto también depende del lenguaje/librería que usemos, usualmente se emplea el blanco para el valor más grande)." ] }, { "cell_type": "code", "execution_count": 68, "id": "60aab1c2", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n", "1 & 2 & 3 & 4 \\\\\n", "1 & 1 & 1 & 1 \\\\\n", "4 & 4 & 4 & 4 \\\\\n", "1 & 1 & 1 & 4\n", "\\end{array}\\right)$$" ], "text/plain": [ "[1 2 3 4]\n", "[1 1 1 1]\n", "[4 4 4 4]\n", "[1 1 1 4]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = matrix([[1,2,3,4],[1,1,1,1],[4,4,4,4],[1,1,1,4]])\n", "show(A)" ] }, { "cell_type": "code", "execution_count": 70, "id": "8782f013", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAHWCAYAAAD6lrl7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAT/ElEQVR4nO3dX4hc9fnA4XdMyCq4MxBsAiGjhBZKhhDBZCgrKsW2C3sh5q5XaUrTi+AakNypFy29WUEoLTgGpcU7iRQbFVpDF2o2iggbMVQMCIKQFE3VgjPrQjcYz+/Cbuj+stF3kpkzZ7PPA3Oxw+z5vmcPyYdz5l+tKIoiAIBvdNOoBwCAtUAwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwh+Tpp5+OHTt2xM033xx79uyJ119/fdQjrWunTp2KBx54ILZt2xa1Wi1eeumlUY+07s3MzES73Y7x8fHYsmVL7Nu3L95///1Rj7WuHT16NHbv3h31ej3q9XpMTEzEq6++OuqxKkMwh+CFF16IRx55JB5//PF455134t57742pqak4d+7cqEdbtxYXF+POO++Mp556atSj8F9zc3MxPT0db731VszOzsaXX34Zk5OTsbi4OOrR1q3t27fHE088EadPn47Tp0/H/fffHw8++GC89957ox6tEmo+fH3wfvCDH8Rdd90VR48evXzfzp07Y9++fTEzMzPCyYiIqNVqcfz48di3b9+oR+F/fPrpp7Fly5aYm5uL++67b9Tj8F+bN2+OJ598Mg4ePDjqUUbOGeaAXbx4Md5+++2YnJxccf/k5GS8+eabI5oKqq/b7UbE1/9BM3qXLl2KY8eOxeLiYkxMTIx6nErYOOoBbjSfffZZXLp0KbZu3bri/q1bt8aFCxdGNBVUW1EUceTIkbjnnnti165dox5nXXv33XdjYmIi/vOf/8Stt94ax48fj1arNeqxKkEwh6RWq634uSiKK+4Dvvbwww/HP/7xj3jjjTdGPcq69/3vfz/OnDkTn3/+ebz44otx4MCBmJubE80QzIG77bbbYsOGDVecTX7yySdXnHUCEYcPH45XXnklTp06Fdu3bx/1OOvepk2b4nvf+15EROzduzfm5+fj97//fTzzzDMjnmz0PIc5YJs2bYo9e/bE7OzsivtnZ2fj7rvvHtFUUD1FUcTDDz8cf/7zn+Pvf/977NixY9QjsYqiKGJpaWnUY1SCM8whOHLkSOzfvz/27t0bExMT8eyzz8a5c+fi0KFDox5t3friiy/igw8+uPzzhx9+GGfOnInNmzfH7bffPsLJ1q/p6el4/vnn4+WXX47x8fHLV2UajUbccsstI55ufXrsscdiamoqms1mLCwsxLFjx+LkyZNx4sSJUY9WDQVD0el0ijvuuKPYtGlTcddddxVzc3OjHmlde+2114qIuOJ24MCBUY+2bq12PCKieO6550Y92rr1i1/84vL/W9/5zneKH/3oR8Xf/va3UY9VGd6HCQAJnsMEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBHJKlpaX49a9/7SOlKsZxqR7HpHock9X54IIh6fV60Wg0otvtRr1eH/U4/JfjUj2OSfU4Jqtb02eYnU6n0tsbhqrvc9W3NwxV3+dh/A2rflzW49+w6sck4gbY59F+Mt/12blzZ2W31+12i4gout3uwLZZFNXe57WwvWEcl6rv86C3N+htrsdjMoxtVv2YFEW19zmj9G8r+eqrr+Kjjz6K8fHx6/5C5UuXLkWv1xvQZIPd3vJ2BjlfRLX3eS1sbxjHper7POjtDXqb6/GYDGObVT8mEdXc56IoYmFhIbZt2xY33fTNF11Lfw7zn//8ZzSbzTKXBIBvdP78+W/9AvPSzzDHx8cj4uvhPJkMwCj1er1oNpuX2/RNSg/m8mXYer0umABUQuYpwjX9KlkAKItgAkCCYAJAgmACQELpL/rJuN73ZzJ4hw4dGvUIrOJnP/vZqEfg/5mYmBj1CAyJM0wASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASCgtmJ1OJ1qtVrTb7bKWBICBKS2Y09PTcfbs2Zifny9rSQAYGJdkASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIKG0YHY6nWi1WtFut8taEgAGprRgTk9Px9mzZ2N+fr6sJQFgYFySBYAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSAhGsK5tNPPx07duyIm2++Ofbs2ROvv/76oOcCgErpO5gvvPBCPPLII/H444/HO++8E/fee29MTU3FuXPnhjEfAFRC38H87W9/GwcPHoxf/vKXsXPnzvjd734XzWYzjh49Ooz5AKAS+grmxYsX4+23347JyckV909OTsabb7656u8sLS1Fr9dbcQOAtaavYH722Wdx6dKl2Lp164r7t27dGhcuXFj1d2ZmZqLRaFy+NZvNa58WAEbkml70U6vVVvxcFMUV9y179NFHo9vtXr6dP3/+WpYEgJHa2M+Db7vtttiwYcMVZ5OffPLJFWedy8bGxmJsbOzaJwSACujrDHPTpk2xZ8+emJ2dXXH/7Oxs3H333QMdDACqpK8zzIiII0eOxP79+2Pv3r0xMTERzz77bJw7dy4OHTo0jPkAoBL6DuZPf/rT+Pe//x2/+c1v4uOPP45du3bFX//617jjjjuGMR8AVELfwYyIeOihh+Khhx4a9CwAUFk+SxYAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABJKC2an04lWqxXtdrusJQFgYEoL5vT0dJw9ezbm5+fLWhIABsYlWQBIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgoLZidTidarVa02+2ylgSAgSktmNPT03H27NmYn58va0kAGBiXZAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgoVYURVHmgr1eLxqNRnS73ajX66sPVauVORIAN7irpS7TpGXOMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgobRgdjqdaLVa0W63y1oSAAbGZ8kCcMPzWbIAUBLBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgITSgtnpdKLVakW73S5rSQAYmFpRFEWZC/Z6vWg0GtHtdqNer68+VK1W5kgA3OCulrpMk5a5JAsACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACX0H89SpU/HAAw/Etm3bolarxUsvvTSEsQCgWvoO5uLiYtx5553x1FNPDWMeAKikjf3+wtTUVExNTQ1jFgCoLM9hAkBC32eY/VpaWoqlpaXLP/d6vWEvCQADN/QzzJmZmWg0GpdvzWZz2EsCwMANPZiPPvpodLvdy7fz588Pe0kAGLihX5IdGxuLsbGxYS8DAEPVdzC/+OKL+OCDDy7//OGHH8aZM2di8+bNcfvttw90OACojKJPr732WhERV9wOHDiQ+v1ut1tERNHtdq/6mNW27+bm5ubmdq2362nSsr7PMH/4wx/G100DgPXD+zABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASChtGB2Op1otVrRbrfLWhIABqZWlPzllr1eLxqNRnS73ajX66sPVauVORIAN7irpS7TpGUuyQJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkBCacHsdDrRarWi3W6XtSQADEytKIqizAV7vV40Go3odrtRr9dXH6pWK3MkAG5wV0tdpknLXJIFgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBICEjaMeYDUlf5YCAHwrZ5gAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkFBaMDudTrRarWi322UtCQADUytK/uDWXq8XjUYjut1u1Ov1MpcGgBX6aZJLsgCQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJBQWjA7nU60Wq1ot9tlLQkAA1MriqIoc8FerxeNRiO63W7U6/UylwaAFfppkkuyAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQ0FcwZ2Zmot1ux/j4eGzZsiX27dsX77///rBmA4DK6CuYc3NzMT09HW+99VbMzs7Gl19+GZOTk7G4uDis+QCgEq7rC6Q//fTT2LJlS8zNzcV9992X+h1fIA1AVfTTpI3Xs1C3242IiM2bN1/1MUtLS7G0tLRiOABYa675RT9FUcSRI0finnvuiV27dl31cTMzM9FoNC7fms3mtS4JACNzzZdkp6en4y9/+Uu88cYbsX379qs+brUzzGaz6ZIsACM39Euyhw8fjldeeSVOnTr1jbGMiBgbG4uxsbFrWQYAKqOvYBZFEYcPH47jx4/HyZMnY8eOHcOaCwAqpa9gTk9Px/PPPx8vv/xyjI+Px4ULFyIiotFoxC233DKUAQGgCvp6DrNWq616/3PPPRc///nPU9vwthIAqmJoz2Fex1s2AWBN81myAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJBQWjA7nU60Wq1ot9tlLQkAA9PX92EOgu/DBKAq+mmSS7IAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkFBaMDudTrRarWi322UtCQADUyuKoihzwV6vF41GI7rdbtTr9TKXBoAV+mmSS7IAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkLBx1AOsplarjXoEgGtS8lvbKZEzTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIKC2YnU4nWq1WtNvtspYEgIGpFSV/8GGv14tGoxHdbjfq9frqQ/ksWWCN8lmya0umSctckgWABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBICE0oLZ6XSi1WpFu90ua0kAGJhaURRFmQv2er1oNBrR7XajXq+vPlStVuZIAANT8n+pXKdMk5a5JAsACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAl9BfPo0aOxe/fuqNfrUa/XY2JiIl599dVhzQYAldFXMLdv3x5PPPFEnD59Ok6fPh33339/PPjgg/Hee+8Naz4AqITr/gLpzZs3x5NPPhkHDx5MPd4XSAM3Ml8gvbb08wXSG691kUuXLsWf/vSnWFxcjImJias+bmlpKZaWllYMBwBrTd8v+nn33Xfj1ltvjbGxsTh06FAcP348Wq3WVR8/MzMTjUbj8q3ZbF7XwAAwCn1fkr148WKcO3cuPv/883jxxRfjD3/4Q8zNzV01mqudYTabTZdkgRuSS7JrSz+XZK/7Ocwf//jH8d3vfjeeeeaZgQ0nmMBaJZhrSz/BvO73YRZFseIMEgBuRH296Oexxx6LqampaDabsbCwEMeOHYuTJ0/GiRMnhjUfAFRCX8H817/+Ffv374+PP/44Go1G7N69O06cOBE/+clPhjUfAFRCX8H84x//OKw5AKDSfJYsACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAklBbMTqcTrVYr2u12WUsCwMDUiqIoylyw1+tFo9GIbrcb9Xp99aFqtTJHAhiYkv9L5TplmrTMJVkASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIKC2YnU4nWq1WtNvtspYEgIGpFUVRlLlgr9eLRqMR3W436vX66kPVamWOBDAwJf+XynXKNGmZS7IAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJCwcdQDrMYbfwGoGmeYAJAgmACQIJgAkCCYAJAgmACQIJgAkFD620qW3zLS6/XKXhoAVlhuUebtjKUHc2FhISIims1m2UsDwKoWFhai0Wh842NqRcmfEvDVV1/FRx99FOPj41Gr1a5rW+12O+bn5wc02WC31+v1otlsxvnz57/1W7z7UeV9XgvbG8Zxqfo+D3p7g97mejwmw9hm1Y9JRDX3uSiKWFhYiG3btsVNN33zs5Sln2HedNNNsX379oFsa8OGDQM9mIPeXkREvV6v9IzrbXvLBnlcqr7Pw/gbVv3fynr8G1b9mERUd5+/7cxy2Zp+0c/09HSltzcMVd/nqm9vGKq+z8P4G1b9uKzHv2HVj0nE2t/n0i/Jrhe9Xi8ajUZ0u92hnCFxbRyX6nFMqscxWd2aPsOssrGxsfjVr34VY2Njox6F/+G4VI9jUj2OyeqcYQJAgjNMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABI+D9RYLlwQGPeBgAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matrix_plot(-A)" ] }, { "cell_type": "markdown", "id": "2f8a521e", "metadata": {}, "source": [ "## Matrices dispersas\n", "Muchos sistemas lineales tienen gran parte de los coeficientes nulos. En ese tipo de sistemas, no merece la pena guardar todos esos ceros, sino la lista de valores no nulos. Para indicarle al ordenador que trabaje de ese modo, debemos indicar que la matriz es dispersa (sparse)." ] }, { "cell_type": "code", "execution_count": 71, "id": "08505bd0", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n", "1 & 0 & 0 & 1 \\\\\n", "0 & 2 & 0 & 0 \\\\\n", "1 & 0 & 3 & 0 \\\\\n", "0 & 0 & 0 & 4\n", "\\end{array}\\right)$$" ], "text/plain": [ "[1 0 0 1]\n", "[0 2 0 0]\n", "[1 0 3 0]\n", "[0 0 0 4]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = matrix(QQ,[[1,0,0,1],[0,2,0,0],[1,0,3,0],[0,0,0,4]],sparse=true)\n", "show(A)" ] }, { "cell_type": "markdown", "id": "3670f067", "metadata": {}, "source": [ "En ese caso, al representarla únicamente marca los valores no nulos, sin discriminar el valor. " ] }, { "cell_type": "code", "execution_count": 72, "id": "60f9c0da", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAHWCAYAAAD6lrl7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAWq0lEQVR4nO3dX4hc9fn48eckkhhiZkhIDYRda6jgz6msYHaQtVba2obmQgwEtTchbb2xNUJJL76NXrR3ayqlFDqGSIvtjSaCf28MDWiiRQRXDJGKhX5JyQa1Gkhn4pKumJzvRWvo/syfZ3TmzGz29YJhcTN7Ps/scec9c2ZnT1GWZRkAwAUtGvQAADAfCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCGafPPLII7Fu3bq4/PLLY/369fHKK68MeqQF7eWXX47bb7891q5dG0VRxLPPPjvokRa8ycnJaDabsWLFirjyyitj06ZN8de//nXQYy1ou3btirGxsajValGr1WJiYiJeeOGFQY81NASzD/bu3Rs/+clP4sEHH4w333wzvv71r8fGjRvj6NGjgx5twZqZmYkbbrghfvvb3w56FP7j4MGDcd9998Vrr70W+/fvj08++SQ2bNgQMzMzgx5twRoZGYmHHnoopqamYmpqKr71rW/FHXfcEX/5y18GPdpQKPzx9d676aab4sYbb4xdu3ad/dx1110XmzZtisnJyQFORkREURTxzDPPxKZNmwY9Cv/lww8/jCuvvDIOHjwYt95666DH4T9WrVoVDz/8cNxzzz2DHmXgPMPssY8//jjeeOON2LBhw5zPb9iwIV599dUBTQXDr91uR8S/76AZvNOnT8eePXtiZmYmJiYmBj3OULhs0ANcao4fPx6nT5+ONWvWzPn8mjVr4v333x/QVDDcyrKM7du3xy233BLXX3/9oMdZ0N56662YmJiIf/3rX3HFFVfEM888E41GY9BjDQXB7JOiKOb8d1mWn/kc8G/btm2Lw4cPx5///OdBj7LgXXvttXHo0KH45z//GU899VRs3bo1Dh48KJohmD23evXqWLx48WeeTX7wwQefedYJRNx///3x/PPPx8svvxwjIyODHmfBW7JkSVxzzTURETE+Ph6vv/56/OY3v4ndu3cPeLLB8xpmjy1ZsiTWr18f+/fvn/P5/fv3x8033zygqWD4lGUZ27Zti6effjpefPHFWLdu3aBH4hzKsozZ2dlBjzEUPMPsg+3bt8eWLVtifHw8JiYm4tFHH42jR4/GvffeO+jRFqyPPvoo/va3v5397yNHjsShQ4di1apVcdVVVw1wsoXrvvvui8cffzyee+65WLFixdmjMvV6PZYtWzbg6RamBx54IDZu3Bijo6Nx8uTJ2LNnTxw4cCD27ds36NGGQ0lftFqt8stf/nK5ZMmS8sYbbywPHjw46JEWtJdeeqmMiM9ctm7dOujRFqxz7Y+IKB977LFBj7Zg/fCHPzx7v/WlL32pvO2228o//elPgx5raHgfJgAkeA0TABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEs09mZ2fjF7/4hT8pNWTsl+Fjnwwf++Tc/OGCPul0OlGv16PdbketVhv0OPyH/TJ87JPhY5+c27x+htlqtYZ6e/0w7Ld52LfXD8N+m/vxPRz2/bIQv4fDvk8iLoHbPNi/zPfFXHfddUO7vXa7XUZE2W63e7bNshzu2zwftteP/TLst7nX2+v1NhfiPunHNod9n5TlcN/mjMrPVnLmzJl49913Y8WKFV/4hMqnT5+OTqfTo8l6u71Pt9PL+SKG+zbPh+31Y78M+23u9fZ6vc2FuE/6sc1h3ycRw3mby7KMkydPxtq1a2PRogsfdK38Ncxjx47F6OholUsCwAVNT09f9ATmlT/DXLFiRUT8ezgvJgMwSJ1OJ0ZHR8+26UIqD+anh2FrtZpgAjAUMi8RzuvfkgWAqggmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkVBbMVqsVjUYjms1mVUsCQM8UZVmWVS7Y6XSiXq9Hu912AmkABqqbJjkkCwAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmVBbPVakWj0Yhms1nVkgDQM0VZlmWVC3Y6najX69Fut6NWq1W5NADM0U2THJIFgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBICEzxXMRx55JNatWxeXX355rF+/Pl555ZVezzWvHTk+Ezv3vRP3P/Fm7Nz3Thw5PjPokQBS3H+dX9fnw9y7d29s2bIlHnnkkfja174Wu3fvjt/97nfx9ttvx1VXXXXRr7/Uz4f55NR0/Oypw1EURZRlefbjzs1jcef46KDHAzivhXj/1U2Tug7mTTfdFDfeeGPs2rXr7Oeuu+662LRpU0xOTvZ0uPnmyPGZuO1XB+LMOb6ji4qIF3/6jbh69fLqBwO4iIV6/9W3E0h//PHH8cYbb8SGDRvmfH7Dhg3x6quvnvNrZmdno9PpzLlcqp6cmo6iKM75b0VRxN6p6YonAshx/3VxXQXz+PHjcfr06VizZs2cz69Zsybef//9c37N5ORk1Ov1s5fR0UvzaX1ExLETp+J8T9jLsoxjJ05VPBFAjvuvi/tcv/Tz/z8K+fRY97ns2LEj2u322cv09KX7KGVk5bILPkIbWbms4okActx/XVxXwVy9enUsXrz4M88mP/jgg8886/zU0qVLo1arzblcqu4aH73gI7S7L9EXzYH5z/3XxXUVzCVLlsT69etj//79cz6/f//+uPnmm3s62Hy0bvXy2Ll5LBYVEYsXFXM+7tw8dkm+YA5cGtx/Xdxl3X7B9u3bY8uWLTE+Ph4TExPx6KOPxtGjR+Pee+/tx3zzzp3jo9G8elXsnZqOYydOxcjKZXH3+Kj/2YCh5/7rwrp+W0nEv/9wwS9/+ct477334vrrr49f//rXceutt6a+9lJ+WwkA80tf34f5RQkmAMOib+/DBICFSjABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIKGyYLZarWg0GtFsNqtaEgB6xvkwAViwnA8TAHpMMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIKGyYLZarWg0GtFsNqtaEgB6pijLsqxywU6nE/V6PdrtdtRqtSqXBoA5ummSQ7IAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQUFkwW61WNBqNaDabVS0JAD1TlGVZVrlgp9OJer0e7XY7arValUsDwBzdNMkhWQBIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIqCyYrVYrGo1GNJvNqpYEgJ4pyrIsq1yw0+lEvV6PdrsdtVqtyqUBYI5umuSQLAAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkdB3Ml19+OW6//fZYu3ZtFEURzz77bB/GAoDhclm3XzAzMxM33HBD/OAHP4jNmzf3YyboiyPHZ+LJqek4duJUjKxcFneNj8a61csHPRYwT3QdzI0bN8bGjRv7MQv0zZNT0/Gzpw5HURRRlmUURRG7D/5v7Nw8FneOjw56PGAe8Boml7wjx2fiZ08djjNlxOkz5ZyP//PU4fj78ZlBjwjMA30P5uzsbHQ6nTkXqNKTU9NRFMU5/60oitg7NV3xRMB81PdgTk5ORr1eP3sZHXX4i2odO3EqyrI857+VZRnHTpyqeCJgPup7MHfs2BHtdvvsZXrao3mqNbJy2QWfYY6sXFbxRMB81PdgLl26NGq12pwLVOmu8dELPsO82y/9AAldB/Ojjz6KQ4cOxaFDhyIi4siRI3Ho0KE4evRor2eDnli3enns3DwWi4qIxYuKOR93bh6Lq721BEgoyvM99D6PAwcOxDe/+c3PfH7r1q3xhz/84aJf3+l0ol6vR7vd9myTSv39+Ezs/a/3Yd49PiqWsMB106Sug/lFCSYAw6KbJnkfJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJFQWzFarFY1GI5rNZlVLAkDPOB8mAAuW82ECQI8JJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJFQWzFarFY1GI5rNZlVLAkDPFGVZllUu2Ol0ol6vR7vdjlqtVuXSADBHN01ySBYAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEioLZqvVikajEc1ms6olAaBnirIsyyoX7HQ6Ua/Xo91uR61Wq3JpAJijmyY5JAsACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJlQWz1WpFo9GIZrNZ1ZIA0DNFWZZllQt2Op2o1+vRbrejVqtVuTQAzNFNkxySBYAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSAhMu6ufLk5GQ8/fTT8c4778SyZcvi5ptvjp07d8a1117br/nmpSPHZ+LJqek4duJUjKxcFneNj8a61csHPRYMHT8rzCddnQ/zu9/9bnzve9+LZrMZn3zySTz44IPx1ltvxdtvvx3Ll+f+J7/Uz4f55NR0/Oypw1EURZRlefbjzs1jcef46KDHg6HhZ4Vh0E2TvtAJpD/88MO48sor4+DBg3Hrrbf2fLj55sjxmbjtVwfizDm+o4uKiBd/+o242qNn8LPC0KjsBNLtdjsiIlatWnXe68zOzkan05lzuVQ9OTUdRVGc89+Kooi9U9MVTwTDyc8K89HnDmZZlrF9+/a45ZZb4vrrrz/v9SYnJ6Ner5+9jI5euodajp04Fed7wl6WZRw7cariiWA4+VlhPvrcwdy2bVscPnw4nnjiiQteb8eOHdFut89epqcv3UeOIyuXXfBR88jKZRVPBMPJzwrz0ecK5v333x/PP/98vPTSSzEyMnLB6y5dujRqtdqcy6XqrvHRCz5qvtsvMkBE+FlhfuoqmGVZxrZt2+Lpp5+OF198MdatW9evuealdauXx87NY7GoiFi8qJjzcefmMb/EAP/hZ4X5qKvfkv3xj38cjz/+eDz33HNz3ntZr9dj2bLcIZRL+bdkP/X34zOx97/eW3b3+Kg7ADgHPysMWt/eVnK+1xwee+yx+P73v9/z4QCgn7ppUld/6ecLvGUTAOY1f0sWABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIqC2ar1YpGoxHNZrOqJQGgZ7o6H2YvOB8mAMOimyY5JAsACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJggkACZUFs9VqRaPRiGazWdWSANAzRVmWZZULdjqdqNfr0W63o1arVbk0AMzRTZMckgWABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBIAEwQSABMEEgATBBICEyoLZarWi0WhEs9msakkA6JmiLMuyygU7nU7U6/Vot9tRq9WqXBoA5uimSQ7JAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQEJlwWy1WtFoNKLZbFa1JAD0TFGWZVnlgp1OJ+r1erTb7ajValUuDQBzdNMkh2QBIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASChq2Du2rUrxsbGolarRa1Wi4mJiXjhhRf6NRsAFTtyfCZ27nsn7n/izdi57504cnxm0CMNjcu6ufLIyEg89NBDcc0110RExB//+Me444474s0334yvfvWrfRkQgGo8OTUdP3vqcBRFEWVZRlEUsfvg/8bOzWNx5/jooMcbuC98AulVq1bFww8/HPfcc0/q+k4gDTB8jhyfidt+dSDOnKMIi4qIF3/6jbh69fLqB+uzSk4gffr06dizZ0/MzMzExMTEea83OzsbnU5nzgWA4fLk1HQURXHOfyuKIvZOTVc80fDpOphvvfVWXHHFFbF06dK4995745lnnolGo3He609OTka9Xj97GR31tB5g2Bw7cSrOd8CxLMs4duJUxRMNn66Dee2118ahQ4fitddeix/96EexdevWePvtt897/R07dkS73T57mZ72KAVg2IysXHbBZ5gjK5dVPNHw+cKvYX7729+Or3zlK7F79+7U9b2GCTB8vIbZx9cwP1WWZczOzn7RzQAwQOtWL4+dm8diURGxeFEx5+POzWOXZCy71dXbSh544IHYuHFjjI6OxsmTJ2PPnj1x4MCB2LdvX7/mA6Aid46PRvPqVbF3ajqOnTgVIyuXxd3jo2L5H10F8x//+Eds2bIl3nvvvajX6zE2Nhb79u2L73znO/2aD4AKXb16efzPd//foMcYSl0F8/e//32/5gCAoeZvyQJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQEJlwWy1WtFoNKLZbFa1JAD0TFGWZVnlgp1OJ+r1erTb7ajValUuDQBzdNMkh2QBIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgQTABIEEwASBBMAEgobJgtlqtaDQa0Ww2q1oSAHqmKMuyrHLBTqcT9Xo92u121Gq1KpcGgDm6aZJDsgCQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJAgmACQIJgAkCCYAJBwWdULfnr6zU6nU/XSADDHpy3KnBq68mCePHkyIiJGR0erXhoAzunkyZNRr9cveJ2izGS1h86cORPvvvturFixIoqi+ELbajab8frrr/dost5ur9PpxOjoaExPT1/0LN7dGObbPB+214/9Muy3udfb6/U2F+I+6cc2h32fRAznbS7LMk6ePBlr166NRYsu/Cpl5c8wFy1aFCMjIz3Z1uLFi3u6M3u9vYiIWq021DMutO19qpf7Zdhvcz++h8P+s7IQv4fDvk8ihvc2X+yZ5afm9S/93HfffUO9vX4Y9ts87Nvrh2G/zf34Hg77flmI38Nh3ycR8/82V35IdqHodDpRr9ej3W735RkSn4/9Mnzsk+Fjn5zbvH6GOcyWLl0aP//5z2Pp0qWDHoX/Yr8MH/tk+Ngn5+YZJgAkeIYJAAmCCQAJggkACYIJAAmCCQAJggkACYIJAAmCCQAJ/wdIj1eQHvlwNQAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "matrix_plot(A)" ] }, { "cell_type": "markdown", "id": "d09ac686", "metadata": {}, "source": [ "
\n", " Ejercicio 2. \n", "Consideremos la matriz\n", "$$A=\\left(\\begin{array}{rrrrr} \n", " 4 & -1 & 0 & 0 & -1 \\\\\n", " 0 & 4 & -1 & -1 & 0 \\\\\n", "-1 & 0 & 4 & 0 & 0 \\\\\n", "-1 & 0 & 0 & 4 & 0 \\\\\n", "-1 & 0 & 0 & 0 & 4 \\\\\n", "\\end{array}\\right)$$\n", "\n", "a) Calcula la factorización LU de $A$. Cuenta cuántos valores no nulos hay en la matriz $A$ y cuántos hay en las matrices $L$ y $U$. \n", " \n", "b) Aplica tres pasos del método de Jacobi para resolver el sistema $Ax = b$ con $b=(1,1,1,1,1)$,\n", " usando matrices dispersas. \n", " \n", "c) Idem para el método de Gauss-Seidel. \n", "
" ] }, { "cell_type": "code", "execution_count": 1, "id": "8e42302f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrrr}\n", "4 & -1 & 0 & 0 & -1 \\\\\n", "0 & 4 & -1 & -1 & 0 \\\\\n", "-1 & 0 & 4 & 0 & 0 \\\\\n", "-1 & 0 & 0 & 4 & 0 \\\\\n", "-1 & 0 & 0 & 0 & 4\n", "\\end{array}\\right)$$" ], "text/plain": [ "[ 4 -1 0 0 -1]\n", "[ 0 4 -1 -1 0]\n", "[-1 0 4 0 0]\n", "[-1 0 0 4 0]\n", "[-1 0 0 0 4]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "A = matrix([[4,-1,0,0,-1],\n", " [0,4,-1,-1,0],\n", " [-1,0,4,0,0],\n", " [-1,0,0,4,0],\n", " [-1,0,0,0,4]])\n", "show(A)" ] }, { "cell_type": "code", "execution_count": 2, "id": "12313435", "metadata": {}, "outputs": [], "source": [ "def es_no_nulo(k):\n", " if k ==0:\n", " return 0\n", " return(1)" ] }, { "cell_type": "code", "execution_count": 3, "id": "2cf63635", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum([sum([es_no_nulo(k) for k in fila]) for fila in A])" ] }, { "cell_type": "code", "execution_count": 12, "id": "49dd45d1", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1000 x 1000 sparse matrix over Integer Ring (use the '.str()' method to see the entries)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = 4*matrix.identity(1000,sparse=true)\n", "B[1:,0] = -1\n", "B[0,1:] = -1\n", "B" ] }, { "cell_type": "code", "execution_count": 13, "id": "fdebb751", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2998" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum([sum([es_no_nulo(k) for k in fila]) for fila in B])" ] }, { "cell_type": "code", "execution_count": 5, "id": "b0ae62fa", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(\n", "[ 1 0 0 0 0] [ 4 -1 0 0 -1]\n", "[ 0 1 0 0 0] [ 0 4 -1 -1 0]\n", "[ -1/4 -1/16 1 0 0] [ 0 0 63/16 -1/16 -1/4]\n", "[ -1/4 -1/16 -1/63 1 0] [ 0 0 0 248/63 -16/63]\n", "[ -1/4 -1/16 -1/63 -1/62 1], [ 0 0 0 0 116/31]\n", ")" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "_,L,U = A.LU()\n", "L,U" ] }, { "cell_type": "code", "execution_count": 6, "id": "5fee5a39", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "14" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum([sum([es_no_nulo(k) for k in fila]) for fila in L])" ] }, { "cell_type": "code", "execution_count": 7, "id": "221ec644", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum([sum([es_no_nulo(k) for k in fila]) for fila in U])" ] }, { "cell_type": "code", "execution_count": null, "id": "a3eac119", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "fd1a3c00", "metadata": {}, "source": [ "## Difuminando y perfilando una imagen (opcional) \n", "\n", "Vamos a ver un ejemplo de cómo se puede aplicar el método iterativo de Gauss-Seidel. Concretamente, veremos cómo funciona (de modo básico) un modelo de difuminado y perfilado. Estos modelos se aplican por ejemplo a imágenes de satélite para eliminar los efectos de difusión de la atmósfera, al análisis de imágenes por ordenador, para retocar fotografías, etc.\n", "\n", "Cada imagen viene cargada en una matriz (por simplicidad asumiremos que está en blanco y negro). Vamos a necesitar algunas funciones que nos pasen los datos de la matriz a un vector y viceversa. Para ello, basta ejecutar la siguiente celda." ] }, { "cell_type": "code", "execution_count": null, "id": "e81a2b19", "metadata": {}, "outputs": [], "source": [ "def matrix_to_vector(A):\n", " n,m=A.dimensions();\n", " v=vector(RDF,n*m)\n", " for i in range(0,n):\n", " for j in range(0,m):\n", " v[m*i+j]=A[i,j];\n", " return v;\n", "def vector_to_matrix(v):\n", " n=sqrt(v.degree());\n", " A=matrix(RDF,n);\n", " for i in range(0,n):\n", " for j in range(0,n):\n", " A[i,j]=v[n*i+j];\n", " return A;" ] }, { "cell_type": "markdown", "id": "298b5008", "metadata": {}, "source": [ "Para pasar la matriz a un vector, numeramos las celdas comenzando por la primera fila:\n", "\n", "$\\left(\\begin{array}{rrr}\n", "1 & 2 & 3 \\\\\n", "4 & 5 & 6 \\\\\n", "7 & 8 & 9\n", "\\end{array}\\right)\\to(1,2,3,4,5,6,7,8,9)$\n", "\n", "Definimos la matriz en la que tenemos la imagen. Cuando creamos la matriz estará a cero (color negro). Vamos a poner algunas posiciones a uno (blanco)." ] }, { "cell_type": "code", "execution_count": null, "id": "37d9ab7f", "metadata": {}, "outputs": [], "source": [ "n=5\n", "B = matrix.zero(RDF,n)\n", "B[1,1:4]=Matrix([1,1,1]) \n", "B[3,1:4]=Matrix([1,1,1]) \n", "B[2,1]=1\n", "show(matrix_plot(B), figsize=5)" ] }, { "cell_type": "markdown", "id": "d93a1ce8", "metadata": {}, "source": [ "Creamos la matriz de difuminado (está libremente basado en el difuminado gaussiano). Para ello cada posición pasa a tomar \"un poco del color de los vecinos\". Estamos considerando cuatro vecinos (arriba, abajo, derecha e izquierda). Si A(i,j) es el color del punto (i,j), después de la transformación será: $$A(i,j)\\to (1-4d) A(i,j)+d\\,A(i-1,j)+d\\,A(i+1,j)+d\\,A(i,j-1)+d\\,A(i,j+1)$$\n", "\n", "Esto puede hacerse escribiendo la matriz $A$ en forma de vector y multiplicando el vector por cierta matriz (dada por la expresión anterior). Guardaremos en $N$ la matriz por la que hay que multiplicar (matriz de difuminado)." ] }, { "cell_type": "code", "execution_count": null, "id": "7c457499", "metadata": {}, "outputs": [], "source": [ "A=matrix(RDF,n*n,sparse=true);\n", "d=1/10;\n", "for i in range(0,n-1):\n", " for j in range(0,n):\n", " A[i+j*n,i+1+j*n]=d; # vecino a la derecha (i+1,j)\n", " A[i+1+j*n,i+j*n]=d; # vecino a la izquierda (i-1,j)\n", " A[i*n+j,(i+1)*n+j]=d; # vecino superior (i,j+1)\n", " A[(i+1)*n+j,i*n+j]=d; # vecino inferior (i,j-1)\n", "for i in range(0,n*n):\n", " A[i,i]=1-4*d; # posición original (i,j)\n", "matrix_plot(A)" ] }, { "cell_type": "markdown", "id": "34a680a3", "metadata": {}, "source": [ "Para difuminar la imagen, escribimos la matriz como vector y multiplicamos ese vector por la matriz de difuminado" ] }, { "cell_type": "code", "execution_count": null, "id": "37e9e05c", "metadata": {}, "outputs": [], "source": [ "Bv=matrix_to_vector(B) # Necesitamos considerar la imagen como un vector\n", "b=A*Bv # difuminamos\n", "show(matrix_plot(vector_to_matrix(b)), figsize=5) # Para dibujar necesitamos volver a pasar a la forma de matriz" ] }, { "cell_type": "markdown", "id": "fc2d6ab1", "metadata": {}, "source": [ "Y podemos hacer el proceso inverso resolviendo el sistema $Nx=b$, para lo cual Gauss-Seidel es un método apropiado (a cada paso se perfila un poco más la imagen). Véamoslo en el ejemplo.\n", "\n", "Construimos las matrices de Gauss-Seidel" ] }, { "cell_type": "code", "execution_count": null, "id": "5925f1db", "metadata": {}, "outputs": [], "source": [ "N = matrix(RDF,n*n);\n", "for i in range(0,n*n):\n", " for j in range(i+1,n*n):\n", " N[i,j] =- A[i,j];\n", "M = A + N;" ] }, { "cell_type": "markdown", "id": "3d0ea2dd", "metadata": {}, "source": [ "Y aplicamos el método" ] }, { "cell_type": "code", "execution_count": null, "id": "f7453d2b", "metadata": {}, "outputs": [], "source": [ "x0 = vector(RDF,n*n)\n", "m = 3 # Número de iteraciones del método\n", "for i in range(1,m+1):\n", " x0 = M.solve_right(N*x0+b)\n", "matrix_plot(vector_to_matrix(vector(x0)))" ] }, { "cell_type": "code", "execution_count": null, "id": "45c41e22", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.5", "language": "sage", "name": "sagemath" }, "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.11.2" } }, "nbformat": 4, "nbformat_minor": 5 }