{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy import sparse\n", "from scipy.sparse.linalg import dsolve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Equação de flexura por diferenças finitas (versão 2)\n", "\n", "A matriz $\\textbf A$ é um tipo de matriz esparsa, ou seja, a maior parte dos elementos da matriz são nulos. Para observar essa característica, vamos gerar a matriz $\\textbf A$ e visualizar os elementos da matriz que são diferentes de zero. Nesse exercício as 5 diagonais foram preenchidas com 1 (apenas para indicar um valor diferente de zero) e usou-se a função `imshow` do `matplotlib` para possibilitar a visualização dos elementos não-nulos." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD7CAYAAABOrvnfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC+tJREFUeJzt3V/IZAUZx/Hvk2sXZmBS7S6yNV1URAgugQQWvUSFEZje\nGEKwhIUXZdJV2kVqXVSCIt0I4RprhRWJm11UWvRuQvRHWXVLywIHNNx3jYx27yyeLuZszr573pnX\n+XfO7vP9wLBnz5yZeTjsb8+Z+b1n3shMJNXwmq4HkLQ6Bl4qxMBLhRh4qRADLxVi4KVCVhL4iLg8\nIv4cEX+NiC+u4jWniYhhRDwZEYcj4vcdzXBPRGxExJGxdRdGxMMR8UxEPBQRF/Rgplsi4vlmXx2O\niMtXPNOeiPhVRPwpIv4YEZ9v1ne9r7aaq9P9NUksu4ePiHOAvwAfAv4O/AG4JjOfXuoLT5/rWeA9\nmfnPDmd4P3ACuDczL27W3Qb8IzNva/5zfENm3tjxTDcDxzPzjlXNsWmmXcCuzHw8Is4HHgOuBD5F\nt/tqq7mupsP9NckqjvCXAn/LzGFmvgx8H/j4Cl53O6LLF8/MR4CXNq2+AjjQLB9g9A+o65mgw32V\nmUcz8/Fm+QTwNHAR3e+rreaCjv9tbWUVgb8IeG7s78/zyk7pUgK/iIhHI+IzXQ8zZmdmbjTLG8DO\nLocZc31EPBER+1d96jwuIgbAXuB39Ghfjc3122ZVL/bXZqsIfF9/dveyzNwLfBT4bHMq2ys5er/V\nh/13F/A24BLgBeD2LoZoTpvvB27IzOPj93W5r5q5ftTMdYKe7K82qwj834E9Y3/fw+go36nMfKH5\n80XgAUZvPfpgo3lvSETsBo51PA+ZeSwbwN10sK8i4lxGYf9OZh5sVne+r8bm+u7Jufqwv7ayisA/\nCrw9IgYR8VrgE8CDK3jdLUXEeRHx+mb5dcBHgCOTH7UyDwL7muV9wMEJ265EE6aTrmLF+yoiAtgP\nPJWZd47d1em+2mqurvfXJEv/lB4gIj4K3AmcA+zPzK8t/UUnz/M2Rkd1gB3A97qYKSLuAz4AvJHR\ne9AvAz8Gfgi8BRgCV2fmvzqc6WZgjdHpaQLPAteNvXdexUzvA34NPMkrp+03Ab+n233VNteXgGvo\ncH9NspLAS+oHf9JOKsTAS4UYeKkQAy8VMnPg+3hBjKTJZvqUfjsXxESEH/9LHcrM036ef8eMz/X/\nC2IAIuLkBTGnXAF389jyOqMyd163nvKs81pnMVMt2jr9m2sdZ9qOdfox062ta2c9pe/rBTGSJpj1\nCL+t0/X1seWV/fiTVNKwuU02a+C3dUHM2qZx+mfQ9QBbGHQ9QItB1wO0GHQ9QItBh687/tqHWrea\n9ZT+VV8QM5h0Z2cGXQ+whUHXA7QYdD1Ai0HXA7QYdD3ARDMd4TPzPxHxOeDnvHJBzEq+surmLT6M\nOGmxH+pJZ5dZT+nJzJ8CP13gLJKWzJ+0kwox8FIhBl4qxMBLhRh4qRADLxUycy3XV5N6ejt6VecR\nXirEwEuFGHipEAMvFWLgpUIMvFTI0n7V1OhLLNtrsGmXuHbF2k5nj1tbv8TSI7xUiIGXCjHwUiEG\nXirEwEuFGHipEAMvFdJJDz9NH3t6O3qdWezhpfIMvFSIgZcKMfBSIQZeKsTAS4XMVctFxBD4N/Bf\n4OXMvHTsvplruUn6WNmBtZ36pr2Wm/drqhNYy8x/zvk8klZgEaf0p/0vIqmf5g18Ar+IiEcj4jOL\nGEjS8sx7Sn9ZZr4QEW8CHo6IP2fmI4sYTNLizRX4zHyh+fPFiHgAuBQYC/z62NaD5iZp8YbNbbKZ\nAx8R5wHnZObxiHgd8BHY/BH62qxPL+lVGXDqAfVQ61bzHOF3Ag9ExMnn+V5mPjTH80lasl5eHjuP\nPvb0dvRaPS+Plcoz8FIhBl4qxMBLhRh4qRADLxVy1tVyk/SxsgNrOy2DtZxUnoGXCjHwUiEGXirE\nwEuFGHipEAMvFVKqh5+mjz29Hb1mYw8vlWfgpUIMvFSIgZcKMfBSIQZeKsRabpv6WNmBtZ22Yi0n\nlWfgpUIMvFSIgZcKMfBSIQZeKmRq4CPinojYiIgjY+sujIiHI+KZiHgoIi5Y7piSFmFqDx8R7wdO\nAPdm5sXNutuAf2TmbRHxReANmXnjpsedVT38NH3s6e3oK5uxh8/MR4CXNq2+AjjQLB8Arpx7PklL\nN+t7+J2ZudEsbwA7FzSPpCWa+0O7HL0nWM7P50paqB0zPm4jInZl5tGI2A0ca99sfWx50NwkLd6w\nuU02a+AfBPYB32j+PNi+2dqMTy/p1Rlw6gH1UOtW26nl7gN+A7wzIp6LiE8BXwc+HBHPAB9s/i6p\n57w8dgX6WNmBtd3ZzctjpfIMvFSIgZcKMfBSIQZeKsTAS4UYeKkQe/ge6GNPb0d/prOHl8oz8FIh\nBl4qxMBLhRh4qRADLxViLddzfazswNqu/6zlpPIMvFSIgZcKMfBSIQZeKsTAS4UYeKkQe/gzXB97\nejv6PrCHl8oz8FIhBl4qxMBLhRh4qRADLxUytZaLiHuAjwHHMvPiZt0twKeBF5vNbsrMn216nLVc\nx/pY2YG13WrMXst9G7h807oE7sjMvc3tZy2Pk9QzUwOfmY8AL7Xcddr/HpL6bZ738NdHxBMRsT8i\nLljYRJKWZseMj7sL+Eqz/FXgduDa0zdbH1seNDdJizdsbpPNFPjMPHZyOSLuBn7SvuXaLE8v6VUb\ncOoB9VDrVjOd0kfE7rG/XgUcmeV5JK3Wdmq5+4APAG8ENhh1bWvAJYw+rX8WuC4zNzY9zlqu5/pY\n21nZLUp7LTf1lD4zr2lZfc9CZpK0Uv6knVSIgZcKMfBSIQZeKsTAS4UYeKkQv7VWrfrY0YM9/fb5\nrbVSeQZeKsTAS4UYeKkQAy8VYuClQqzlNJM+1nZWduOs5aTyDLxUiIGXCjHwUiEGXirEwEuFGHip\nEHt4LVwfO3qo1tPbw0vlGXipEAMvFWLgpUIMvFSIgZcKmVjLRcQe4F7gzYx+U+y3MvObEXEh8APg\nrYx+C/3VmfmvTY+1llOrPtZ2Z19lN1st9zLwhcx8N/Be4LMR8S7gRuDhzHwH8Mvm75J6bmLgM/No\nZj7eLJ8AngYuAq4ADjSbHQCuXOaQkhZj2+/hI2IA7AV+B+zMzI3mrg1g58Ink7RwO7azUUScD9wP\n3JCZxyNeeWuQmTl6v95mfWx50NwkLd6wuU02NfARcS6jsH8nMw82qzciYldmHo2I3cCx9kevbWtU\nSfMacOoB9VDrVhNP6WN0KN8PPJWZd47d9SCwr1neBxzc/FhJ/TPtCH8Z8EngyYg43Ky7Cfg68MOI\nuJamllvahJIWxstj1St97OjhTOzpvTxWKs/AS4UYeKkQAy8VYuClQgy8VIi1nM4ofazt+lnZWctJ\n5Rl4qRADLxVi4KVCDLxUiIGXCjHwUiH28Dpr9LGjh656ent4qTwDLxVi4KVCDLxUiIGXCjHwUiHW\nciqjj7Xd8io7azmpPAMvFWLgpUIMvFSIgZcKMfBSIdN+XfSeiPhVRPwpIv4YEZ9v1t8SEc9HxOHm\ndvlqxpU0j4k9fETsAnZl5uMRcT7wGHAlo18PfTwz75jwWHt4nTH62NHDPD19ew8/8ffDZ+ZR4Giz\nfCIingYuau4+7ckk9du238NHxADYC/y2WXV9RDwREfsj4oIlzCZpwSYe4U9qTud/BNzQHOnvAr7S\n3P1V4Hbg2tMfuT62PGhukhZv2Nwmmxr4iDgXuB/4bmYeBMjMY2P33w38pP3Ra9PnlLQAA049oB5q\n3Wrap/QB7Aeeysw7x9bvHtvsKuDIjFNKWqFpR/jLgE8CT0bE4Wbdl4BrIuISIIFngeuWN6KkRfHy\nWGkb+ljbTa7svDxWKs/AS4UYeKkQAy8VYuClQgy8VIiBlwqxh5fm1M+OHnt4qToDLxVi4KVCDLxU\niIGXCjHwUiHb+oorSVub9s2yfartPMJLhRh4qRADLxVi4KVCDLxUiIGXCjHwUiH28NKSTerpV93R\ne4SXCjHwUiErDPxwdS+1bcOuB9jCsOsBWgy7HqDFsOsBWgy7HmAiA99Lw64HaDHseoAWw64HaDHs\neoCJPKWXCjHwUiFL/tZaSV1p+9bapQVeUv94Si8VYuClQgy8VIiBlwox8FIh/wOx9L53JNy0aQAA\nAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 30\n", "\n", "A = np.zeros((N,N))\n", "\n", "A[range(N),range(N)]= 1\n", "A[range(0,N-1),range(1,N)]=1\n", "A[range(1,N),range(0,N-1)]=1\n", "A[range(0,N-2),range(2,N)]=1\n", "A[range(2,N),range(0,N-2)]=1\n", "\n", "plt.matshow(A, interpolation='nearest')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nesta figura os elementos não-nulos são representados por quadrados vermelhos. Podemos ver que a maior parte dos elementos são nulos (região em azul). O efeito fica ainda mais evidente com o aumento do número de pontos $N$. (_Verifique!_)\n", "\n", "Assim há um grande desperdício de memória para alocar a matriz $\\textbf A$. Adicionalmente, a eficiência da solução do sistema linear fica comprometida, ou seja, o algoritmo gasta um tempo significativo multiplicando e somando termos nulos que não influenciam no resultado final. \n", "\n", "Podemos melhorar a eficiência do algoritmo alocando apenas os termos não-nulos da matriz $\\textbf A$. Para isso vamos ter que importar alguns métodos para manipulações algébricas de matrizes esparsas da biblioteca Scipy, como indicado no início deste notebook:\n", "\n", "```Python\n", "from scipy import sparse\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para alocar a matriz esparsa, teremos que criar 3 vetores auxiliares: um para armazenar os valores dos elementos não-nulos e outros dois para armazenar os índices $(i,j)$ dos respectivos elementos.\n", "\n", "Por exemplo, considere a seguinte matriz esparsa (onde os pontos representam os elementos nulos, indicados desta forma para facilitar a visualização dos termos não-nulos) \n", "$$\n", "\\begin{bmatrix}\n", "1 & 2 & . & . \\\\\n", ". & 3 & 4 & 5 \\\\\n", ". & . & 6 & . \\\\\n", ". & . & . & 7 \\\\\n", "\\end{bmatrix}\n", "$$ \n", "\n", "Para representá-la de forma compacta vamos criar um vetor `vals` para armazenar os valores e dois vetores `indi` e `indj` para armazenar os índices das linhas e das colunas, respectivamente. O resultado seria o seguinte:\n", "\n", "```Python\n", "vals = [1, 2, 3, 4, 5, 6, 7]\n", "indi = [0, 0, 1, 1, 1, 2, 3]\n", "indj = [0, 1, 1, 2, 3, 2, 3]\n", "```\n", "Para criar a matriz esparsa usamos a classe `csr_matrix` (_Compressed Sparse Row matrix_):\n", "```Python\n", "ni,nj = 4,4\n", "Asp = sparse.csr_matrix((vals,(indi,indj)), shape=(ni,nj))\n", "```\n", "onde `ni` e `nj` são respectivamente o número de linhas e colunas da matriz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Podemos obter a matriz \"cheia\" original, com os zeros, usando-se o método `todense()`:\n", "```Python\n", "Aoriginal = sparse.csr_matrix((vals,(indi,indj)), shape=(ni,nj)).todense()\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A matriz esparsa pode ser multiplicada com um vetor ou matriz usando-se o método `dot()`:\n", "```Python\n", "v = [1,1,1,1]\n", "b = Asp.dot(v)\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Para resolver o sistema linear `Asp v = b` usa-se o método `spsolve`:\n", "```Python\n", "v = dsolve.spsolve(Asp,b)\n", "```\n", "e o seguinte módulo deve ser importado:\n", "```Python\n", "from scipy.sparse.linalg import dsolve\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.7.4" } }, "nbformat": 4, "nbformat_minor": 1 }