{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD4: Systèmes linéaires : méthodes directes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Programmer une fonction qui prend en entrée une matrice $A$ triangulaire inférieure (resp.\n",
    "triangulaire supérieure) et un vecteur $b$, qui résout le système $Ax = b$ par une méthode de\n",
    "descente (resp. remontée) et qui donne le vecteur solution $x$ en sortie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trisup(U, b):\n",
    "    \"\"\"Résout un système triangulaire supérieur Ux = b.\n",
    "\n",
    "    Args:\n",
    "        U (np.ndarray): Matrice U.\n",
    "        b (np.ndarray): Vecteur b.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Vecteur x solution de Ux = b.\n",
    "    \"\"\"\n",
    "    n = U.shape[0]\n",
    "    x = np.zeros(n)\n",
    "    for i in reversed(range(n)):\n",
    "        x[i] = (b[i] - sum(x[j]*U[i,j] for j in range(i+1, n)))/U[i,i]\n",
    "    return x\n",
    "\n",
    "def triinf(L, b):\n",
    "    \"\"\"Résout un système triangulaire inférieur Lx = b.\n",
    "\n",
    "    Args:\n",
    "        L (np.ndarray): Matrice L.\n",
    "        b (np.ndarray): Vecteur b.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Vecteur x solution de Lx = b.\n",
    "    \"\"\"\n",
    "    n = L.shape[0]\n",
    "    x = np.zeros(n)\n",
    "    for i in range(n):\n",
    "        x[i] = (b[i] - sum(x[j]*L[i,j] for j in range(i+1)))/L[i,i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Tester vos fonctions sur les systèmes suivants :\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    1 & 0 & 0 \\\\\n",
    "    2 & 3 & 0 \\\\\n",
    "    1 & 4 & -1\n",
    "\\end{bmatrix} x = \n",
    "\\begin{bmatrix}\n",
    "    1 \\\\\n",
    "    8 \\\\\n",
    "    10\n",
    "\\end{bmatrix} \\qquad \n",
    "\\begin{bmatrix}\n",
    "    1 & 2 & 3 \\\\\n",
    "    0 & 4 & 8 \\\\\n",
    "    0 & 0 & 5\n",
    "\\end{bmatrix} x = \n",
    "\\begin{bmatrix}\n",
    "    6 \\\\\n",
    "    16 \\\\\n",
    "    15\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Vérifiez vos résultats avec la fonction `np.linalg.solve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = np.array([\n",
    "    [1, 2, 3],\n",
    "    [0, 4, 8],\n",
    "    [0, 0, 5]\n",
    "])\n",
    "\n",
    "b = np.array([6, 16, 15])\n",
    "\n",
    "print(f\"Résultats de l'algorithme trisup {trisup(U, b)}\")\n",
    "print(f\"Résultats avec solve {np.linalg.solve(U, b)}\")\n",
    "\n",
    "L = np.array([\n",
    "    [1, 0, 0],\n",
    "    [2, 3, 0],\n",
    "    [1, 4, -1]\n",
    "])\n",
    "\n",
    "b = np.array([1, 8, 10])\n",
    "\n",
    "print(f\"Résultats de l'algorithme triinf {triinf(L, b)}\")\n",
    "print(f\"Résultats avec solve {np.linalg.solve(L, b)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Programmer la méthode de Gauss (sans permutation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pivot_gauss(A, b):\n",
    "    \"\"\"Méthode du pivot de Gauss pour résoudre Ax = b.\n",
    "\n",
    "    Args:\n",
    "        A (np.ndarray): Matrice A.\n",
    "        b (np.ndarray): Vecteur b.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Solution du système Ax = b.\n",
    "    \"\"\"\n",
    "    n = A.shape[0]\n",
    "    U = np.copy(A)\n",
    "    y = np.copy(b)\n",
    "    for k in range(n-1):\n",
    "        for i in range(k+1, n):\n",
    "            pivot = U[i, k]/U[k, k]\n",
    "            for j in range(k, n):\n",
    "                U[i,j] -= pivot*U[k, j]\n",
    "            y[i] -= pivot*y[k]\n",
    "    x = trisup(U, y)\n",
    "    return x\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. L'appliquer à la résolution du système $Ax = b$, où $A$ est définie par \n",
    "\n",
    "$$\n",
    "\\left \\{\n",
    "        \\begin{array}{ll}\n",
    "            A_{i,i} = 2 & 1 \\leq i \\leq n \\\\\n",
    "            A_{i, i+1} = A_{i+1, i} = -1 & 1 \\leq i \\leq n-1 \\\\\n",
    "            A_{i,j} = 0 & sinon\n",
    "        \\end{array}\n",
    "        \\right.\n",
    "$$\n",
    "\n",
    "et \n",
    "\n",
    "$$\n",
    "b = (1, 0, ..., 0, 1)^T \\in \\mathbb{R}^n\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_matrix_and_vector(n):\n",
    "    \"\"\"Construit les matrices A et b de l'exemple.\n",
    "\n",
    "    Args:\n",
    "        n (int): Taille de la matrice.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Matrice A.\n",
    "        np.ndarray: Vecteur b.\n",
    "    \"\"\"\n",
    "    A = 2*np.diag(np.ones(n)) - np.diag(np.ones(n-1), k= -1) - np.diag(np.ones(n-1), k= 1)\n",
    "    b = np.zeros(n)\n",
    "    b[0] = 1\n",
    "    b[-1] = 1\n",
    "    return A, b\n",
    "\n",
    "n = 5\n",
    "A, b = create_matrix_and_vector(n)\n",
    "\n",
    "x = pivot_gauss(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapter le programme précédent pour obtenir la factorisation $LU$ de la matrice $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lu_decomposition(A):\n",
    "    \"\"\"Calcul la décomposition LU de la matrice A.\n",
    "\n",
    "    Args:\n",
    "        A (np.ndarray): Matrice A.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Matrice L.\n",
    "        np.ndarray: Matrice U.\n",
    "    \"\"\"\n",
    "    n = A.shape[0]\n",
    "    L = np.zeros((n, n))\n",
    "    U = np.copy(L)\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            U[i, j] = A[i, j] - np.sum([L[i, k]*U[k, j] for k in range(i)])\n",
    "        for j in range(i, n):\n",
    "            L[j, i] = (A[j, i] - np.sum([L[j, k]*U[k, i] for k in range(i)]))/U[i, i]\n",
    "    return L, U\n",
    "\n",
    "A = (np.array([2, 4, -6, -1, -1, 2, 0, 2, 0]).reshape((3, 3))).T\n",
    "\n",
    "L, U = lu_decomposition(A)\n",
    "L, U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Écrire un algorithme permettant de calculer la décomposition de Cholesky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cholesky(A):\n",
    "    \"\"\"Calcul la décomposition de Cholesky de la matrice A.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: Matrice L telle que A = LL^T\n",
    "    \"\"\"\n",
    "    n = A.shape[0]\n",
    "    L = np.zeros((n, n))\n",
    "    for i in range(n):\n",
    "        for j in range(i+1):\n",
    "            sum_coefs = np.sum([L[i][k] * L[j][k] for k in range(j+1)])\n",
    "\n",
    "            if i == j:\n",
    "                L[i,i] = np.sqrt(A[i,i] - sum_coefs)\n",
    "            else:\n",
    "                L[i,j] = (A[i,j] - sum_coefs)/L[j,j]\n",
    "    return L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Tester sur la matrice $A$ de l'exercice précédent, et sur la matrice $B$ suivante\n",
    "\n",
    "$$\n",
    "B = \\begin{bmatrix}\n",
    "    1 & 1 & 1 & 1 \\\\\n",
    "    1 & 5 & 5 & 5 \\\\\n",
    "    1 & 5 & 14 & 14 \\\\\n",
    "    1 & 5 & 14 & 15\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.array([\n",
    "    [1,1,1,1],\n",
    "    [1,5,5,5],\n",
    "    [1,5,14,14],\n",
    "    [1,5,14,15]\n",
    "])\n",
    "\n",
    "L_A = cholesky(A)\n",
    "print(L_A@(L_A.T))\n",
    "\n",
    "L_B = cholesky(B)\n",
    "L_B@(L_B.T)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "method_num",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}