sistemas_lineares(triang_LU_Gauss).ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/EduardoWS/Calculo-Numerico/blob/main/sistemas_lineares(triang_LU_Gauss).ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Triangular Superior\n",
        "\n",
        "def TriangularSup(A,b):\n",
        "  n = len(b)\n",
        "  x = [0]*n\n",
        "  x[n-1] = b[n-1]/A[n-1][n-1]\n",
        "\n",
        "  for i in list(range(n-1,0,-1)):\n",
        "    s = 0\n",
        "    for j in list(range(i+1,n+1,1)):\n",
        "      s = s + A[i-1][j-1]*x[j-1]\n",
        "\n",
        "    x[i-1] = (b[i-1]-s)/(A[i-1][i-1]) \n",
        "\n",
        "  return x\n",
        "  \n",
        "# Dados\n",
        "A = [[3,2,4],[0,(1/3),(2/3)],[0,0,1]]\n",
        "b = [1,(5/3),0]\n",
        "x = TriangularSup(A,b)\n",
        "\n",
        "print(x)\n",
        "\n",
        "################################################################################\n",
        "\n",
        "# Triangular Inferior\n",
        "\n",
        "def TriangularInf(A,b):\n",
        "  n = len(b)\n",
        "  x = [0]*n\n",
        "  x[0] = b[0]/A[0][0]\n",
        "\n",
        "  for i in list(range(1,n,1)):\n",
        "    s = 0\n",
        "    for j in list(range(0,i,1)):\n",
        "      s = s + A[i][j]*x[j]\n",
        "\n",
        "    x[i] = (b[i]-s)/(A[i][i]) \n",
        "\n",
        "  return x\n",
        "  \n",
        "# Dados\n",
        "import numpy as np\n",
        "A = np.array([[3,2,4],[0,(1/3),(2/3)],[0,0,1]])\n",
        "A = A.T\n",
        "b = [1,(5/3),0]\n",
        "x = TriangularInf(A,b)\n",
        "\n",
        "print(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "hmTnzFZSoyKT",
        "outputId": "637ee5db-d329-41fb-b720-e4df68d9dfb1"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[-3.0000000000000004, 5.000000000000001, 0.0]\n",
            "[0.3333333333333333, 3.0, -3.333333333333333]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#Decomposição LU com SymPy\n",
        "\n",
        "import sympy\n",
        "A = sympy.Matrix([[1, 2,0], [1, 3,1],[-2, 0,1] ])\n",
        "b = sympy.Matrix([3, 5,-1])\n",
        "L, U, _ = A.LUdecomposition()\n",
        "print(L)\n",
        "print(U)\n",
        "x = A.solve(b); \n",
        "print(x)\n",
        "\n",
        "#Decomposição LU\n",
        "import numpy as np\n",
        "\n",
        "def decompLU(A):\n",
        "     n = np.shape(A)[0]\n",
        "     L = np.eye(n) \n",
        "     U = np.zeros((n,n))\n",
        "     for k  in np.arange(n):\n",
        "        for j in np.arange(k,n):\n",
        "           U[k,j]=A[k,j]\n",
        "           for s in np.arange(k):\n",
        "               U[k,j] = U[k,j] - L[k,s]*U[s,j]\n",
        "        for i in np.arange(k+1,n):\n",
        "           L[i,k]=A[i,k]\n",
        "           for s in np.arange(k):\n",
        "               L[i,k] = L[i,k] - L[i,s]*U[s,k]\n",
        "           L[i,k] = L[i,k]/U[k,k]\n",
        "     return L, U   \n",
        "\n",
        "def LUsolve(A,b):\n",
        "   n = np.shape(A)[0]  \n",
        "   L,U = decompLU(A)\n",
        "   for k in range(1,n):\n",
        "      b[k] = b[k] - L[k,0:k].dot(b[0:k])\n",
        "   for k in range(n-1,-1,-1):\n",
        "      b[k] = (b[k] - U[k,k+1:n].dot(b[k+1:n]))/U[k,k]\n",
        "   return b\n",
        "\n",
        "def LUsolve2(A,b):\n",
        "   n = np.shape(A)[0]  \n",
        "   L,U = decompLU(A)\n",
        "   y = TriangularInf(L,b)\n",
        "   x = TriangularSup(U,y)\n",
        "   return x\n",
        "\n",
        "#Exemplo \n",
        "A = np.array([[1, 2,0], [1, 3,1],[-2, 0,1]])\n",
        "b = np.array([3,5,-1])\n",
        "x = LUsolve2(A,b)\n",
        "print(L)\n",
        "print(U)\n",
        "print(x)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yslFLMNiofrK",
        "outputId": "c9c28631-6411-46ae-f6db-4eb910988cbd"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Matrix([[1, 0, 0], [1, 1, 0], [-2, 4, 1]])\n",
            "Matrix([[1, 2, 0], [0, 1, 1], [0, 0, -3]])\n",
            "Matrix([[1], [1], [1]])\n",
            "Matrix([[1, 0, 0], [1, 1, 0], [-2, 4, 1]])\n",
            "Matrix([[1, 2, 0], [0, 1, 1], [0, 0, -3]])\n",
            "[1.0, 1.0, 1.0]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WtZxUW5-oeNO",
        "outputId": "d187b717-166d-4aa0-8400-239850afe674"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[0.0, 1.0, -1.0]\n"
          ]
        }
      ],
      "source": [
        "#Eliminação de Gauss\n",
        "\n",
        "import numpy as np\n",
        "\n",
        "def ElimGauss(a,b):\n",
        " n = len(b)\n",
        " for k in range(0,n-1):\n",
        "   for i in range(k+1,n):\n",
        "     if a[i,k] != 0.0:\n",
        "        lam = a[i,k]/a[k,k]\n",
        "        a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]\n",
        "        b[i] = b[i] - lam*b[k]        \n",
        " return a, b\n",
        "\n",
        "def SolveSis(a,b):\n",
        " n = len(b)\n",
        " a, b = ElimGauss(a,b)\n",
        " x = TriangularSup(a,b)\n",
        " return x\n",
        " #for k in range(n-1,-1,-1):\n",
        " #   b[k] = (b[k] - a[k,k+1:n].dot(b[k+1:n]))/a[k,k]\n",
        " #return b\n",
        "\n",
        "#Exemplo \n",
        "A = np.array([[2, 1,2], [4, 3,3],[6, 5,-1]])\n",
        "b = np.array([-1,0,6])\n",
        "x = SolveSis(A,b)\n",
        "print(x)"
      ]
    }
  ]
}