{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7c8c115e-dc7b-41b3-907a-700b3503fc61",
   "metadata": {},
   "source": [
    "## MT09 - TP3 - Automne 2025\n",
    "### Factorisation de Choleski. Méthode de Newton"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c071e68-a309-4a2c-96ad-ecc661d6891a",
   "metadata": {},
   "source": [
    "### 1. Programmation de la factorisation de Choleski sous forme vectorielle\n",
    "Pour une matrice $A$ ***symétrique définie positive***, programmez la factorisation de Choleski :\n",
    "\n",
    "\\begin{align*}\n",
    "& b_{11} = \\sqrt{a_{11}},\\quad b_{i1}=\\frac{a_{i1}}{b_{11}}, \\ \\text{ pour } i>1, \\\\[1.5ex]\n",
    "& b_{jj} = \\sqrt{a_{jj}-\\sum_{k=1}^{j-1} b_{jk}^2},\\ \\text{ pour } j\\geq 2,\\quad\n",
    "b_{ij}  =\\frac{1}{b_{jj}} \\left( a_{ij}-\\sum_{k=1}^{j-1} b_{ik} b_{jk}\\right),\\ 1<j<i.\n",
    "\\end{align*}\n",
    "\n",
    "dans une fonction\n",
    "\n",
    "```\n",
    "def choleski(A):\n",
    "    ...\n",
    "    return B\n",
    "```\n",
    "\n",
    "en remplaçant les sommations par des produits scalaires dans le code (attention aux indices démarrant à 0 en python). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cce744c6-bbce-4d53-a077-55b7d6a10aac",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def choleski(A):\n",
    "    # A compléter ...\n",
    "    return "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "410ae58f-39ff-4930-be8a-fa85dda9478d",
   "metadata": {},
   "source": [
    "Appliquez la factorisation de Choleski à la matrice\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix}\n",
    "2 & -1 & 0 \\\\\n",
    "-1 & 2 & -1 \\\\\n",
    "0 & -1 & 2\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "et vérifiez que $B B^T = A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d1a734-2c49-401e-b6b6-cda6a5b3bf9a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "874bff42-7151-48ab-86eb-fa9d5cd87e92",
   "metadata": {},
   "source": [
    "\n",
    "Question subsidiaire : peut-on s'arranger pour se limiter à deux boucles imbriquées, voire à une seule boucle ? (indication : en vectorisant, puis en travaillant par bloc). A faire s'il vous reste du temps, à la fin du TP par exemple)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b0cc252-3c3f-469b-8bd7-a7043cefa539",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef3e8e0b-ce4e-4675-b65d-495fa303d742",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e971acbf-213f-49d9-a934-8dc77b717e4d",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2. Méthode de Newton-Raphson (méthode de Newton, cas vectoriel)\n",
    "\n",
    "La méthode de Newton-Raphson permet de résoudre numériquement un système d'équations algébriques non-linéaires écrit sous forme vectorielle condensée\n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{x}) = \\mathbf{0}.\n",
    "$$\n",
    "\n",
    "La méthode itérative s'écrit\n",
    "\n",
    "\\begin{align*}\n",
    "& \\mathbf{x}^{(0)} \\text{ donné}, \\\\[1.3ex]\n",
    "& \\mathbf{x}^{(k+1)} = \\mathbf{x}^{(k)} - [J(\\mathbf{x}^{(k)})]^{-1}\\, \\mathbf{F}(\\mathbf{x}^{(k)})\n",
    "\\end{align*}\n",
    "\n",
    "en partant de $\\mathbf{x}^{(0)}$ (initialisation, dit 'initial guess'). La matrice\n",
    "$J(\\mathbf{x}^{(k)})$ est la matrice jacobienne de $\\mathbf{F}$ au point \n",
    "$\\mathbf{x}^{(k)}$ (voir cours/poly).\n",
    "\n",
    "En pratique, on procède en deux étapes: on résout d'abord le système linéaire\n",
    "\n",
    "$$\n",
    "J(\\mathbf{x}^{(k)}) \\, \\mathbf{w}^{(k)} = -\\mathbf{F}(\\mathbf{x}^{(k)}),\n",
    "$$\n",
    "\n",
    "puis on fait l'affectation\n",
    "\n",
    "$$\n",
    "\\mathbf{x}^{(k+1)} = \\mathbf{x}^{(k)} + \\mathbf{w}^{(k)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e48e8663-97bc-4cf7-8557-cb18a0943057",
   "metadata": {},
   "source": [
    "Ecrire une fonction `python`\n",
    "```\n",
    "xout, kout, Fout = newtonRaphson(x0, foncjac, kmax, tol)\n",
    "```\n",
    "qui résout l'équation $F(x)=0$ avec la méthode de Newton-Raphson avec ```x0``` l'initial guess, la fonction `python` générique\n",
    "```\n",
    "F, J = foncjac(x)\n",
    "``` \n",
    "étant supposée retourner la fonction $F$ ainsi que sa matrice jacobienne $J(\\mathbf{x})=DF(\\mathbf{x})=(\\partial F_i/\\partial x_j)_{ij}$ au point trouvé, $k_{max}$ est le nombre max d'itérations de Newton,  $tol$ est la tolérance de convergence : on fera le test\n",
    "\n",
    "$$\n",
    "\\frac{\\|\\mathbf{F}(\\mathbf{x}^{(k)})\\|}{\\|\\mathbf{F}(\\mathbf{x}^{(0)})\\|} > tol\n",
    "$$\n",
    "\n",
    "(si oui et si $k<k_{max}$ on continue, sinon on s'arrête).\n",
    "\n",
    "L'utilisation de ```newtonRaphson()``` suppose que la fonction ```foncjac()``` a été préalablement définie, en dehors de la méthode de Newton.\n",
    "\n",
    "La fonction ```newtonRaphson()``` doit retourner $x_{out}$ la solution trouvée, $k_{out}$ le numéro de dernière itération et $\\mathbf{F}_{out}=\\mathbf{F}(\\mathbf{x}_{out})$ le 'résidu vecteur'.\n",
    "\n",
    "NOTA BENE : on remarquera que ```python``` autorise le passage de fonctions (ici ```foncjac()```) en arguments de fonctions `python`. Il s'agit bien de l'objet 'fonction', et non de l'évaluation de la fonction en un point.\n",
    "La fonction ```foncjac()``` pourra ensuite être évaluée aux itérés de la méthode de Newton-Raphson dans ```newtonRaphson()```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ac627cd6-9fe7-4e63-a5b7-ffae92f03876",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "\n",
    "def newtonRaphson(x0, foncjac, kmax, tol):\n",
    "    k=0\n",
    "    xk = x0\n",
    "    # ...\n",
    "    return xout, kout, Fout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25e5ddbf-e125-4f15-8173-8fc19c821ab3",
   "metadata": {},
   "source": [
    "Vérification : écrire la fonction ```foncjac()``` et tester la fonction ```newtonRaphson()``` pour la fonction à 2 variables\n",
    "$\\mathbf{x}=(x_1,x_2)$ :\n",
    "\n",
    "$$\n",
    "\\mathbf{f}(\\mathbf{x}) = \\begin{pmatrix} x_1^2 - 2 \\\\ x_1^4 + x_2^2 - 5 \\end{pmatrix} \n",
    "$$\n",
    "\n",
    "en prenant $\\mathbf{x}^0=(4,2)$, $k_{max}=16$, $tol=10^{-15}$. Quels sont les zéros de $\\mathbf{f}$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46104d6-a672-430f-9ffa-50d8da16315a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Exemple\n",
    "def foncjac(x):\n",
    "    # ...\n",
    "\n",
    "# ...\n",
    "\n",
    "xout, kout, Fout = newtonRaphson(x0, foncjac, kmax, tol)\n",
    "print('xout = ', xout, \"\\nkout = \", kout, \"\\nRésidu = \", Fout)\n",
    "# Verif\n",
    "print('Erreur =', np.sqrt(2) - xout[0], 1.0 - xout[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eff893a-0cd2-45df-808b-d68e578eae56",
   "metadata": {},
   "source": [
    "### 3. Problème inverse, identification paramétrique, solutions de systèmes d'équations algébriques, \n",
    "\n",
    "Un modèle simplifié de croissance cellulaire avec saturation donne l'évolution de la densité de cellules $\\rho(t)$ au cours du temps $t$ :\n",
    "$$\n",
    "\\rho(t) = \\rho_M\\, \\frac{e^{\\alpha t}}{\\frac{\\rho_M}{\\rho_0}-1 + e^{\\alpha t}},\n",
    "\\qquad t>0,\n",
    "$$\n",
    "\n",
    "où $\\rho_0$ est la densité initiale, $\\rho_M$ est la densité maximale $(\\rho_M\\gg \\rho_0)$, et $\\alpha$ est le taux de croissance cellulaire dans la phase exponentielle. \n",
    "\n",
    "Tracer la fonction $t\\mapsto \\rho(t)$ sur l'intervalle $[0,10]$ pour les paramètres\n",
    "$(\\rho_0,\\rho_,\\alpha)=(3, 1000, 1.2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02bc00f4-e5d3-46d8-a7d0-ccfe8a96fd5b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c1d12d01-4a45-4259-a237-303f12325903",
   "metadata": {},
   "source": [
    "Les quantités $\\rho_M$ et $\\alpha$ sont d'intérêt biologique, mais ne sont pas des observables (ne peuvent pas être mesurées directement).\n",
    "Expérimentalement, on peut compter le nombre de cellules sur un élément de surface, on peut donc mesurer la densité à un temps de prolifération donné.\n",
    "\n",
    "À partir de trois mesures de densité de cellules à rois instants différents, on peut\n",
    "espérer déterminer les trois paramètres $\\rho_0$, $\\rho_M$ et $\\alpha$ du modèle.\n",
    "\n",
    "Supposons que les trois mesures sont les suivantes :\n",
    "\n",
    "| Temps $t^k$        | Densité $\\rho(t^k)$ | \n",
    "| :---------------- | :------: | \n",
    "| 2.0               |   32.104 | \n",
    "| 4.0               |  267.726 | \n",
    "| 8.0               |  977.987 | \n",
    "\n",
    "À partir des mesures, déterminer les paramètres $\\rho_0$, $\\rho_M$ et $\\alpha$.\n",
    "Pour cela, dans un premier temps, écrivez le système d'équations  $\\mathbf{F}(\\mathbf{x})=0$ à résoudre.\n",
    "\n",
    "En utilisant `scipy.optimize.fsolve()` (voir l'aide associée), résoudre le problème. On pourra utiliser l'initialisation :\n",
    "\n",
    "```\n",
    "X0 = np.array([1.0, 2000.0, 1.0]) # X0=(rho0, rhoM, alpha) initial guess \n",
    "```\n",
    "\n",
    "Affichez les valeurs de  $\\rho_0$, $\\rho_M$ et $\\alpha$ ainsi obtenues.\n",
    "Arrondir les valeurs à deux décimales en utilisant `np.round(...,2)`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "a1825b00-e204-4dce-88ff-4efbc8e1f6d7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'nfev': 25, 'fjac': array([[-0.15555545, -0.98174627, -0.10943838],\n",
      "       [-0.0258552 , -0.1067026 ,  0.99395476],\n",
      "       [ 0.98748874, -0.15744463,  0.00878508]]), 'r': array([-6.54560640e+01, -1.61890211e-01, -7.87285652e+02,  9.22863106e-01,\n",
      "        1.00575418e+02, -6.18274814e+01]), 'qtf': array([ 3.00091968e-08,  9.62977997e-08, -8.41232527e-09]), 'fvec': array([ 5.90105742e-11,  3.48677531e-10, -1.95541361e-11])}\n"
     ]
    }
   ],
   "source": [
    " def F(X):\n",
    "    # ....\n",
    "    return np.array([f1,f2,f3])\n",
    "\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "X0 = np.array([1.0, 2000.0, 1.0])\n",
    "\n",
    "Xsol, infodict, ier, mesg = fsolve(F, X0, full_output=True)\n",
    "\n",
    "print(infodict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d587d448-8b7c-4bf3-b020-5e243ab73c45",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho0, rhoM, alpha = Xsol\n",
    "print('rho0 = ', np.round(rho0,2) )\n",
    "print('rhoM = ', np.round(rhoM,2) )\n",
    "print('alpha = ', np.round(alpha,2) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b491b1d9-dd6b-4562-8455-c4faedf3eb75",
   "metadata": {},
   "outputs": [],
   "source": [
    "?fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "379e1ea8-2b07-41d9-8a54-5120d001baf1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "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.11.5"
  },
  "toc-autonumbering": false,
  "toc-showcode": true,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}