{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compléments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing libraries used in this book\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# scipy.linalg.lu : LU decomposition\n",
    "from scipy.linalg import lu\n",
    "# scipy.linalg.cholesky : Cholesky decomposition\n",
    "from scipy.linalg import cholesky\n",
    "from scipy.linalg import hilbert\n",
    "\n",
    "from IterativeMethodsLib import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans le cas d'une matrice mal conditionnée, on peut essayre d'utiliser la factorisation LU ou de Cholesky comme préconditionneur. Mais cela ne marche pas forcément ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 4\n",
    "A = hilbert(n)\n",
    "\n",
    "D = np.diag(np.diag(A))\n",
    "\n",
    "MaxIterations = 1000\n",
    "\n",
    "x = np.ones([n,1])\n",
    "b = A.dot(x)\n",
    "alpha = 0.5\n",
    "\n",
    "# we do not keep track of all the sequence, just the last two entries\n",
    "xk = b\n",
    "rk = b - A.dot(xk)\n",
    "\n",
    "RelativeError = []\n",
    "residualNorm = []\n",
    "\n",
    "# We rewrite Richardson for the sake of the example: \n",
    "for k in range(MaxIterations) :\n",
    "    \n",
    "    zk = np.linalg.solve(D,rk)\n",
    "\n",
    "    xk = xk + alpha*zk\n",
    "    rk = rk - alpha*A.dot(zk)\n",
    "    #rk = b - A.dot(xk)\n",
    "    \n",
    "    \n",
    "    RelativeError.append( np.linalg.norm(x-xk) / np.linalg.norm(x) )\n",
    "    residualNorm.append( np.linalg.norm(rk) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.plot(range(MaxIterations), RelativeError, 'b:.')\n",
    "plt.plot(range(MaxIterations), RelativeError, 'b:.',range(MaxIterations), residualNorm, 'g:.')\n",
    "\n",
    "plt.xlabel('n'); plt.ylabel('err');\n",
    "# plt.xscale('log') \n",
    "plt.yscale('log')\n",
    "plt.grid(True)\n",
    "plt.legend(['Rel err','res'])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "A = hilbert(n)\n",
    "\n",
    "# epsilon = 2\n",
    "#A = np.array([[epsilon, 1, 2],\n",
    "#              [1, 3, 1],\n",
    "#              [2, 1, 3] ])\n",
    "\n",
    "# B = np.diag(np.diag(A, k=1),k=1) + np.diag(np.diag(A, k=0),k=0) + np.diag(np.diag(A, k=-1),k=-1) \n",
    "# L = cholesky(B, lower=True)\n",
    "P, L, U = lu(A)\n",
    "\n",
    "MaxIterations = 30\n",
    "\n",
    "x = np.ones([n,1])\n",
    "b = A.dot(x)\n",
    "alpha = 0.5\n",
    "\n",
    "# we do not keep track of all the sequence, just the last two entries\n",
    "xk = b\n",
    "rk = b - A.dot(xk)\n",
    "\n",
    "RelativeError = []\n",
    "residualNorm = []\n",
    "\n",
    "for k in range(MaxIterations) :\n",
    "    \n",
    "    # y  = np.linalg.solve(L,rk)\n",
    "    # zk = np.linalg.solve(L.T,y)\n",
    "\n",
    "    y  = np.linalg.solve(L,P.T.dot(rk))\n",
    "    zk = np.linalg.solve(U,y)\n",
    "\n",
    "    xk = xk + alpha*zk\n",
    "    rk = rk - alpha*A.dot(zk)\n",
    "    #rk = b - A.dot(xk)\n",
    "    \n",
    "    \n",
    "    RelativeError.append( np.linalg.norm(x-xk) / np.linalg.norm(x) )\n",
    "    residualNorm.append( np.linalg.norm(rk) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.plot(range(MaxIterations), RelativeError, 'b:.')\n",
    "plt.plot(range(MaxIterations), RelativeError, 'b:.',range(MaxIterations), residualNorm, 'g:.')\n",
    "\n",
    "plt.xlabel('n'); plt.ylabel('err');\n",
    "# plt.xscale('log') \n",
    "plt.yscale('log')\n",
    "plt.grid(True)\n",
    "plt.legend(['Rel err','res'])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}