SebastianKusnierzLab10.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "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.8.2"
    },
    "colab": {
      "name": "lab10.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VyOEcJKEVDdz",
        "colab_type": "text"
      },
      "source": [
        "##  Rozwiązywanie układu równań liniowych metodami rozkładu LU \n",
        "\n",
        "Podaj zasadę działania metod opartych o dekompozycję LU. \n",
        "\n",
        "#### Zadanie 1\n",
        "Zaimplementuj metody rozkładu LU:\n",
        "- Rozkład Crouta \n",
        "- Rozkład Doolitla\n",
        "- Rozkład Choleskyego \n",
        "\n",
        "Dla każdej z metod podaj warunki niezbędne aby można ją było zastosować. Sprawdź poprawność działania tych metod. \n",
        "\n",
        "Przetestuj wydajność algorytmów dla kilku rozmiarów macierzy (podobnie jak w ćwiczeniu 9). \n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pgHpSY6fj7UI",
        "colab_type": "text"
      },
      "source": [
        "Dekompozycja LU - jest to jedna z metod rozwiązywania układu równań nieliniowych, nazwa ta pochodzi od użytych w tej metodzie macierzy trójkątnych tj. dolnotrójkątnej(L) oraz górnotrójkątnej(U).   \n",
        "Niech dany będzie układ równań liniowych:    \n",
        "$A * x = y$      \n",
        ",gdzie $A$ – macierz współczynników, $x$ – wektor niewiadomych, $y$ – wektor danych(wyrazów wolnych).\n",
        "\n",
        "W metodzie LU macierz współczynników $A$ zapisywana jest jako iloczyn pewnych macierzy dolnej $L$ i górnej $U$\n",
        "\n",
        "$A = L*U$   \n",
        ", gdzie :\n",
        "![alt text](https://wikimedia.org/api/rest_v1/media/math/render/svg/0365e58fcda12c0e33a5509d37fd0658cefba13d) \n",
        "\n",
        "Układ równań przyjmuje wówczas postać :\n",
        "\n",
        "$L*U*x = y$         \n",
        "a jego rozwiązanie sprowadza się do rozwiązania dwóch układów równań z macierzami trójkątnymi:\n",
        "\n",
        "$L * z = y$   \n",
        "$U * x = z$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H53HWeW-mRQ7",
        "colab_type": "text"
      },
      "source": [
        "- **Rozkład Crouta** - aby zastosować rozkład Crouta musi zostać spełniony następujący warunek, na diagonalnej macierzy trójkątnej dolnej $U$ muszą występować wyłacznie 1 ($u_{ii} = 1$).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FZLCaompjdGq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import math\n",
        "import time\n",
        "import numpy as np\n",
        "\n",
        "\n",
        "def crout(A):\n",
        "    n = len(A)\n",
        "    L = [[0] * n for i in range(n)]\n",
        "    U = [[0] * n for i in range(n)]\n",
        "\n",
        "    for j in range(n):\n",
        "        U[j][j] = 1\n",
        "        for i in range(j, n):\n",
        "            alpha = float(A[i][j])\n",
        "            for k in range(j):\n",
        "                alpha -= L[i][k] * U[k][j]\n",
        "            L[i][j] = alpha\n",
        "        for i in range(j + 1, n):\n",
        "            tmpU = float(A[j][i])\n",
        "            for k in range(j):\n",
        "                tmpU -= L[j][k] * U[k][i]\n",
        "            if int(L[j][j]) == 0:\n",
        "                L[j][j] = 0\n",
        "            U[j][i] = tmpU / L[j][j]\n",
        "    return [L, U]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vbtTjd5toQb-",
        "colab_type": "text"
      },
      "source": [
        "- **Rozkład Doolitla** - aby zastosować rozkład Dollitla musi zostać spełniony następujący warunek, na diagonalnej macierzy trójkątnej górnej $L$ muszą występować wyłacznie 1 ($l_{ii} = 1$)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iE4FvAzkonbR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def dollitle(A):\n",
        "    n = len(A)\n",
        "    L = [[0] * n for i in range(n)]\n",
        "    U = [[0] * n for i in range(n)]\n",
        "\n",
        "    for j in range(n):\n",
        "        L[j][j] = 1.0\n",
        "        for i in range(j + 1):\n",
        "            s1 = sum(U[k][j] * L[i][k] for k in range(i))\n",
        "            U[i][j] = A[i][j] - s1\n",
        "        for i in range(j, n):\n",
        "            s2 = sum(U[k][j] * L[i][k] for k in range(j))\n",
        "            L[i][j] = (A[i][j] - s2) / U[j][j]\n",
        "\n",
        "    return [L, U]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZZd8-pX0ooas",
        "colab_type": "text"
      },
      "source": [
        "- **Rozkład Choleskyego** - aby zastosować rozkład Choleskyego muszą zostać spełniony następujący waruneki,   \n",
        " - wartości na diagonalnej macierzy trójkątnej górnej $L$ oraz diagonalnej   macierz trójkątnej dolnej $U$  muszą być sobie równe ($u_{ii} = l_{ii}$),\n",
        " - macierz $A$ musi być symetryczna oraz dodatnie określona."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jzDp6dmnpgJT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def cholesky(A):\n",
        "    n = len(A)\n",
        "    L = [[0.0] * n for i in range(n)]\n",
        "\n",
        "    for i in range(n):\n",
        "        for k in range(i + 1):\n",
        "            tmp_sum = sum(L[i][j] * L[k][j] for j in range(k))\n",
        "            if (i == k):\n",
        "                L[i][k] = math.sqrt(A[i][i] - tmp_sum)\n",
        "            else:\n",
        "                L[i][k] = (1.0 / L[k][k] * (A[i][k] - tmp_sum))\n",
        "    return [L, np.array(L).transpose()]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9Pd0TxzwpvJU",
        "colab_type": "text"
      },
      "source": [
        "Poprawności działania tych metod zbadam na macierzy 3x3 oraz stworzyłem dodatkową metode create_rand_symmetric_and_positive_matrix, aby zapewnić spełnienie warunków dla metody cholesky "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "To-YH-hUp9GA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def solve_equtition(A, B, fun):\n",
        "    [L, U] = fun(A)\n",
        "    z = np.linalg.solve(L, B)\n",
        "    x = np.linalg.solve(U, z)\n",
        "    return x\n",
        "\n",
        "\n",
        "def create_rand_matrix(min, max, cols, rows):\n",
        "    return np.random.randint(low=min, high=max + 1, size=(rows, cols))\n",
        "\n",
        "\n",
        "def create_rand_vector(min, max, size):\n",
        "    return np.random.randint(low=min, high=max + 1, size=(size,))\n",
        "\n",
        "\n",
        "def create_rand_symmetric_and_positive_matrix(min, max, cols, rows):\n",
        "    A = create_rand_matrix(min, max, cols, rows)\n",
        "    B = np.dot(A, A.transpose())\n",
        "    return (B + B.T) / 2\n",
        "\n",
        "\n",
        "def print_matrix(A):\n",
        "    for i in A:\n",
        "        print(i)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P7xlltUhpufl",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 493
        },
        "outputId": "0ee7241a-548f-454b-c116-244ba5c07dff"
      },
      "source": [
        "A = create_rand_symmetric_and_positive_matrix(3, 10, 3, 3)\n",
        "Y = create_rand_vector(3, 10, 3)\n",
        "[L_crout, U_crout] = crout(A)\n",
        "[L_doll, U_doll] = dollitle(A)\n",
        "[L_chol, U_chol] = cholesky(A)\n",
        "\n",
        "print(\"Macierz A ma postać :\")\n",
        "print_matrix(A)\n",
        "print(\"Macierz L dla rozkładu Crouta:\")\n",
        "print_matrix(L_crout)\n",
        "print(\"Macierz U dla rozkładu Crouta:\")\n",
        "print_matrix(U_crout)\n",
        "print(\"Macierz L dla rozkładu Dollitla:\")\n",
        "print_matrix(L_doll)\n",
        "print(\"Macierz U dla rozkładu Dollitla:\")\n",
        "print_matrix(U_doll)\n",
        "print(\"Macierz L dla rozkładu Choleskyego:\")\n",
        "print_matrix(L_chol)\n",
        "print(\"Macierz U dla rozkładu Choleskyego:\")\n",
        "print_matrix(U_chol)"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Macierz A ma postać :\n",
            "[ 74. 123.  48.]\n",
            "[123. 245.  99.]\n",
            "[48. 99. 43.]\n",
            "Macierz L dla rozkładu Crouta:\n",
            "[74.0, 0, 0]\n",
            "[123.0, 40.55405405405406, 0]\n",
            "[48.0, 19.216216216216225, 2.7594135288237247]\n",
            "Macierz U dla rozkładu Crouta:\n",
            "[1, 1.662162162162162, 0.6486486486486487]\n",
            "[0, 1, 0.4738420526491167]\n",
            "[0, 0, 1]\n",
            "Macierz L dla rozkładu Dollitla:\n",
            "[1.0, 0, 0]\n",
            "[1.662162162162162, 1.0, 0]\n",
            "[0.6486486486486487, 0.4738420526491167, 1.0]\n",
            "Macierz U dla rozkładu Dollitla:\n",
            "[74.0, 123.0, 48.0]\n",
            "[0, 40.55405405405406, 19.216216216216225]\n",
            "[0, 0, 2.7594135288237283]\n",
            "Macierz L dla rozkładu Choleskyego:\n",
            "[8.602325267042627, 0.0, 0.0]\n",
            "[14.298459565489772, 6.368206502152238, 0.0]\n",
            "[5.579886659703326, 3.01752404067327, 1.6611482561239743]\n",
            "Macierz U dla rozkładu Choleskyego:\n",
            "[ 8.60232527 14.29845957  5.57988666]\n",
            "[0.         6.3682065  3.01752404]\n",
            "[0.         0.         1.66114826]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jCFbreMSrb73",
        "colab_type": "text"
      },
      "source": [
        "Analizując powyższe wyniki widzimy, że każdy warunek dla danej metody został spełniony.         \n",
        "Następnie przeanalizuje czy wszystkie powyższe metody dają poprawnę rozwiązanie : "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b4kp_ylew85v",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 119
        },
        "outputId": "116fa17b-2e09-432d-e7e3-a9ec466f2558"
      },
      "source": [
        "def solve_equtition(A, Y, fun):\n",
        "    [L, U] = fun(A)\n",
        "    z = np.linalg.solve(L, Y)\n",
        "    x = np.linalg.solve(U, z)\n",
        "    return x\n",
        "\n",
        "print(\"Metoda Crouta\")\n",
        "print(solve_equtition(A, Y, crout))\n",
        "print(\"Metoda Dollitla\")\n",
        "print(solve_equtition(A, Y, dollitle))\n",
        "print(\"Metoda Choleskyego\")\n",
        "print(solve_equtition(A, Y, cholesky))"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Metoda Crouta\n",
            "[ 0.66731071 -1.17485811  2.16930322]\n",
            "Metoda Dollitla\n",
            "[ 0.66731071 -1.17485811  2.16930322]\n",
            "Metoda Choleskyego\n",
            "[ 0.66731071 -1.17485811  2.16930322]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N-r0F6qUxixm",
        "colab_type": "text"
      },
      "source": [
        "Widzimy również, że wszystkie 3 metody dają identyczne poprawne rozwiązanie."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tYxaL1tAwmLd",
        "colab_type": "text"
      },
      "source": [
        "Do testowania wydajności powyższych funkcji zastowałem poniższą metodę : "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E2jVms7HwlSH",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compare_exec_time(min, max, N):\n",
        "    A = create_rand_symmetric_and_positive_matrix(min, max, N, N)\n",
        "    B = create_rand_vector(min, max, N)\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, crout)\n",
        "    print(\"---LU Crouta took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, dollitle)\n",
        "    print(\"---LU Doolitla took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, cholesky)\n",
        "    print(\"---LU Choleskyego took %s seconds ---\" % (time.time() - start_time))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y7gW0_qhx6wh",
        "colab_type": "text"
      },
      "source": [
        "Wyniki testów wykonałem poza środowiskiem jupytera, a same wyniki są nastepujące :     \n",
        "- Dla macierzy 100x100 zakres losowanych liczb [0, 500] :    \n",
        "  ---LU Crouta took 0.059658050537109375 seconds ---     \n",
        "  ---LU Doolitla took 0.1942441463470459 seconds ---     \n",
        "  ---LU Choleskyego took 0.14393281936645508 seconds ---    \n",
        "- Dla macierzy 200x200 zakres losowanych liczb [0, 500] :    \n",
        "  ---LU Crouta took 0.29308199882507324 seconds ---     \n",
        "  ---LU Doolitla took 0.7315621376037598 seconds ---     \n",
        "  ---LU Choleskyego took 0.40914058685302734 seconds ---   \n",
        "- Dla macierzy 300x300 zakres losowanych liczb [0, 500] :     \n",
        "  ---LU Crouta took 0.8489158153533936 seconds ---    \n",
        "  ---LU Doolitla took 2.247528553009033 seconds ---     \n",
        "  ---LU Choleskyego took 1.131791353225708 seconds ---  \n",
        "- Dla macierzy 500x500 zakres losowanych liczb [0, 500] :      \n",
        "  ---LU Crouta took 3.9642786979675293 seconds ---     \n",
        "  ---LU Doolitla took 9.914392709732056 seconds ---      \n",
        "  ---LU Choleskyego took 4.811501741409302 seconds ---     \n",
        "- Dla macierzy 1000x1000 zakres losowanych liczb [0, 500] :     \n",
        "  ---LU Crouta took 37.514655351638794 seconds ---    \n",
        "  ---LU Doolitla took 88.99602842330933 seconds ---      \n",
        "  ---LU Choleskyego took 38.43804574012756 seconds ---     "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yBfIvp1vulIZ",
        "colab_type": "text"
      },
      "source": [
        "#### Zadanie 2 (dla chętnych)\n",
        "\n",
        "#### Zadanie 3 \n",
        "Zapoznaj się z funkcją rozwiązywania układów równań liniowych dostarczoną przez bibliotekę numpy/scipy. Porównaj jej wydajność z własnymi implementacjami.  \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UP3lddkq10j3",
        "colab_type": "text"
      },
      "source": [
        "Biblioteki numpy oraz scipy dostarczają nam funkcję $linalg.solve()$ służącą do rozwiązywania równań.  \n",
        "Do porównania wydajności tych funkcji z moją implementacją użyłem następującej funkncji :    "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E3IT1ToqVDd0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def compare_exec_time_2(min, max, N):\n",
        "    A = create_rand_symmetric_and_positive_matrix(min, max, N, N)\n",
        "    B = create_rand_vector(min, max, N)\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, crout)\n",
        "    print(\"---LU Crouta took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, dollitle)\n",
        "    print(\"---LU Doolitla took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    solve_equtition(A, B, cholesky)\n",
        "    print(\"---LU Choleskyego took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    np.linalg.solve(A, B)\n",
        "    print(\"---Numpy solver took %s seconds ---\" % (time.time() - start_time))\n",
        "    start_time = time.time()\n",
        "    scipy.linalg.solve(A, B)\n",
        "    print(\"---Scipy solver took %s seconds ---\" % (time.time() - start_time))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PDZDtsx62cIy",
        "colab_type": "text"
      },
      "source": [
        "Wyniki testów wykonałem poza środowiskiem jupytera, a same wyniki są nastepujące :\n",
        "\n",
        "- Dla macierzy 100x100 zakres losowanych liczb [0, 500] :    \n",
        "  ---LU Crouta took 0.05288290977478027 seconds ---           \n",
        "  ---LU Doolitla took 0.19524312019348145 seconds ---           \n",
        "  ---LU Choleskyego took 0.1520702838897705 seconds ---         \n",
        "  ---Numpy solver took 0.0009083747863769531 seconds ---        \n",
        "  ---Scipy solver took 0.0009958744049072266 seconds ---        \n",
        "- Dla macierzy 200x200 zakres losowanych liczb [0, 500] :         \n",
        "  ---LU Crouta took 0.2916877269744873 seconds ---          \n",
        "  ---LU Doolitla took 0.7662258148193359 seconds ---           \n",
        "  ---LU Choleskyego took 0.41582250595092773 seconds ---       \n",
        "  ---Numpy solver took 0.0016617774963378906 seconds ---          \n",
        "  ---Scipy solver took 0.0023343563079833984 seconds ---       \n",
        "- Dla macierzy 300x300 zakres losowanych liczb [0, 500] :              \n",
        "  ---LU Crouta took 0.9223763942718506 seconds ---             \n",
        "  ---LU Doolitla took 2.279259443283081 seconds ---           \n",
        "  ---LU Choleskyego took 1.1543560028076172 seconds ---              \n",
        "  ---Numpy solver took 0.0030765533447265625 seconds ---            \n",
        "  ---Scipy solver took 0.004083871841430664 seconds ---           \n",
        "- Dla macierzy 500x500 zakres losowanych liczb [0, 500] :         \n",
        "  ---LU Crouta took 4.130779504776001 seconds ---           \n",
        "  ---LU Doolitla took 10.669626712799072 seconds ---           \n",
        "  ---LU Choleskyego took 4.949918985366821 seconds ---                    \n",
        "  ---Numpy solver took 0.006085872650146484 seconds ---       \n",
        "  ---Scipy solver took 0.009285449981689453 seconds ---           \n",
        "- Dla macierzy 1000x1000 zakres losowanych liczb [0, 500] :       \n",
        "  ---LU Crouta took 36.34651827812195 seconds ---        \n",
        "  ---LU Doolitla took 86.19828629493713 seconds ---                   \n",
        "  ---LU Choleskyego took 40.54396080970764 seconds ---             \n",
        "  ---Numpy solver took 0.016066551208496094 seconds ---              \n",
        "  ---Scipy solver took 0.02889227867126465 seconds ---                    \n"
      ]
    }
  ]
}