Lab3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MOwNiT Laboratorium 3\n",
    "## Zadanie 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Napisać kod (w ulubionym** języku) do eliminacji Gaussa z i bez pivotingu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zwykła eliminacja Gaussa bez pivotingu:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def gauss(A, b):\n",
    "    n = len(A)\n",
    "    A = np.copy(A)\n",
    "    b = np.copy(b)\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if i != j:\n",
    "                scale = A[j, i]/A[i, i]\n",
    "                A[j] -= A[i]*scale\n",
    "                b[j] -= b[i]*scale\n",
    "\n",
    "    for i in range(n):\n",
    "        b[i] /= A[i, i]\n",
    "        \n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jeżeli współczynniki znacząco różnią się wielkością to powyższa metoda może doprowadzić do rozpropagowania się niedokładności wynikającej z błędu zaokrąglenia. W takich przypadkach lepiej zastosować wersję algorytmu, w którym pivotem jest największa wartość bezwzględna z kolumn macierzy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_pivoting(A, b):\n",
    "    n = len(A)\n",
    "    A = np.copy(A)\n",
    "    b = np.copy(b)\n",
    "    \n",
    "    for i in range(n-1):\n",
    "        maxidx = abs(A[i:,i]).argmax() + i\n",
    "        if maxidx != i:\n",
    "            A[[i,maxidx]] = A[[maxidx, i]]\n",
    "            b[[i,maxidx]] = b[[maxidx, i]]\n",
    "            \n",
    "        for row in range(i+1, n):\n",
    "            scale = A[row,i]/A[i,i]\n",
    "            A[row, i:] = A[row, i:] - scale*A[i, i:]\n",
    "            b[row] = b[row] - scale*b[i]\n",
    "\n",
    "    for i in range(n-1, -1, -1):\n",
    "        b[i] = (b[i] - np.dot(A[i,i+1:], b[i+1:])) / A[i,i]\n",
    "        \n",
    "    return b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poprawność obu algorytmów można przetestować porównując ich wynik z wynikiem zwracanym przez funkcję biblioteczną."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(True, True)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.matrix([[0.0001, -5.0300, 5.8090, 7.8320],\n",
    "               [2.2660, 1.9950,  1.2120, 8.0080],\n",
    "               [8.8500, 5.6810,  4.5520, 1.3020],\n",
    "               [6.7750, -2.253,  2.9080, 3.9700]])\n",
    "\n",
    "b = np.matrix([9.5740, 7.2190, 5.7300, 6.2910]).transpose()\n",
    "\n",
    "x = np.linalg.solve(A, b)\n",
    "y = gauss(A, b)\n",
    "z = gauss_pivoting(A, b)\n",
    "\n",
    "np.allclose(x, y), np.allclose(x, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zadanie 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zaimplementować algorytm faktoryzacji LU macierzy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mult_matrix(M, N):\n",
    "    return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in np.transpose(N)] for row_m in M]\n",
    "   \n",
    "\n",
    "def pivot_matrix(M):\n",
    "    n = len(M)\n",
    "    identity_matrix = np.identity(n)\n",
    "\n",
    "    for j in range(n):\n",
    "        row = max(range(j, n), key=(lambda x: abs(M[x][j])))\n",
    "        if j != row:\n",
    "            identity_matrix[j], identity_matrix[row] = identity_matrix[row], identity_matrix[j]\n",
    "    return identity_matrix\n",
    "\n",
    "\n",
    "def agh_superfast_lu(A):\n",
    "    n = len(A)\n",
    "                                                                                                                                                                                                               \n",
    "    L = np.zeros((n, n))\n",
    "    U = np.zeros((n, n))\n",
    "                                                                                                                                                                                            \n",
    "    P = pivot_matrix(A)\n",
    "    PA = mult_matrix(P, A)\n",
    "    \n",
    "    for i in range(n):\n",
    "        L[i][i] = 1.0\n",
    "        for j in range(i+1):\n",
    "            s = sum(U[k][i] * L[j][k] for k in range(j))\n",
    "            U[j][i] = PA[j][i] - s\n",
    "        for j in range(i, n):\n",
    "            s = sum(U[k][i] * L[j][k] for k in range(i))\n",
    "            L[j][i] = (PA[j][i] - s) / U[i][i]\n",
    "\n",
    "    return (P, L, U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zadanie 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zaimplementować funkcję sprawdzającą, czy dana macierz jest symetryczna i dodatnio określona."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def agh_superfast_check_spd(A):\n",
    "    return np.all(np.linalg.eigvals(A) > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zadanie 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zaimplementować algorytm faktoryzacji Cholesky'ego macierzy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def agh_superfast_cholesky(A):\n",
    "    n = len(A)\n",
    "    L = np.zeros((n, n))\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(i+1):\n",
    "            tmp = sum(L[i][k] * L[j][k] for k in range(j))\n",
    "            if (i == j):\n",
    "                L[i][j] = np.sqrt(A[i][i] - sum)\n",
    "            else:\n",
    "                L[i][j] = (1.0 / L[j][j] * (A[i][j] - tmp))\n",
    "    return L"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mownit",
   "language": "python",
   "name": ".venv"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}