sistemas_lineares(Gauss_LU_Pivo_LU_banda_Cholesky).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(Gauss_LU_Pivo_LU_banda_Cholesky).ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2JwuuB4ypy-g",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "cb39b6ac-4e1e-4011-8cae-567baf05e457"
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "def sub_regressiva(U,y):\n",
        "    n = np.shape(y)[0]\n",
        "    x = np.zeros(n)\n",
        "\n",
        "    for i in np.arange(n-1,-1,-1):\n",
        "        x[i] = (y[i]-U[i,i+1:n].dot(x[i+1:n]))/U[i,i]\n",
        "\n",
        "    return x\n",
        "\n",
        "def sub_progressiva(L,b):\n",
        "    n = np.shape(b)[0]\n",
        "    z = np.zeros(n)\n",
        "\n",
        "    for i in np.arange(n):\n",
        "        z[i] = (b[i]-L[i,0:i].dot(z[0:i]))/L[i,i]\n",
        "\n",
        "    return z\n",
        "\n",
        "# Gauss com pivô\n",
        "def eliminacao_gauss_pivo(A,b):\n",
        "    n = np.shape(A)[0]\n",
        "    newA = np.copy(A)\n",
        "    newb = np.copy(b)\n",
        "\n",
        "    for k in np.arange(n-1):\n",
        "        p = np.argmax(np.abs(newA[k:n,k]))\n",
        "        newA[[k,p],k:n] = newA[[p,k],k:n]\n",
        "        newb[[k,p]] = newb[[p,k]]\n",
        "\n",
        "        for i in np.arange(k+1,n):\n",
        "            m = -newA[i,k]/newA[k,k]\n",
        "            newA[i,k:n] = newA[i,k:n] + m*newA[k,k:n]\n",
        "            newb[i] = newb[i] + m*newb[k]\n",
        "\n",
        "    return newA,newb\n",
        "\n",
        "# LU com pivô\n",
        "def fatoraLUpivo(A):\n",
        "    U = np.copy(A)  \n",
        "    n = np.shape(U)[0]  \n",
        "    L = np.eye(n)  \n",
        "    P = np.copy(L)\n",
        "\n",
        "    for j in np.arange(n-1):\n",
        "        p = np.argmax(np.abs(U[j:n,j]))\n",
        "        P[[j,p],:] = P[[p,j],:]\n",
        "        U[[j,p],j:n] = U[[p,j],j:n]\n",
        "        L[[j,p],0:j] = L[[p,j],0:j]\n",
        "\n",
        "        for i in np.arange(j+1,n):\n",
        "            L[i,j] = U[i,j]/U[j,j] \n",
        "            U[i,j+1:n] = U[i,j+1:n] - L[i,j]*U[j,j+1:n]\n",
        "            U[i,j] = 0 \n",
        "    return L,U,P\n",
        "\n",
        "# Exemplo\n",
        "A = np.array([[1, 1, 1],  \n",
        "              [1, 0,10],  \n",
        "              [0,10, 1]], dtype='double')\n",
        "b = [0,-48,25]\n",
        "\n",
        "(newA,newb) = eliminacao_gauss_pivo(A,b)\n",
        "x_pivo_gauss = sub_regressiva(newA,newb)\n",
        "print(x_pivo_gauss)\n",
        "\n",
        "(L,U,P) = fatoraLUpivo(A)\n",
        "bpivo = P.dot(b)\n",
        "ypivo = sub_progressiva(L,bpivo)\n",
        "x_pivo_lu = sub_regressiva(U,ypivo)\n",
        "print(x_pivo_lu)\n"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[ 2.  3. -5.]\n",
            "[ 2.  3. -5.]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#LU banda\n",
        "\n",
        "import numpy as np\n",
        "import time\n",
        "\n",
        "def LU( A ):\n",
        "    \n",
        "    n = A.shape[ 0 ]\n",
        "    \n",
        "    U = A.copy()\n",
        "    L = np.eye( n )\n",
        "    \n",
        "    for j in range( n - 1 ):\n",
        "        \n",
        "        for i in range( j + 1, n ):\n",
        "                \n",
        "            L[ i, j ] = U[ i, j ] / U[ j, j ]\n",
        "            U[ i, j : n ] = U[ i, j : n ] - L[ i, j ] * U[ j, j : n ]\n",
        "            \n",
        "    return ( L, U )\n",
        "\n",
        "\n",
        "def LU_banda( A, p ):\n",
        "    \n",
        "    n = A.shape[ 0 ]\n",
        "    \n",
        "    U = A.copy()\n",
        "    L = np.eye( n )\n",
        "    \n",
        "    for j in range( n - 1 ):\n",
        "        \n",
        "        v = min( n, j + p + 1 )\n",
        "        for i in range( j + 1, v ):\n",
        "                \n",
        "            L[ i, j ] = U[ i, j ] / U[ j, j ]\n",
        "            U[ i, j : v ] = U[ i, j : v ] - L[ i, j ] * U[ j, j : v ]\n",
        "            \n",
        "    return ( L, U )\n",
        "\n",
        "\n",
        "n = 2000\n",
        "p = 2\n",
        "\n",
        "A = np.zeros( ( n, n ) )\n",
        "for i in range( n ):\n",
        "    for j in range( max( 0, i - p ), min( n, i + p + 1 ) ):\n",
        "        A[ i, j ] = np.random.normal()    \n",
        "\n",
        "\n",
        "start_time = time.time()\n",
        "( L, U ) = LU( A )\n",
        "end_time = time.time()\n",
        "print( end_time - start_time )\n",
        "\n",
        "\n",
        "start_time = time.time()\n",
        "( L_banda, U_banda ) = LU_banda( A, p )\n",
        "end_time = time.time()\n",
        "print( end_time - start_time )\n",
        "\n",
        "\n",
        "print( np.linalg.norm( L @ U - A ) )\n",
        "print( np.linalg.norm( L_banda @ U_banda - A ) )"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6eKuFJ6MsYDO",
        "outputId": "0e7b6ee4-5edc-4adf-e404-865a966ca221"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "35.78927946090698\n",
            "0.039710044860839844\n",
            "2.594287423993927e-13\n",
            "2.594287423993927e-13\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from math import sqrt\n",
        "import numpy as np\n",
        "\n",
        "def cholesky(A):\n",
        "    n = np.shape(A)[0]\n",
        "\n",
        "    L = np.zeros((n,n))\n",
        "\n",
        "    for i in np.arange(n):\n",
        "        for j in np.arange(i+1):\n",
        "            tmp_sum = sum(L[i][k]*L[j][k] \n",
        "            for k in np.arange(j))           \n",
        "            if (i == j):\n",
        "                L[i][j] = sqrt(A[i][i] - tmp_sum)\n",
        "            else:\n",
        "                L[i][j] = ((A[i][j] - tmp_sum)/L[j][j])\n",
        "    return L\n",
        " \n",
        "A = [[3,2,0], [2,5,1], [0,1,3]]\n",
        "L = cholesky(A)\n",
        "\n",
        "print(L)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Bn2JF_jgyVMm",
        "outputId": "60661865-fd5e-4ba4-9789-996767294061"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[1.73205081 0.         0.        ]\n",
            " [1.15470054 1.91485422 0.        ]\n",
            " [0.         0.52223297 1.65144565]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "rhsogcpRvcoT"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}