Skip to content

Instantly share code, notes, and snippets.

@patonelli

patonelli/03-SistemasLineares-II

Forked from anonymous/file1.py
Last active August 29, 2015 13:57
Show Gist options
  • Save patonelli/9746377 to your computer and use it in GitHub Desktop.
Save patonelli/9746377 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
03-SistemasLineares-II
{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exemplos do M\u00e9todo da Elimina\u00e7\u00e3o de Gauss\n",
"\n",
"Veremos como funciona o algoritmo nas situa\u00e7\u00f5es mais simples, usando a pivota\u00e7\u00e3o para estabilizar os c\u00e1lculos.\n",
"\n",
"## Um sistema simples $Ax=b$.\n",
"Vamos considerar o exemplo simples:\n",
"$$ \\begin{pmatrix}3 & 2& 1 \\\\\n",
"1 & 4 & 2 \\\\\n",
"6 & 1& 0 \n",
"\\end{pmatrix}\\begin{pmatrix} x \\\\ y \\\\ z \\end{pmatrix}=\\begin{pmatrix} 1 \\\\ 2 \\\\ 1\\end{pmatrix}\n",
"$$\n",
"Para executar as opera\u00e7\u00f5es elementares, consideremos as matrizes $A$ e $b$ agregadas da forma $[A|b]$. Aplica-se o algoritmo \u00e0 matriz\n",
"$$ \\begin{bmatrix}3 & 2& 1 & 1\\\\\n",
"1 & 4 & 2& 2 \\\\\n",
"6 & 1& 0 &1 \n",
"\\end{bmatrix}$$\n",
"\n",
"1. Pivota\u00e7\u00e3o na coluna 1 transforma o sistema em: $$\\begin{bmatrix}6 & 1& 0 & 1\\\\ 1 & 4 & 2& 2 \\\\ 3 & 2& 1 &1 \\end{bmatrix}$$\n",
"\n",
"2. A seguir aplicamos o primeiro passo da elimina\u00e7\u00e3o de Gauss. Obtemos a matriz $$\\begin{bmatrix}6 & 1& 0 & 1\\\\ 0 & 3.833 & 2& 1.833 \\\\ 0 & 1.5& 1 &0.5 \\end{bmatrix}$$\n",
"\n",
"3. Agora novamente a pivota\u00e7\u00e3o e o segundo passo do m\u00e9todo da elimina\u00e7\u00e3o temos a matriz $$\\begin{bmatrix}6 & 1& 0 & 1\\\\ 0 & 3.833 & 2& 1.833 \\\\ 0 & 0& 0.271 & -0.271 \\end{bmatrix}$$\n",
"\n",
"O sistema est\u00e1 escalonado e resolvemos com o algoritmo *BackStep* e a solu\u00e7\u00e3o \u00e9 $(0,1,-1)$.\n",
"\n"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"# A is type array from numpy\n",
"# remember: A[0.0] is the first element\n",
"# Elementary row operations:\n",
"\n",
"def TrocaLinha(A,i,j):\n",
" \"\"\"Troca as linhas i e j da matriz A\"\"\"\n",
" B=A.copy()\n",
" buffer = B[i].copy()\n",
" B[i] = B[j]\n",
" B[j] = buffer\n",
" return B\n",
"\n",
"\n",
"def SubsLinha(A,i,k,alfa, j=0):\n",
" '''soma a linha i com um multiplo da linha k a partir da coluna j'''\n",
" B=A.copy()\n",
" L,C=shape(B)\n",
" B[i,j:C] = B[i,j:C] + alfa*B[k,j:C]\n",
" return B\n",
"\n",
"def BackStep(A,b):\n",
" \"\"\" Se a matriz for triangular superior resolve Ax=b \n",
" A \u00e9 uma matriz triangular.\n",
" b \u00e9 uma lista de n\u00fameros reais de mesma dimens\u00e3o de A. \"\"\"\n",
" #i Vamos testar se A \u00e9 TS e se a a diagonal n\u00e3o tem zero\n",
" # Esta parte consome muito tempo e n\u00e3o deve ser usada depois\n",
" #n = len(A)\n",
" #for i in range(n):\n",
" # teste1 = A[i,i]==0\n",
" # if i==n:\n",
" # teste2 = False\n",
" # else: teste2 = False in (A[i+1:n,i]==0)\n",
" # if teste1 or teste2:\n",
" # raise ValueError(\" Matriz A n\u00e3o pode ser parametro do BackStep \")\n",
" # Esta foi s\u00f3 a parte de teste.\n",
" n=shape(A)[0]\n",
" x=np.zeros(n)\n",
" k=n-1\n",
" while k>=0:\n",
" x[k] = (1/A[k,k])*(b[k]-sum([A[k,j]*x[j] for j in range(k+1,n)]))\n",
" k-=1\n",
" return x\n",
"\n",
"# step k of gauss elimination:\n",
"def GaussStep(A,k):\n",
" '''O k-esimo passo na eliminacao de gauss'''\n",
" L=shape(A)[0] # numero de linhas da matriz A\n",
" B=A.copy()\n",
" for j in range(k+1,L):\n",
" m=-B[j,k]/B[k,k]\n",
" B=SubsLinha(B,j,k,m,k)\n",
" B[j,k]=-m # Aproveitamos a matriz A pra guardar o multiplicador(ta errado!_)\n",
" return B\n",
"\n",
"# triangle form, without care and pivoting\n",
"def TriangleForm(A):\n",
" '''Retorna a matriz escalonada,\n",
" com os multiplicadores na parte triangular inferior'''\n",
" B=A.copy()\n",
" L=shape(B)[0] # numero de linhas = numero de incognitas\n",
" for k in range(L-1):\n",
" p=pivo(B,k)\n",
" B=TrocaLinha(B,k,p)\n",
" B = GaussStep(B,k)\n",
" \n",
" return B\n",
" \n",
"def pivo(A,n):\n",
" ''' retorna o indice de pivotacao da matriz A na coluna n,\n",
" a partir da linha n'''\n",
" tamanho = shape(A)\n",
" l=tamanho[0] # numero de linhas de A\n",
" c= tamanho[1] # numero de colunas\n",
" if (l<=n) or (c<=n):\n",
" raise ValueError(\"n deve ser menor que dimensao de A\")\n",
" p=n # inicio o pivo como n, e vou aumentando\n",
" for k in range(n,l):\n",
" if abs(A[k,n]) > abs(A[p,n]):\n",
" p=k\n",
" return p"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 16
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A = np.array([[3.,2.,1.,1],[1,4,2,2],[6,1,0,1]])\n",
"A"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 17,
"text": [
"array([[ 3., 2., 1., 1.],\n",
" [ 1., 4., 2., 2.],\n",
" [ 6., 1., 0., 1.]])"
]
}
],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A1 = TrocaLinha(A,0,pivo(A,0))\n",
"A1\n",
" "
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 18,
"text": [
"array([[ 6., 1., 0., 1.],\n",
" [ 1., 4., 2., 2.],\n",
" [ 3., 2., 1., 1.]])"
]
}
],
"prompt_number": 18
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A2 = GaussStep(A1,0)\n",
"A2"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 19,
"text": [
"array([[ 6. , 1. , 0. , 1. ],\n",
" [ 0.16666667, 3.83333333, 2. , 1.83333333],\n",
" [ 0.5 , 1.5 , 1. , 0.5 ]])"
]
}
],
"prompt_number": 19
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A3 = GaussStep(A2,1)\n",
"A3"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 20,
"text": [
"array([[ 6. , 1. , 0. , 1. ],\n",
" [ 0.16666667, 3.83333333, 2. , 1.83333333],\n",
" [ 0.5 , 0.39130435, 0.2173913 , -0.2173913 ]])"
]
}
],
"prompt_number": 20
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A4= A3[:,0:3]\n",
"A4"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 21,
"text": [
"array([[ 6. , 1. , 0. ],\n",
" [ 0.16666667, 3.83333333, 2. ],\n",
" [ 0.5 , 0.39130435, 0.2173913 ]])"
]
}
],
"prompt_number": 21
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"B1= A3[:,3]\n",
"B1"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 22,
"text": [
"array([ 1. , 1.83333333, -0.2173913 ])"
]
}
],
"prompt_number": 22
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x=BackStep(A4,B1)\n",
"x"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 23,
"text": [
"array([ 9.25185854e-17, 1.00000000e+00, -1.00000000e+00])"
]
}
],
"prompt_number": 23
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Invers\u00e3o de matrizes\n",
"S = np.array([[3.,2,1,1,0,0],[1,4,2,0,1,0],[6,1,0,0,0,1]])\n",
"S"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 24,
"text": [
"array([[ 3., 2., 1., 1., 0., 0.],\n",
" [ 1., 4., 2., 0., 1., 0.],\n",
" [ 6., 1., 0., 0., 0., 1.]])"
]
}
],
"prompt_number": 24
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"S1 = TriangleForm(S)\n",
"print(S1)\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[[ 6. 1. 0. 0. 0. 1. ]\n",
" [ 0.16666667 3.83333333 2. 0. 1. -0.16666667]\n",
" [ 0.5 0.39130435 0.2173913 1. -0.39130435 -0.43478261]]\n"
]
}
],
"prompt_number": 26
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Ac= S1[:,0:3]\n",
"Ac"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 27,
"text": [
"array([[ 6. , 1. , 0. ],\n",
" [ 0.16666667, 3.83333333, 2. ],\n",
" [ 0.5 , 0.39130435, 0.2173913 ]])"
]
}
],
"prompt_number": 27
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"I1 = S1[:,3]\n",
"I2 = S1[:,4]\n",
"I3 = S1[:,5]\n",
"print(I2)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[ 0. 1. -0.39130435]\n"
]
}
],
"prompt_number": 29
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"X1 = BackStep(Ac,I1)\n",
"X2 = BackStep(Ac,I2)\n",
"X3 = BackStep(Ac,I3)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 30
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"X1\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 31,
"text": [
"array([ 0.4, -2.4, 4.6])"
]
}
],
"prompt_number": 31
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Ain = np.array([X1,X2,X3])\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 32
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Ain\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 33,
"text": [
"array([[ 4.00000000e-01, -2.40000000e+00, 4.60000000e+00],\n",
" [ -2.00000000e-01, 1.20000000e+00, -1.80000000e+00],\n",
" [ 3.70074342e-17, 1.00000000e+00, -2.00000000e+00]])"
]
}
],
"prompt_number": 33
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Ain.T"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 35,
"text": [
"array([[ 4.00000000e-01, -2.00000000e-01, 3.70074342e-17],\n",
" [ -2.40000000e+00, 1.20000000e+00, 1.00000000e+00],\n",
" [ 4.60000000e+00, -1.80000000e+00, -2.00000000e+00]])"
]
}
],
"prompt_number": 35
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 36,
"text": [
"array([[ 3., 2., 1., 1.],\n",
" [ 1., 4., 2., 2.],\n",
" [ 6., 1., 0., 1.]])"
]
}
],
"prompt_number": 36
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"A4"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 37,
"text": [
"array([[ 6. , 1. , 0. ],\n",
" [ 0.16666667, 3.83333333, 2. ],\n",
" [ 0.5 , 0.39130435, 0.2173913 ]])"
]
}
],
"prompt_number": 37
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"An = S[:,0:3]\n",
"An"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 38,
"text": [
"array([[ 3., 2., 1.],\n",
" [ 1., 4., 2.],\n",
" [ 6., 1., 0.]])"
]
}
],
"prompt_number": 38
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"An*(Ain.T)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 39,
"text": [
"array([[ 1.20000000e+00, -4.00000000e-01, 3.70074342e-17],\n",
" [ -2.40000000e+00, 4.80000000e+00, 2.00000000e+00],\n",
" [ 2.76000000e+01, -1.80000000e+00, -0.00000000e+00]])"
]
}
],
"prompt_number": 39
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"inv(An)\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 42,
"text": [
"array([[ 4.00000000e-01, -2.00000000e-01, 5.55111512e-17],\n",
" [ -2.40000000e+00, 1.20000000e+00, 1.00000000e+00],\n",
" [ 4.60000000e+00, -1.80000000e+00, -2.00000000e+00]])"
]
}
],
"prompt_number": 42
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"An*inv(An)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 43,
"text": [
"array([[ 1.20000000e+00, -4.00000000e-01, 5.55111512e-17],\n",
" [ -2.40000000e+00, 4.80000000e+00, 2.00000000e+00],\n",
" [ 2.76000000e+01, -1.80000000e+00, -0.00000000e+00]])"
]
}
],
"prompt_number": 43
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment