path: root/PVCM/cama/fr/ma61 Système matriciel non linéaire.ipynb

{
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Système matriciel non linéaire\n",
    "\n",
    "Si la matrice A dépend de ${\\bf x}$, ce qu'on écrit $A({\\bf x})$, alors le système\n",
    "matriciel \n",
    "\n",
    "$$\n",
    "A({\\bf x}) \\, {\\bf x} = {\\bf b}\n",
    "$$\n",
    "\n",
    "n'est pas linéaire.\n",
    "\n",
    "#### Exemple\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 &  x_1 \\\\\n",
    "2x_1 & -x_2 \\\\\n",
    "\\end{bmatrix}\n",
    "\\;\n",
    "\\begin{bmatrix}\n",
    "x_{1} \\\\\n",
    "x_{2} \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "b_{1} \\\\\n",
    "b_{2} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "donne le système suivant qui n'est clairement pas linéaire puisqu'on a des\n",
    "carrés :\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x_1 + x_1 \\, x_2 &= b_1 \\\\\n",
    "2 \\, x_1^2  - x_2^2 &= b_2 \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Comment résoudre un tel système ? A-t-il une solution ?\n",
    "\n",
    "On voit bien qu'il y a toujours autant d'équations que d'inconnues donc il y a\n",
    "de l'espoir. Par contre on sait que les équations du 2e degré ont plusieurs solutions donc il devrait en être de même dans certains cas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### La méthode du point fixe\n",
    "\n",
    "La méthode du point fixe consiste à appliquer l'algorithme itératif suivant :\n",
    "\n",
    "$$\n",
    "{\\bf x}^{k+1} = {\\bf g}({\\bf x}^k)\n",
    "$$\n",
    "\n",
    "dans le but de résoudre l'équation ${\\bf g}({\\bf x}) = {\\bf x}$.\n",
    "\n",
    "Si vous ne connaissez pas cette méthode voici un exemple d'utilisation en 1D avec\n",
    "la courbe $g$ dessinée pour bien comprendre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEBAAIDAQAAAAAAAAAAAAAAAQMEAgYHBf/EAEoQAAEDAgIGBwUEBgcIAwEAAAEAAhEDIQQxBRJBUWHRExdTcZPS4RQikZKxMlSBoSMkNFLB8BUzQkRyc5QGYoKDorKzwjVDYwf/xAAXAQEBAQEAAAAAAAAAAAAAAAAAAQID/8QAHREBAQEAAgIDAAAAAAAAAAAAAAECAxExURIhIv/aAAwDAQACEQMRAD8A8/REQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERARdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lTq50x2+C8R3lQdQRdv6udMdvgvEd5U6udMdvgvEd5UHUEXb+rnTHb4LxHeVOrnTHb4LxHeVB1BF2/q50x2+C8R3lU6udMdvgvEd5UHqBNkJsUJEZoSIzVRZUBsEkb0BEC6gA2QHNARvQEXuqE3SbpInNJE5oBNwhKEiRdCRvQCbFWVCRGaSN6gA2CA2QEQLoCN6oxVcTSoR0hILjAAaSfgFyq1qdBhfVdqtEXhYaBD8biqjs2FtNt9mqHfV35BYNLP1W0YN9Y6o3kscAPiQoN8OBgi4IsVSbhcaYbTYxgNmtgLkSJF1QJQmxQkb18z22q7EPaKha9tYM6GBBaXRM5zqgu+CD6kqA2CSN6AiBdQAbIDmgI3r5w0hWqY/GYWjSk0SwCoSNVstm95Ko+jN0m608PjA3B4Z1d+vWqU2uIY2S4xcwNiuKxvRB3Rhuu12r75hv2dbPuUG2TcISuNOoKlNj4jWAMHYseLrOpUmmnqlzntaNbK5CDMTYqyoSIzSRvQAbBAbLDhcS3EMe4CAyo6nnnBhZgRvQAc0m6Ai91NZutEiTxVFm6E3C4uqMa9gc4AuMN4mJ/gVyJEi6AShNitLD4upVxAlzOicwvAGbRIDSTxufwW7rAtkEEEKCyoDYJI3oCIF0AGyA5oCN6Ai91Qm6TdJE5pInNAJuEJQkSLoSN6ATYqyoSIzSRvUAGwQGyAiBdARvVAHNJugIvdJE5oE3Qm4SROaEiRdAJXz9P46ro/RFbE0A3pGFsawkXIC+gSN6wY6izE4V1J7G1GmJa4Ag3UvhrFk1Lrw6von/aXSGNdV6QUQGNkarDnc7+C2tIacx+EwlR/Rsa9lhrNkH7EHPc5fQo6Kw9Ik08JRYXCCWsblKzVcIK9Isq0GPa65DmtO7kFmS/Hp23ycd5ZqT8+nxv9mv9ocbpTSLsPiRS1BSL/caQZkceK7ROa+do/R9DCVzUp4alSOrq6zWgH8l9GRe6uZZPtnn3je++OdQOSpyKhAjJCBBsq4qgyCQNwUAECyKoyQbVABuCAC9giLtTapAnIJAnIIKcwhUIEiyEDcEFORRQgQbKwNwQBkEGSgAgWQAbgg0m16WExOIbiKjaYqPD2FxjW90CO+QbdywaTeX1MBXYJArQxr2ke+WuDSeF19QAXso5jHOGs0HVMiRkUHzcXj6lDBF1R2q6nWbSq1GAe40ke/fIQQeCx4PSLq+h2YoVw/WrloeWgSwVdXLuX1OhpdManRt1y2C6LwtPSdDD0sHiMR0DNYN13ENu4C8fkgy4DFuxtJ9Tow1mvDDMyIH5rXrQ7TlF1tZjS3Vi8EEl3xa0fHet+kymyixtJjWsAAa0CAAuRa25gTGcIMdZ9droo0mPtm9+qPoVlGQSBuCACBZBRkvmYGi9umNKvc0htR1MtO+GAL6QA3BABewQfFptr0sZh9am8GlQY3UZMvIDhG6Ntzb8VNINcMWKr6bqh6WmG05JZrEESd94vshfbgTkELQSLBB8yvj6goUJcWk1TTr1KbZ6OAdkGJIGe9aWH0pia7qdPEt1XUaw6YdGZDZgOnIXjft2BdgIEiyaoGwINPA4t+JNcVejYW1HNYwH3gAYk32rJpGt7PgatTX1AAJf+4CYLvwF/wAFnLG3OqJ3wqWtLSC0EHMEIrr+hsY4YagGObSoufVrVHViZDS8kD3jMwQfxC7A0gtBFwdq49Gw6pLGki4tkuQAjJEfIp6RxILw99Iu13MMNgUj0ga0m+0XjguTNIk0nVX06b8Qxs0x9mWloJcc4C2cHSpV8NWY9jKjHVqkhwkH3ytvUYZBaCIjLYg+SMcdIVMEaBbRqa7nOZUEuaNV0SJGYuvp4mo+lh3vpsD3tEht7/AFcyxjnNloOrcWyKroGyUHw8G8UtHYPDOpVP1oBlSrVgCNQnfuEBZzjn0XYqmzVdRwrdf3Rm0CNUcQQZOyw7tzBUtXR+GZVZD202SCMiAFslojIINXAYirVafaH0g9xljGm4bxuVtjILg2lTZOoxrZzgQuQAgWRVGSDaoANwQAXsiLtTapAnIJAnIIKcwhUIEiyEDcEFORRQgQbKwNwQBkEGSgAgWQAbggo2ptUAF7BIE5BBdqHMKQJyCECRZBShyKhA3BCBBsgqDIJA3BQAQLIqjJTegA3BIF7IhAj1Qix5oRbMoRY3QWAoBYc1Y4lQCwugAWQDNAOJQDO5QIv6pF0i+ZXEPaQXB4gTJkWQciLjmhAWjiNKUKGkKWCJca1QawA3TE/EhbVesyhRqVHkkU26zgM4QZCLHmrAWthcZTxjahp6w1DF4vxsswexzi0PBLTBAIsUHICw5oBZYn12UqlOm94aakxJAyWUDiUADNIv6rhUcabHOA1jIgSBKlCoaocSNUgxGsD9FRki641KbajDTeJa4EETsXKL5lYcRVdTdTZTbrve7ImIAzOX8yFBl1AGaokRYLA7D1i0/rdQcQ1vJY6WPp1qrGs1jTfZr8pOcQRwPwW4RbMqjWODqHPHYn/o8qrMK4NvisQ7vI/gFsxxKx1KtOi0GrUDBGbjGQlQYxhjH7RX+YclPZJzr1yP8AHH0WwLiZWngsb7TiMay2ph6optcDIdLGu+phBk9jZP8AW1/FcnsgDvdr1x/zCfqtgXNisPS1DiCxtJ2q0wXuIA/DeghwxkfrFb5hyXF2EfMjGYgcJb/ELZNrkoRbNBrHC1o/bq/ys8q5ez1o/a3/ACt5LK+oxrgx1QBzvsgkSUe9lPV13husdUTtO5BiFCtA/Wn/ACjkgoVvvT/lHJWpiaVKoym95DnRADSYkwJ3X3rMBxKow0aD2VC52IqVB+6Q0D8hKzRf1QDO5SL5lAi6ECRzSL5lcOlpl4YKg1tbVidsT9FBzICEWPNCOJQi2aCwFALDmo9zWN1nvDRvJhYn19Ws2kL+7rPMgBg3/wA7igzAWQDNcab2v1gx4dqmHQQYO5YxVe3F9C8ANc3WY4HOIkG3FUZov6pF0i+a1q+NpUcR0PvvqBmu5rACWtmJKg2SLjmhAXFjm1GMex0tcJB3hcMVXp4ag+tVdDabdY5ZBBlIseasBYaeIo1gOhrNqBwJBYQQf5lZo4lBALDmgFkAsLrWo4guNR9TUZQBIY4ugmJkxusg2QM0i/qsOJfUp0H1KIDi33oO0bhxWVpDgHNMggEFBYuhFxzWI16YxLaGuekImAMu/cssZXQCAhFjzQ96jyGsLnOgDMlBygKAWHNauNq16TqXQlkOdEOElxtYfhJnZC2gLC6ABZIF0A4lIzuqBmM0MwhmEMwVBYKgmArfgoJgIAlBN0ExsQTfJBCCZE7F8rRuFZozFuwTapqdKwPuLiIbfvAHwK+teViOHpnFNxGr+kDdWZOX8k/FBpVNDsraRp4yvU1qgoOouAEAyQQRuIhcn6Ka9h1q1Q1XMLHPJJ1pG4nK8wvoGZCGeCDSwGjzgGBramuOjDXEtguIJg27/otavoQ4ioOlxP6IOcdRrIJkznOeyd1l9YzBS/BB83B6IbhegDaxqNpT/WNBJke9fiQCvpCYQTAQTGxBjr0RiKFSi4w14gkC4WOlheiwhwwcC3VLQS3Zx3rYE3yS8qj5uNwOrox1Ntf2YMmoXUaYAPDV2rno6i44GgH1Huq0S9ge/MwS2THct17BUaWPa1zXAgg5EJTptpMaymxrWjIBQatHR7cPhmUKVQ6rageJAnOTlvM/FbhmNiGeCGYVFgrXxeEbi6TGvc4ajg8au8ZfndbF+CgmAoPngVtKYKtQxFF+FDjF76wnI/xjfmsf9A4c0TSNWtDnB7yCAXOAiTbaDFtwiIX1BMbEE3QfJp6Ao0i4UMTiKLHVDULaboEmf4W/AbVt0MC6jinVva67y+NZriNUwIyi34LbvKXlBjxFFuIouo1LseII3hatDRxo4L2R1UVaIaA3XbrOB3ySZve4W8ZkIZ4IND+jB03SOxD3DVYNQtaB7hJGQtn+QVOCqV20alWqadem9zg4NaSAZ92TOw/kt4zCX4INF9HFGlSquqNc+lrPNPUmdwFxcZTxWxSp1hSpB9aXNaA86o942k/X4rMJgIJjYgjA4CC6TvhW8oJvkl5QLyvm4zCYmjiHYzAFr677OZU+yRGQ3XAvyEfSvOxDMhB89uNxhqODsC9rGsaZMXP9rInJbOMwoxdDonuIbIJhoMx3ghZzPBDMKj5ml8DQrYFtN7yypSY7oYdq31SMhwWycMauu5tTVZWphrmhtxbYfxWWvh6eJaG1WawH+8R+HdwWUTGxQfPfSxmDo0KOj6dJ9JrocKhu1t7D6fzK2dSo/EUahAaGtdIJuCYt9VnExsQTfJULytTF6OpYqq57nOb0lPo6gaB+kbMgG3E5bytu8pedig+ezR9dtWm72mBTPugCQ4AEXHcfiJ4LgNDwGziqp1dYA7YcZOfGD+AyX0zMhDPBB89+iW1KNKnVxFV7aTy8CwknfC2cFhn4WiKJrOqsb9gvkujiZus5mFb8EEEwFpP0cKtBlGrVJbTfr0yGiQRlO9bomAgmNiDTGHxs4nWxTHteCKbTTjUta/x/JWlVr08b7O7DkYcMGpWBmTuI2ZLbE3yS87EGjjdGnFNrMbiqtFtaC8MDTMQNoNoGS1cHoEYLD1MPSxlbo6hm5u0yXW4STIX2LzsQzIQfMo6H6BxezGVw51XpjcQXEQfw4LXp6MdiKeLpvxLqhGI1tZ7LO90WIBEi+/ZwX2zPBDMIMXR1CaRc9ssB1tVka1u+w9FlEwFb8FBMBAEpBugmNiXuqBJjJCbZFJt6ITY8kC+5AbCxVnv+CgNhyUAHggOdkBt6IDnyQNuSTfJJv6JN1RrsxtKrVLKWtU1Xapc1pLQd05LYJO4r5uHbWw+Lq0xTq9CausxrGNDYIuSSd82zWOtjNI0gCaEjpw2dW3R3vAvJiPxCg+sTbIrDXxTKD6LHh2tWfqNgbYlfNbpDSjKWKdi8DTpCnTeabmuLgSBadwKOGMxFSniujpvfRZVFKAQ3XIEEzeMx+e2wfYBsLFAeCw4So9+Hb0jXhwABc5urrWzjYswNtvwQAc7JtyQHPkk39ECb5ITcWSb+iE3HJUDO5CbZFCe/4ITY8kC+5AbCxVnv+CgNhyUAHgsNbEtosc43DHAPg/ZB2lZgbei+XorBHDY3HVXVC6pVqA1LWPutj4XHcg3cNjKWKfVFB2u2k4NLgZExMD4hZ5vksVCn0IeJkucXTEZm3wEBZZv6IBNxZDO5CbjkhPf8EAm2RS+5CbHkrPf8EEBsLFAeCA2HJAbeioA52TbkgOfJJv6IE3yQm4sk39EJuOSAZ3LC/GUGg6zw0AuBJtECSrRr9Mavu+6x+q1370RP5yPwWFuFZ7dXxBpgFzAwEZkbT9Pgg5VNI4WnT1zVBadrRrfTZcI7F6ldzXarKVOA57jHvG8fzvWppXR9KpgOjp0HVngFjNZxcW622+cWW5WweHxBJqtc4OADmydV0ZSMioAqvGMY0O1qVWmXDKxBH1n8lsA52WJtFgfTcC6abS0Amc4zm+xZQc1Q25JN8km/ok39FAJuLIZ3ITcckJ7/AIIBNsil9yE2PJWe/wCCCA2FigNskBsOS41C3on6wJEGRqz+SoCo0uLQQXbpC5bcl87RlMjDtfWwzKdW0Um0w3UixP8AH+b7mIdWFJxw7QamwET/ABH1UGWb5ITcWK0cRisVSxWHij+gdqh7rWJMRnvjKc1q1Nekyg0U3tNStFd3vNIP70g3Gz8RuQfXJyCpNslo1ejfSoVqQc9zKjQx5BJIJAN9ogrZxTntw7jSBLhBiJkTf8kGW+5AbCy1qPS0q3Rw51DVljtX7PAkmT8Fh0fin1HVBiHgPe/9HTDSIaGjeJzm6DfB4JJvZaWjMSa7KhfW13axIYQA5jdgIGS3Zz5KgSIzQkRmhyVORUEkb0BEC6qDIIqAjegIvdUZINqIkic0kTmrtTaghIkXWOv0haOhexrgf7TZBHxWU5hCg0MXQxdbB1qZxNP32ObDKUZjiStylIpMDyC6BPeuZyKIICIF0BG9UZBBkggIvdJE5qjam1BJE5oSJF1dqHMIISN6EiM1ShyKCSN6AiBdVBkEVARGa16JHtuJuMmn8itkZLVo/wDyeK/y6f8A7IjZkTmkic1dqbUEJEi6EjeqcwhQQkRmkjeqciiCAiBdARvVGQQZIICL3SROao2ptQSROawYuqRqUqToq1DAP7o2u/D6ws7nBoLnEAASSVrYUGo92JeCDUswHMMGXxz/AB4IM1NjKNJlNlmtAAuuOIrto05+091mMGbimIxDaWqwDXqu+ywZnjwHFcaNAtca1Yh9ZwidjRuHD6oMWO6Sk04qi5gqNplpaROtOX4z9VG0a9PEUWU3RhmMAJ6SLjhq3+K2qtJtXU1iYa4OjfGSyDIIrRdRrjSHtLHDUA1DTDvtj948Rs/HhG6CL3VGSDaiJInNJE5q7U2oISJF0JG9U5hCghIjNNYAZhU5FaukWg4Ulz6bWNILxUMNcNxO5BsNc0tBDgRwKoI3rT0fiG12kU6TaIb9unEOkgEHuzzC54vGeyvosFPXNZxaIcBB2fhx7t6DK6m19ZlXWIcyQIOYOYPwHwWSROYWLE1jQouqNpPqkEDVYJKwNxRbg6mJqua8C+qz+zsi9yUHLHUqtd9DonN1Gva54Jiwc0z+X5raJEi4Wrgw2q1mIp1Kxa4SA90yO7YtvaEGOmynRpNp0xDG2AmYXMkRmqUORQSRvXENZrB8DWiJ4LmgyCKgI3hJF7qjJTeiBAjJCBBskCPVCLHmgsDcFABAsrAUAsOaAANwQAXsEAsgGaBAnIJAnIJF/VIugECRZCBuCEXHNCAgECDZWBuChFjzVgIIAIFkAG4IBYc0AsgAC9gkCcggGaRf1QIE5BCBIskXQi45oBA3BCBBshAQix5oLA3BQAQLKwFALDmgACMlr0QDjcTYWDR+R5rYAstfDj9Zxf8AjH/a1BsQJyCQJyCRf1SLoBAkWQgbghFxzQgIBAg2VgbgoRY81YCCACBZABuCAWHNALIAAvYJAnJAM1ixNXoaRc0azzDWNnMnJBixAGJrjDNHuD3qp4bG/j9O9cq1f3+hw7Q+rtJ+yzv5Z/VYKDX1G9HQedUkmpiNr3bdXnsyE7NylRp0WNZTbqtHH80HGjh20QTJfUcfee7M+nBZSBBshAQix5oLA3BQAQLKwFALDmgADcEAF7BALIBmgQJyCQJyCRf1SLoBAkWQgbghFxzQgIBAg2SBtAQix5qwEEDRnAXCpSZVpuY5oMiLiVzAsOaAWQYMFRq0qOriK3TPm7tWFmDGgmGgTc2VAzSL+qDTOBLdJDE0XU2NLSHtLCSZIm82yGxbhAkWSLoRcc0AgbghAg2QgIRY80FgbgoAIFlYCgFhzQABuCQL2QCyQLoBFsyhFjdDMZoZhBY4lQCwurBUEwEADiUAzuUEoJugRfMpF8yl5S8oBFxdCOJQzIQgoBFjdWOJUMwrBQQCwugHEoJgIJQAM7lIvmUE3S8oEXzKEXF0vKGZCoEcShFjdCChmFBY4lQCwurBUEwEADitfDj9Zxd//sH/AGtWwJha2EJOKxot7tUD/oaf4oNmL5lIvmUvKXlAIuLoRxKGZCEFAIsbqxxKhmFYKCAWF0A4lBMBYq1cUYF31HTqsaLu/neqOb3tpMc+o8Na25J2L5bA/SmOJqtqU8NRsxpEdIdpO7dHfO5Z6zaoAq1S12IcdWjTH2WO38YznhZbNCm+i5tFoHRNYPeOZKgyhsQBYAZKkXF0vKGZCARxKEWN0IKGYQWOJUAsLqwVBMBAA4lAM7lBKCboEXzKRfMpeUvKARcXQjiUMyEIKARY3VjiVDMKwUEAsLoBxKCYCCUADO5SL5lBN0vKBF8yhFxdLyhmQqBHEoRY3QgoZhQWOJUAsLqwVBMBAA4lIzuglIN1QMwhmChJjJCbZFBb8FBMBL7kBsLFQBMbEE3yQHggOdlQvKXnYm3JJvkgGZCGeCE3FkM7kAzBVvwUJtkUvuUATAQTGxAbCxQHgqAm+SXlAc7JtyQLzsQzISb5ITcWQDPBDMFDO5CbZFBb8FBMBL7kBsLFQBPBamDn2vH5f1w/8bFtg8CtXBz7Xj7G9Yf+Nio2ryl52JtySb5IBmQhnghNxZDO5AMwVb8Fwq1G06bnvIa0CSSbBa/6bFi+tRoHZk9/lH59ygr8Q97jSwwa54s55+yzv3nh9Fyp06eFpvqvdLol9V5uY/hwXPWpYXDiQKdNg2WAWFlN+JeKtdhbTaZp0jv/AHncdw2d+QXDMfUquxNRsEiKbTm1vM8lWk/0lVH/AOLPq5bA22XAFnTmw6QNE74v6qjnediGZCTfJCbiyAZ4IZgoZ3ITbIoLfgoJgJfcgNhYqAJjYgm+SA8EBzsqF5S87E25JN8kAzIQzwQm4shncgGYKt+ChNsil9ygCYCCY2IDYWKA8FQE3yS8oDnZNuSBediGZCTfJCbiyAZ4IZgoZ3ITbIoLfgoJgJfcgNhYqAJjYl7oDwSTeyoTb0Qmx5ISIzQkRmgs9/wUBsOSSN6AiBdQAbeiA58kBG9ARe6BN/RJv6JInNJE5oBNxyQnv+CEiRdCRvQCbHkrPf8ABQkRmkjegA2HJAbeiAiBdARvVAHPkk39EBF7pInNAm/ohNxySROaEiRdAJ7/AIITY8kJG9CRGaCNqMd9lwPcuNWs2jRdUdMNE2Ga+fUwHROcyiGOpYmtrVWuZIbmQY3SGg/HatfTmBp1W4GocTUpvoVGhjWfZdcEyO5p+m1QfUOJLMOyq6k9usWhzDEskxeLWlTCH9Ni/wDO/wDRq+fV0fjGCq3D1qL21w/pNdpBBMwBnaSbcfwW3htZ1fHgVQ39KIt9k6jfig3Zv6JN/RYg1/TF3Tyz92B9VAyqGBpxEumdbVGW5UZibjksVfEtpFrA0vqu+yxuZ48BxWGrinVZbhXMESDWf9hp4bz/ADKw0KuHo1C+niXYl7oDgwB5J3kgW7phQbDaDnuFXEu13i7WAe6zu3nifyXOrimsf0dNpq1onUbs4k7AsQZiaxcalU0aZuGAgu7p2fn3rJg6VGhhwyg4OaCZdMkmbydpQSlhyXNq4h3SVR9m3us7h/HNbANtqAiBdARvVAHPksYcfaXDo4bqCH77m3871kBF7rGHH2lw6RurqD3doMm/87kGSb+iE3HJJE5oSJF0Anv+CE2PJCRvQkRmgs9/wUBsOSSN6AiBdQAbeiA58kBG9ARe6BN/RJv6JInNJE5oBNxyQnv+CEiRdCRvQCbHkrPf8FCRGaSN6ADYckBt6ICIF0BG9UAc+STf0QEXukic0Cb+iE3HJJE5oSJF0Anv+CE2PJCRvQkRmgs9/wAFAbDkkjegIgXUAG3ok5oCN6SL3VA5KnIqECMkIEGygqDIJA3BQAQLIqjJBtUAG4IAL2CIu1NqkCcgkCcggpzCFQgSLIQNwQU5FFHQGkxsWsMfQ/cxH+mqeVBtDIIMlqDHUIHuV/8ATVPKuTcW0g6uHrn/AJZH1QbI2ptWqKuJfPR4QN/zagH/AGympi3O+3QpjaAwuPxkfRBtbU2havstQul2Lq9zWsA+kq+xi36av86DZKHJazsEw/aq1z/zSPohwGGgh1Mv/wAx5d9SgzVa1KiJq1GMG9zgFr+2YGs1p6WlWAMjV9+DvsstLCYah/U0KVP/AAsAWUAQEVrjHUSLNrnuoP5LSw9XWxeLcKOKP6cOAaC2f0bBeY3L6oA3LRbXbSxWKY1vSVXPGrTbn9htzuHFEcdd+GYajaDKLN9WqS6TwEyfxXD2fG422LqilQebU6XuOI4m/wAAR/BbdHDRU6avD62wxZg3N57UxEDFYQQ27zn/AITkg4f0fhGU46EOABA1yXkd05LLgjOBw5302n8gsr4Amyx4Z3SYSi8wS5jTMRsQZTkVjw1mOGoGQ91gIm5v+OayECDZY6OqQ8jUJ1yDq3/koMoyCDJQAQLIANwQUbViB/Wnfa+wNlsztWQAXsFjEe0uGqbMBnZcnkgy7UOYUgTkEIFrIKUORUIEZBCBBsEFQZBcHVKTPtvY3vICxe2YNtjiaAPF4RWwMkG1awx2C+9YfxAp7fgr/rVD5wiNram1a3tuEkxWYe66e24Wft/9B5INk5hCtb2ygYjXI3ik4/wVOLobqngv5INg5FFo0dIUqtXEM1CGURJqapjuMgQeC26dSnVZr03Nc2YkZIOYyCDJQAQLIANwQUbU2qAC9gkCcggu1DmFIE5BCBIsgpQ5FQgbghAg2QVBkEgbgoAIFkVRkpvQAbgkC9kQgR6oRY80ItmUIsboLAUAsOascSoBYXQALIBmgHEoBncoEX9Ui6RfMpF8ygEXHNCAhFxdCOJQCLHmrAUIsbqxxKCACAgAhALC6AcSgAZpF/VAM7lIvmUCLoRcc0i+ZQi4uqBCw1qD3uLm4mtTEfZbqx+YKzEcShFjdBr+yOi+LxB/FvJBgwQJxFc/8cfRbMcSoBYXUGhWwnv9HR6cucLvfXfqNHdNzwUwuiKFB1Ul1V5e4HWNV0mwzuvoAcSgGdygwDBYcW1Ce9xP8Vr16GDZWaH0qxcBrDU1yBs2Lfi+ZWpXwQfjmYpztYMb9gi1pId3jmgvsVAtDi/ENaRN8RUEfmuNHBYYMLKFWtqsIbArvIblbPuWdjhisKx4lrarJvFpCxYalTwLKODp65bDi0mNhGfxQcjgxH7RiPEXCngGsBjE4gySZNTetsixuvmYdh0XVZhm031adZxIe1v2MhB4XtwB3XDb9ldH7ViPiOSDCEi+JxB/4o+gWwBYXQDiUGuMG2818QR/mELCMNR9pLQMVrRBf0j43xM8fzW8Bncr5uIa3R+LfjKgL6byGwCAaZMCbmL2/JBs+wUZ/rMT/qanmT2CjI9/Ef6ip5lnpvZVaH03h7SLOaZBXIi4ugwew4cf2XH/ABPcfqVxdo/BuB1sNSdP7zQfqtkjiUIsboMTMHhqYinh6TB/usAXMU2QPdC5xxKgFhdAa0AWsgAugHEoBncoECUgSkXzKRfMoBAkIQEIuLoRxKDhVpMfSewyGuF4MKUMPTw9IUqQ1abfstmwG4cFkIsbqxxKCAWHNALIBYXQDiUADNIv6oBncpF8ygRdCLjmkXzKEXF1QICEWPNCOJQixuoLAUAsOascSoBYXQALJAugHEpGd1QMxmhmEMwhmCgsFQTAVvwUEwFAEoJugmNiCb5IF5S8peUvOxAMyEIKGZCGeCAZhWCoZgq34IIJgIJQTAQTGxAE3S8oJvkl5QLyhmQl52IZkKgQUMwhnghmCoLBUEwFb8FBMBAEoJugmNiCb5IF5SDOxLyl52IIGhoa1oAAsABko+mHOY45sMj6fxXIzIQzwQDMZqwVDMFW/BBBMBBKCYCCY2IAm6hbrSHQRuIVE3yS8oEGUMyEvOxDMhUCChmEM8EMwVBYKgmArfgoJgIAlBN0ExsQTfJAvKXlLyl52IBmQhBQzIQzwQDMKwVDMFW/BBBMBBKCYCCY2IAm6XlBN8kvKBeUMyEvOxDMhUCChmEM8EMwVBYKgmArfgoJgIAlIN0ExsS91QJMZITbIpNvRCbHkgX3IDYWKs9/wUBsOSgA8EBzsgNvRAc+SBtySb5JN/RJv6IBNxZDO5CbjkhPf8EAm2RS+5CbHkrPf8EEBsLFAeCA2HJAbeioA52TbkgOfJJv6IE3yQm4sk39EJuOSAZ3ITbIoT3/AAQmx5IF9yA2FirPf8FAbDkoAPBAc7IDb0QHPkgbckm+STf0Sb+iATcWQzuQm45IT3/BAJtkUvuQmx5Kz3/BBAbCxQHggNhyQG3oqAOdk25IDnySb+iBN8kJuLJN/RCbjkgGdyE2yKE9/wAEJseSBfcgNhYqz3/BQGw5KADwQHOyA29EBz5IG3JJvkk39Em/ogE3FkM7kJuOSE9/wQCbZFL7kJseSs9/wQQGwsUB4IDYckBt6KgDnZNuSA58km/ogTfJCbiyTf0Qm45IBnchNsihPf8ABCbHkgX3IDYWKs9/wUBsOSgA8Ek3sgNvRJzVAkRmhIjNDkqcioJI3oCIF1UGQRUBG9ARe6oyQbURJE5pInNXam1BCRIuhI3qnMIUEJEZpI3qnIoggIgXQEb1RkEGSCAi90kTmqNqbUEkTmhIkXV2ocwghI3oSIzVKHIoJI3oCIF1UGQRUBG9ARe6oyQbURJE5pInNXam1BCRIuhI3qnMIUEJEZpI3qnIoggIgXQEb1RkEGSCAi90kTmqNqbUEkTmhIkXV2ocwghI3oSIzVKHIoJI3oCIF1UGQRUBG9ARe6oyQbURJE5pInNXam1BCRIuhI3qnMIUEJEZpI3qnIoggIgXQEb1RkEGSCAi90kTmqNqbUEkTmhIkXV2ocwghI3oSIzVKHIoJI3oCIF1UGQRUBG9JF7qjJTeiBAjJCBBskCPVCLHmgsDcFABAsrAUAsOaAANwQAXsEAsgGaBAnIJAnIJF/VIugECRZCBuCEXHNCAgECDZWBuChFjzVgIIAIFkAG4IBYc0AsgAC9gkCcggGaRf1QIE5BCBIskXQi45oBA3BCBBshAQix5oLA3BQAQLKwFALDmgADcEAF7BALIBmgQJyCQJyCRf1SLoBAkWQgbghFxzQgIBAg2VgbgoRY81YCCACBZABuCAWHNALIAAvYJAnIIBmkX9UCBOQQgSLJF0IuOaAQNwQgQbIQEIseaCwNwUAECysBQCw5oAA3BABewQCyAZoECcgkCcgkX9Ui6AQJFkIG4IRcc0ICAQINlYG4KEWPNWAggAgWQAbggFhzQCyAAL2CQJyCAZpF/VAgTkEIEiyRdCLjmgEDcEIEGyEBCLHmgsDcFABAsrAUAsOaAANwSBeyAWSBdAItmUIsbrzDrG0x93wPyP8ydY2mD/d8D8j/Mg9QjiVALC68w6x9Mfd8D8j/MnWNpiP2fA/I/zIr08C2ZQDO5XmHWNpj7vgfkf5k6xtMfd8D8j/MiPT4vmkXzK8w6xtMfd8D8j/MnWNpif2fA/I/zIPTyLi6EcSvMOsbTH3fA/I/zIf8A+jaYP93wPyP8yD08ixurHEry/rG0wf7vgfkf5k6xtMfd8D8j/Mg9PAsLoBxK8w6xtMR+z4H5H+ZOsbTH3fA/I/zIPTwM7lIvmV5h1jaY+74H5H+ZOsbTH3fA/I/zIPT4vmUIuLrzDrG0xP7Pgfkf5k6xtMfd8D8j/Mg9PI4lCLG68w6xtMfd8D8j/MnWNpg/3fA/I/zIPUI4lQCwuvMOsbTH3fA/I/zJ1jaYj9nwPyP8yK9PA4lAM7leYdY2mPu+B+R/mTrG0x93wPyP8yI9Pi+ZSL5leYdY2mPu+B+R/mTrG0xP7Pgfkf5kHp5FxdCOJXmHWNpj7vgfkf5k6xtMfd8D8j/Mg9PIsbqxxK8v6xtMH+74H5H+ZOsfTH3fA/I/zIr08CwugFsyvMOsbTEfs+B+R/mTrG0x93wPyP8AMiPTwM7lIvmvMOsbTH3fA/I/zJ1jaY+74Lw3eZB6fF8yhFxdeYdY2mOwwXhu8ydY2mOwwXhu8yD08jiUIsbrzDrG0x2GC8N3mTrG0x2GC8N3mQeoRxKgFhdeYdY2mPu+C8N3mTrG0x2GC8N3mRXp4FsygGdyvMOsbTHYYLw3eZOsbTHYYLw3eZEenxfNIvmV5h1jaY7DBeG7zJ1jaY7DBeG7zIPTyLi6EcSvMOsbTHYYLw3eZOsbTHYYLw3eZB6eRY3VjiV5f1jaY7DBeG7zJ1jaY+74Lw3eZFengWF0AtmV5h1jaY7DBeG7zJ1jaY7DBeG7zIj08DO5SL5rzDrG0x2GC8N3mTrG0x2GC8N3mQenxfMoRcXXmHWNpjsMF4bvMnWNpjsMF4bvMg9PI4lCLG68w6xtMdhgvDd5k6xtMdhgvDd5kHqEcSoBYXXmHWNpj7vgvDd5k6xtMdhgvDd5kV6eBbMpGd15h1jaY7DBeG7zKdYumOwwXhu8yI6iiIiiIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIg/9k=\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/ufPZOgrdHK4\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f852a19f240>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('ufPZOgrdHK4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On voit (si on a regardé la vidéo), on sait (sinon) que suivant le point de départ, la méthode converge ou diverge.\n",
    "\n",
    "\n",
    "Regardons comment on peut résoudre notre système non linéaire avec la méthode du point fixe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### La méthode du point fixe pour résoudre $A({\\bf x}) \\, {\\bf x} = {\\bf b}$\n",
    "\n",
    "La première chose est de définir la fonction $\\bf g$ telle que la solution de \n",
    "${\\bf g}({\\bf x}) = {\\bf x}$ soit aussi la solution du système matriciel non linéaire.\n",
    "Cela implique que ${\\bf g}$ est :\n",
    "\n",
    "$$\n",
    "{\\bf g}({\\bf x}) = A^{-1}({\\bf x}) \\, {\\bf b}\n",
    "$$\n",
    "\n",
    "Bien sûr on sait qu'inverser A est trop cher aussi on préfère écrire notre\n",
    "algorithme itératif ${\\bf x^{k+1}} = {\\bf g}({\\bf x^k})$  sous forme de problème matriciel linéaire à résoudre :\n",
    "\n",
    "$$\n",
    "A({\\bf x}^k) \\, {\\bf x}^{k+1} = {\\bf b}\n",
    "$$\n",
    "\n",
    "Si on connait ${\\bf x}^k$ alors on peut évaluer $A({\\bf x}^k)$ et donc il\n",
    "s'agit d'une matrice de réels. On voit qu'alors le système est bien linéaire\n",
    "et qu'il permet de trouver ${\\bf x}^{k+1}$.\n",
    "\n",
    "Cet algorithme itératif peut donc marcher mais on sent bien que cela va\n",
    "dépendre du type de la matrice A (on a l'habitude) et de la valeur initiale\n",
    "${\\bf x}^0$ comme souligné dans la vidéo. De plus on voit qu'il va être\n",
    "très cher puisqu'on a un système matriciel à résoudre à chaque itération !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test\n",
    "\n",
    "Regardons comment marche cet algorithme pour résoudre\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "x_0 - 2 x_1 &  x_1 \\\\\n",
    "x_0 & 2 x_0 + x_1 \\\\\n",
    "\\end{bmatrix}\n",
    "\\;\n",
    "\\begin{bmatrix}\n",
    "x_0 \\\\\n",
    "x_1 \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 \\\\\n",
    "9 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Il y a deux solutions qui sont [1, 2] et [2, 1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^1 =  [1.5 2.5]  erreur² =  0.48999999999999977\n",
      "x^2 =  [0.73913043 1.43478261]  erreur² =  0.18534274820344412\n",
      "x^3 =  [1.37617066 2.74037461]  erreur² =  0.6387945198307682\n",
      "x^4 =  [0.72846505 1.45602065]  erreur² =  0.18097695808599654\n",
      "x^5 =  [1.37323934 2.74623359]  erreur² =  0.6430523663786714\n",
      "x^6 =  [0.72824232 1.45646606]  erreur² =  0.1808855461166179\n",
      "x^7 =  [1.37317932 2.74635363]  erreur² =  0.6431399506048754\n",
      "x^8 =  [0.72823777 1.45647516]  erreur² =  0.1808836795920618\n",
      "x^9 =  [1.37317809 2.74635608]  erreur² =  0.6431417383272344\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as lin\n",
    "\n",
    "def A(x):\n",
    "    return np.array([[x[0] - 2*x[1], x[1]] ,\n",
    "                     [x[0] , x[1] + 2*x[0]]]) / 10 \n",
    "\n",
    "b = np.array([1, 9]) / 10         # avec normalisation grossière\n",
    "\n",
    "x = np.array([1, 1])\n",
    "for i in range(1, 10):\n",
    "    x = lin.solve(A(x), b)\n",
    "    print(f'x^{i} = ', x, end='  ')\n",
    "    print(f'erreur² = ', np.square(A(x) @ x - b).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La solution oscille sans converger. On peut tester d'autre ${\\bf x}^0$ \n",
    "mais on obtient souvent le même résultat.\n",
    "\n",
    "La méthode du point fixe de base a un petit rayon de convergence à savoir il\n",
    "faut que la valeur initial soit proche de la solution pour converger.\n",
    "\n",
    "### Appliquons l'inertie\n",
    "\n",
    "Ou la méthode de surrelaxation qui consite à avancer moins vite vers le prochain\n",
    "point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^1 =  [1.99601592 1.00199403]  erreur² =  2.845730924209918e-06\n",
      "x^2 =  [1.98435526 1.00784963]  erreur² =  4.338186795289511e-05\n",
      "x^3 =  [1.94167764 1.02955991]  erreur² =  0.00057735040271093\n",
      "x^4 =  [1.81575881 1.09680461]  erreur² =  0.005006527863139465\n",
      "x^5 =  [1.58849624 1.23874527]  erreur² =  0.017837680043912996\n",
      "x^6 =  [1.3654867  1.43296028]  erreur² =  0.023567622207356913\n",
      "x^7 =  [1.21494033 1.62266463]  erreur² =  0.015938576382362332\n",
      "x^8 =  [1.12471347 1.76667765]  erreur² =  0.007550226984782101\n",
      "x^9 =  [1.07194609 1.86111664]  erreur² =  0.0030000984268709436\n"
     ]
    }
   ],
   "source": [
    "µ = 0.5  # j'avance de moitié vers le prochain x\n",
    "\n",
    "x = np.array([1.998, 1])\n",
    "for i in range(1, 10):\n",
    "    x_old = x\n",
    "    x = lin.solve(A(x), b)\n",
    "    x = µ * x + (1-µ) * x_old\n",
    "    print(f'x^{i} = ', x, end='  ')\n",
    "    print(f'erreur² = ', np.square(A(x) @ x - b).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et là la convergence est rapide. L'inertie augmente significativement le rayon\n",
    "de convergence de la méthode du point fixe (plus µ est petit et plus le rayon\n",
    "de convergence est grand mais moins on avance vite).\n",
    "\n",
    "Bien sûr l'algorithme converge vers l'une des solutions et y reste. Il ne peut\n",
    "pas trouver toutes les solutions. Si on veut trouver au autre solution il y\n",
    "faut choisir un autre point initial et avoir de la chance d'être dans le \n",
    "rayon de convergence de cette autre solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## La méthode de Newton-Raphson\n",
    "\n",
    "Il s'agit d'appliquer la méthode de Newton pour une fonction vectorielle.\n",
    "\n",
    "Voici une vidéo qui rappelle ce qu'est la méthode de Newton pour une fonction\n",
    "en 1D. Vous pouvez vous arrêter à la 4e minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAEAAgMBAQAAAAAAAAAAAAAAAQQCAwUGB//EAEsQAAEDAwEEBwUGAgYIBgMBAAEAAhEDBCExBRJBURMXYXGR0dIiU4GT4RQyVJKhsTPBFSM0QlLwBkRic4Kis8IkQ2NysvElg5QH/8QAFwEBAQEBAAAAAAAAAAAAAAAAAAECA//EAB0RAQEBAAICAwAAAAAAAAAAAAABAgMRMVESISL/2gAMAwEAAhEDEQA/APn6IiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6U6udse/svmO9KDyCL1/Vztj39l8x3pTq52x7+y+Y70oPIIvX9XO2Pf2XzHelOrnbHv7L5jvSg8gi9f1c7Y9/ZfMd6VHVztj39l8x3pQfUCcITgoSI1QkRqqiZUA4CSOaAiBlQAcIDqgI5oCM5VCcpOUkTqkidUAnIQlCRIyhI5oBOCplQSI1SRzUAHAQHCAiBlARzVGutXbQaCQ5xc7da1oy48goNwGUDVqsdTgZackdmOK1yH7SAJ/hU5jtcYn/lPitW1a3RUqRA3iarQG8zkj9QEFulVbWY2o2YI4hZk5HktVuwUaFOnv7262C48e1bd4GMhAJQnBVTaVZ9JlItO7SdUArVAcsbBz4wOwGVlZEm2Li5xa4ks3zJ3eHnnmoLUqAcBN4cwudtoPdZ0RRfuv+0UoOo++NRIkIOiCgOq51Gba9uqjjUrvNKnMcTvPwBoAttS9Nu+j9pbTpMqOLS41PumJHLWCgtte1xMGYwY4KZyudsf2mVa4BayoZE43jJJd8Z/RdGROqATkISqtm+o81qlSrvNdVcGNj7oBLY/SfirRI5qgTgqZVKveAbRtram9p3w9zxqYaB/NwV2RzUEA4CA4VGwvqte4r0LiiKTmH2IdO83H65HiFYs6prWlKq8iXtDvFBuB1ScoCM5Wm1rGtSL3wDvuAjkHED9AqN05QnISROqEiRlAJQnBQkc0JEaoJlQDgJI5oCIGVABwgOqAjmgIzlUJyk5SROqSJ1QCchCUJEjKEjmgE4KmVBIjVJHNQAcBAcICIGUBHNUAdUnKAjOUkTqgTlCchJE6oSJGUAlc/b99V2fsitc0A3pGFsbwkZIC6BI5rRfUWXNq6k9jajTEtcAQcqXw1iyal14eX2T/AKS7QvXVekFEBjZG6w65PPsVraG3L+0tKj+jY17MDebIP3IOvJy6FHZVvSJNO0osLhBLWN0lbqtoK9Isq0GPa7JDmtPLyCzJfj07b5OO8s1J+fTjf6Nf6Q3u1Nout7kUtwUi/wBhpBmR29q9ROq52z9n0LSualO2pUju7u81oB/RdGRnKuZZPtnn3je++OdQOik6FQQI0QgQcKuKUGgSByCgAQMIqRog4qAByCADOAiJ4pxUQJ0CQJ0CCTqEKggSMIQOQQSdCiggQcKYHIIA0CDRQAIGEAHIIKtQmheOrFr3U3sDTutLi0gmMDOd4+Cxvbc31OhuiGio1zg8RLcgiO4q4AM4CQJ0QcW5o3lS2pW9vQ9m3qe2x5htVokBvaIgrU256DZuz7cXDKFQVNy4NKPYIY9xEHtbxXfgToqG0HUxdbOl7RFwTr/6b0GN+7f2ZSbdncbV3W1jpGJI+JEfFXbcvNpSNURU3BvDtjKOrUAM1aY/4gsX3Vq1pLq9EDteEGQtqDaxrCizpTq/dzy1Wq9pvqNt2sExWY53YBn+Sn7fY/i7f5jVH9IWIAm6ofnCK13TLhv2mpbt3nPogMgwQ4E+eO5aranW/wDCCpSdu06hO890nLXfpkATlWBtCx/FUPzhZi8tTkVmEdiIq2NnWtX03OIc8tc2oQcR/d8AI+Ko06N+XUHV6dxTYxwFw0PDjXO6ZcIOGzHI9mF2PtdtP8Vifa7af4rUHM+z3IsGW7GVaLBcPLo9p3Ru3nDjzI8OIXWa8mhvBjwYwHarA3dtj+tahu7b3rEHEt6V5Wey5t7ZtKpTtHUi+phxqOgkjnkfr8F2tnNfTs2MqNqtc3BNUgk8Z1OPiodf2QkG5og8i4BSL+ydhtxSceTXSiq4o1emuKzGnfbWloON9ppsBA8PEKo+jdP2ds4MfVptZTDatNrZdvACNDwIPYutRr0qx3WBxIEyWEDxIW0AckRxqL726rWtKuyvTDWtdUe0bgL91wcJB0kt7+C01P8AwlKlUt21h01aoyp0AkudLonlkarvgDOAkCdEHLtKu0KTzVviXUw4s3adLu9qNYmR3QV02uD2tcAQDnIg+CmBOiECRhBJQ6FQQOQQgQcIJQaBIHIKABAwipGiDioAHIIAM4CIninFRAnQJAnQIJOoQqCBIwhA5BBJ0KKCBBwpgcggDQINFAAgYQAcggkcU4qABnASBOgQTxQ6hRAnQIQJGEElDoVBA5BCBBwglBoEgcgoAEDCKkaKOaADkEgZwiECPqhGD5oRjUoRg5VEwFAGB5qY7SoAwMqABhANUA7SgGuSgRn6pGVpua4t2bxBd/xNEDicnRYXF0Lc25qHdbVduknh7JP8kFkjI80IC1069Gs4tpVmPc37wa4EjvSrVFNzWkPJd/hbMINhGD5qYCg6HK1V63RMYWjfc9wa0AxP+RJQbQMDzQDCps2hRqBjqRL6ZjedukFs/dMEaK4BjUoAGqRn6rFzm02Oe94a1okk4ACinUbUc8MdO4d098T/ADQS+myo0seN5pGQeK1CytWwG21Jo5BgC3xnUrB1Wm17WGoA5zt0CckxMeAQQLeiz7tJje5qyNNgH3QsiO0oRjVA3W8kAEBCQMF3atFtX6dgcYaHZYN7JbwMcEG8AQgAytNWo6lUpA5pvduk8jw8vitwGuSgQJSBKxe9tNzQ50bx3R2lASahEOAHExBQZECQhAWm8qm3tK1ZuTTY5wB4wFD7hvtspObUrUy0PYDlsnj8MqjeQIUwFTdegVnMI9jLWuyS5w1wBp281mLqbalW3T/WOa0tn7pJj9FBYAEBAMLDpG9I2lvHfLd6I4LGncUqlzUt2PmpTAc8D+7Mx+xQbQNUjP1WJc1rt0vAJ0BOq0W95SrWQu+kAp7m+6CDu4kjHJUWYyhGR5rDf/rGtAdlu9MYUl7OkDN8b/8AhkSoMiAhGPqhHaVVubynT2bWu2PBaykagnsEoLcBQBgea5exr24uKNT7TLixxAIHtfeIggcYAPxVyle21QsaK7Q587rXEAmJ0HwKCwBhANViHNg+2MGNQpEEuAdkaoJjKRlaBcs+3utd4B4pipBIkgkjT4LZ01LpCzpW7wO6RImYmPBBmRkeaEBa6tejSE1KzWid3JGvLvWizvPtFB1Wru0x0m6Bvd2D2yYQWyMHzUwFU2lfUdnWdSvWd91pLWCJeRwC2vrtFN5Y5r3sGWbwBk6A8kG0DA80AwsaTi+k1zsTwBBVG42hu3Vk2jVY6hVNR1R8yA1o58Mwg6AGqQJ+qqt2hauc1rK4cX1eiEf493ejwCrbYqB/RWjakOquO9Do3WgEknlGD/8AaDpwJQjI81roVqVw3eo1Q9owSOag3FLp3US+HtAJntnyKDaQEIwfNV7m6ZbVGCs5rKbwfbc4CCMx4T4LdvNdT3w6WkSDwhBnAUAYHmsKNWncUhUo1A9hmHNMgwYWY4CcwgAYSBlRIgy7TVGkObvNdIOQRxVEmY1QzCGYQzBUEwVAmApz2KBMBAEoJygmOCCc6IK95ZsvGNZUcQGuDwW6yJhZuo7/AEJqOl1M7wIESYj+ZW3MpmeCo51rZssb8NY+o/paRgvdMBrtP+db7izNW7o3IqAOpYALZBnWf5clvfuiowu3d7IbJz8FmZ7FBUq2b3msOkinVbkHeJ/UxHwWmytC2i63NYl1tXLmvDQMkScaR7ZCv1KjabSXlrcE5KilSZSaRSa1oJLjHEnUoKtrs80qVxTq1BUFcuJO7BzONe3CutBDYnRBMBBMcEEQSCMZ7FRtqdXZ7rezpUzVoHemrEFvETGP2V8TnRMygZnVVq1oHXlO6Y4NrMG5JGC3iP8APmt1aqKFJ9WphjGlzj2BZ5MFBzDtZ/RB5sbhpl3subGgJB+IB7lbvLZl9bBjyN2Q4eyHfvK31HBjC55a1oySTgBQ1wfTDmRukY7lRyrixo295QrdM/e9lgpbrd2Puk6SMO0mM6KwdnVGkVaVyRXbRNJj3MBjkT2rK3u7TaFRsBpqsLy1rhkBrt0nulXRMBQUqVG7q2Ip3homu1wO9TndJBBmOGRoronK03Nyy1pB9SYc8MEZMkwtoeN8slu9rE5hBoubQXL2l9R43ILQIwZ17+HxWdKi6nVqO6QlrzvBpH3ea25ngsKVVtZu/TIc2SJ7jBVGF7SdXs61JpzUpuaPiCqVJjKO2S1m+HVaZqVA4yASRp+viuhXqijSfVf91gLj8AtNoftFNl1UosZWc3dJBmBOk8uKDT/Rp6ZlU1s0nvcyGD+8ZMzP8lsNrUbbNphwcRX6T4dJvfsrZmOCxfUFMAvIbJDRJ1JUGllqadwa1KpuioZqNInexw5cOzsWFva1aN3UqNNFtJ+rW086k6z2q2JgLF1QU90OIG87dHaVRUudmU7m56Zz3Ndu7hjl2ctSJWYsKf8ARzbFxJpNpimYxIAhWS4MBLi0CdSYUCo0ueA5vsfezpxyoNNrbPt3Pb0gdTJkAt9qeMmc+CrXdtc07p95bblSq4NptY4QGiRJmddfAcleNVorspE+29pcB2CJz8Qsi4B7WktDjoJ1Qa306lS26N1Xce5oDn0xoeMSuBtLYdGy2RUNO6uZZS3G7zgd4+yMkj/Zb3CYhegqVmsD8hxYAXNbk9mFgx9O8t5dSBbvEFlQAwQf5EIORsOw6a1F06q7fc+sxxIglpqGdIgktmforTdgWm7DzVe0kncLobmeAjnCt3Nzb7NoMdWLaVIvDAdACT5qKe0Ld9AVTUaymXFrC5334xjmgqN/0c2a1z3Nolpe4OdD3CSDI481k7Ydud8Mq1qbXO3oY+I9nd1108IxC22W0xd3ApsouaDT6Qlxgtk4BHCRlWq1dluzfquDW7wbJ5kgD9SgwoW76dd9R9XpC5rQSWwcT5rmUdmU7nat7dPqEvD+jZ7EFhhhmeOgXazK03NzTtAx1YwKlRtMRmXOIAQVamyKTmuDatRheXbzhBMO11H66rSzYFCjTpso1agDKvS7rstLyIJMQTPfC6xmQtBu6UPh09HUbScADhxiB/zDxQar7Z/2x7HGu9ga0gtaBnII17QO9Yu2VTLLlgrVWtuTLwCNeMGJHjjgrxmCoa8PndIO6YOdCg10qG5atoEy0N3cCMfBc6j/AKP2tO1Fuatw5jWlrSahBaDrkRyGvILrCYCCY4IKVHZVrRDdxhljg6S4kkzMntnU8VuubK3u2Pp3FJj2vEOkZOCNfiVvE50TMoKlhs+ls5j6dsXBj3F264zB4mdTPaVWGxKVSrUuLqo99xXjpC1xDRA0aOAXUzPBDMhUc5mxbdlOpTLnva8yN8726YA/7Rr28yrdW3NS0dQc8wWwTAz8NFuM9iGYKg5NTYZrPe6reVpJkGn7B1nMYPIdhPetjtjtfV6Z1zWLzBmeIdvfATGOwLpZ7EEwEHDo7LdcO2gx1aqRVcAX1GRJDiSMRLcga6Y4LsUKbqVBtMlvsiPZbA+AWwTHBM5VAkxohONCk4+iE4PkgZ5IDgYKme/wUA4HkoAPYgOuEBx9EB18kDjok50Sc/RJz9EHK2xa16znPoMLnmkWUzAIa6ZzORwyOSzB2jVFR9MsFOow9Fv+y5hyRIjtA+ErpE5HkhPeg41a12m+1DXP36jBLXNfumd0tzHwd4rYLfaJe11vWNCkGN/q6kPJdPtTrw7V1ScFTPeg0s6WabnNAJbDwHYB7Oaxu21H2r202uLjwaYJE5APAxK3g4CA4+iDiWz9turucKLKVCm4t6F5BLgNId2gDXiTyV22dtQVwLhlu6lMFzCQ6Ocf50Vs1qbajabngPfO606mNYWc5+iDlV6G1q9KpQc+2FN9N7CTMmZg4iMHTs1XRos6KkxgDsc3F36lbJyhOR5IKNxbVa206NRzGm3Ywgje1JIOneB+vYte2aF9W+zvs6pY2lUD6jGmHVIIxMxGszzXSJ70JwfJB5D/AEa+21C+tSFIOgbzXZI9recJ7STw/urqWlPbtFrA7oDvOl01C+PZzkxxHD/F2LpWNhbbPpubbUy3e+8TklWQcBBxbu0v7qvvMbugVWVAKzgWANyAAOM6k/qrdF20Te71alSFD2hDXSRMQZjSAfieyVfBxxQHVBovWPqWtRlNm84j7pMb2cie3RUtl219bPaK/R9FuR0bHYYd4kkYzqB8F1Jz9EnKCttGhVurGtQov6N9VhYHnhI1VPZdrtG2qVhdPYaRAbSDXk7sExiBGCB/w9q6s5HkhPeg4zbTbFMtbRrUWU204G851T2szM5M47ojit9a2v69xaOqOpdHRcHvAJ9twkaR2g9nbqukTg+Sme9BgwuI9psZxBmQtF22o425psktqz2fddqrIOAgOPoqOdZ1rqvWdTvrYUw5oe1jt0wQRyJnUZ7NAtFK2vDVoipTqCgGRVYHgF7yTLjBgt48+xdGnbht5UuC9x3wAGnRukx3wFvnP0UFGjTqUrq1Y4SGUqgkZ/vN3RPcs6VrUpXL6hqGoHgw5+XM7Bwjy4rfVuKVKTUfG62SOMdyxq3VCjVp0qlVrX1PuNOp4fzHigr7Oo1aNCqyqx4qOqOcahcHF0nB+AjgFYtminbNY1j2hoj2tT2/HVbie9CcHyVGi6oms639kFtOrvme4/zIXNfs6/32Flanu06JaySR7RkEkDWZHcV2p/zCgHAUHJsrG8s6lc03scyoDugnIIaACTGcg9gwqzrHbFxVt23T6Jt6VQVHtFSS4gyB90aGMzw0XfBx9EB18kFbZ/2oWrRenerid4gAA54Qqe0LC7vL6k8PaKVF++wFx13cYjg4z4aLqzlJz9EGm1FVtvSbWk1GgBxLpk85wqbrK4Fd9VjmxVrb72uMxAhpHg2R2ePSnI8kJ70HHpWu1uhFG5rUqlMxvODyHCCDyHb+nesH7P2nSY77Fcim4h7oe8uAeSY1BxBHdHFdsnBUz3oNNqa/QxcAdICRLRAI4HXktoPYgOAoFRpcWhw3hkjigkHXCcdEB18knP0VCc6ITkYSc/RCcjyQDPJCcaFCe/wQnB8kDPJAcDBUz3+CgHA8lAB7Ek5wgOPok6qgSI1QkRqh0UnQqCJHNARAypQaBFQCOaAjOVI0QcUREidUkTqp4pxQQSJGUJHNSdQhQQSI1XPr3LhtelbtrOAc0PLd5u7GcaTJ5TwPJdE6FYuYx33mg6HI5aIqtf1HU6VF1MPdFVstY6CQsdmbRp7RoOqMY+nuugteIOgP8/GVcGgWqnSc27q1SRuPa0AcZBdP7hEUqhbU27TD3sa2lTDmy7LnHeEAd38lVp3lW7ruthcsArsLiHMyyZG60yMiNNdSusKAF4a4DRLN12MnSM+K3cUFVt8wUbZ7wR0zRGRgmMf55LVcXzzcto2+7LXgO3hO9nIHcJJKvwCRhY9GwVC8MaHnBdGSgyJHNCRGqkodCgiRzQEQMqUGgRUAjmgIzlSNEHFERInVJE6qeKcUEEiRlCRzUnUIUEEiNUkc1J0KIIBEDKAjmpGgQaIIBGcpInVARJHFTxQVbu2Zc1KG+ym4U375LhkRkR8YWq6taVW9pvcTvvaW72PZjIj45V/itFX+12//ABfsitxI5oSI1UlDoUREjmgIgZUoNAioBHNARnKkaIOKIiROqSJ1U8U4oIJEjKEjmpOoUOIAk4AQVri8bTpXLmtLnURpMbxIkAeI8VXbtCv/AEfa3j7cjpGAvpMcCQT90SY4n9ljZTeOdWP8LpC//wBzv7vgAPj3LoPpMeKYIgU3bwA00RWqjcVOkNKrTjcYC6qMMJ5CVyqVyKNjV2iKlN1Wu5wYQ4ENEk5OhMAdmAO1d0aLXb0W0aDaYyG6IKttdGqbZ7KjnUqjXN9sAOJHE8tDjtV6ROqjdaXAkAkHB5LLiiIkTqhIkZU8UOoQQSOaEiNVJQ6FBEjmgIgZUoNAioBHNJGcqRoo5ogQI0QgQcJAj6oRg+aCYHIKABAwpgKAMDzQAByCADOAgGEA1QIE6BIE6BIz9UjKAQJGEIHIIRkeaEBBruKjaFu+q5stYJMclNCo2vSFQMLWuy2eI4FZPaHNIcJBEEHiqn2KrTvGVLe56K3a3dNvueyc6jlr+yDKvfUrZ7umAbSpsBdUPMmA0DUk+Swp7W2dUpNqNu6G65+4CXR7XLPHIWy42faXebmi2rgCHZGMjCx/oqwIINpSILi4y2ZJJJ/+R8UFijVpVml1JzXtmJCzgTotVC1oW+90FFlOddwRKUmPFas559kkbgnQR5yqIbcUX3TqDc1G/eG7pp5/vyW4gYwFXr2bKhrOpudSrVafRmqw+0Bw8JPiqhtdq/Z6TRfUukawte7o8uMYdrGvZxUFq9t3VQx1MuD2kQA6BkiT2wJVkgRosXM32N3t5pkHDlkRjj4oJgcgoAEDCmAoAwPNAAHIIAM4CAYQDVAgToEgToEjP1SMoBAkYQgcghGR5oQEAgQcKYHIKCMHzUwEEACBhQYDSd2Y4DipAwPNAMfVBzLJvTbTuLms0Uy13R02Ojejdac/vHaVhQv6f294eTvu5vAYxoJER/iwO3PIYt2tjToXNzWFKk01am8HNGT7ImT3yt9a3pXFN9KqwOY8Q4Hig1u/q7+n/grMLdcbwyI+E+Cxu3Np3VmXEDeqFo7TunyWdKzoUejbTaQ2mSWDeJDcRj4HRY3gHSWrsy2sOPNpH80G19Wm2ZDsEAwwnX4IKjXNqQx3sSDLCJ7p1WwgIQIQanVg2kH9DUOYIDcjthZNewv6ODvbu990x46fBbICgAQEGplZjnPaGP8AYmZpkA92M/BGVd9j3Ci+Ro0iCfFbQAgAygxYS7dJpFsiSDHs9hygM1CCyANDjKyjP1SBKCubh4eGmzrAT94bpH7qjtuvWextjbUXipcPaw1MQ1h+87wx3rp16jKNM1HzDeWp7AqlhS6RzruqAalQw0icN5Z7Zz8eKC3Toso0W06bYa0QAtkDkFBGPqpgIIAEDCADkEAwPNAMIAAzgJAnQIBqkZ+qBAnQIQJGEjKEZHmgEDkEIEHCEBCMHzQTA5BQAIGFMBQBgeaAAOQSBnCAYSBlAIxqUIwcoZjVDMIJjtKgDAypgqBMBAA7SgGuSglBOUCM6lIzqUzKZlAIyMoR2lDMhCCgEYOVVuLs0r+1tmt3jW3iTyAjzVozB0XNfSuP6epPNVm50L49nMbzZH7Z/QILAq13VmUoDSXOceMMGB8Tj9eSzuX1G2rqtB4JZ7cRO8Bkgd603Vg65c/+u3GVGNY4Ae0IJIIM415cktad9Se+jV+zvt4im5jS0tEaET+youMO8JBwchTGdStdtSNG3p0t6dxobPOAtdvQfTr1nuM7xke2447iYHwQWIzqVz7urUqbSt7ai97SJqPMECAW+OCR8exdDMqneWxq3dB5FQta14O48t1gwYPYoMTeEXYBfLXP6OnTAEuIPtOnkM47Fuq1S2pRe14dRqeyYjU6Gf0+IVd1vvvNR2z6ZJM+08TPONJWqjY/Z6HRUbNzKYIcGivIEGRAOBog6Jr0em6HpmdLrubw3vBaLO5qV3DfDQx9MVGRqAef6LLpbv8Au2rQe2pH7BLKgaLXTRpUpj7jy7HLIEdyCyB2lANclBKCcoEZ1KRnUpmUzKARkZQjtKGZCEFAIwcqY7SoMwpgoIAwMoB2lBMBBKABrkpGdSgnKZlAjOpVa8Ht2wnJrD9if5KzmVWu/wCNaZH8b/scgskdpQjByhnmhmFRMdpUAYGVMFQJgKAB2lANclBKCcoEZ1KRnUpmdVSuarqznUqZPRjD3N1cf8Le3meH7Brc77feBgcehp5w6Cdfa8RA+J5K5csa6m3efuRUYQe3eGPjp8VFvb9AyJBJMkxHcO4DCw2kSLZpB/8AOpf9RqCyRjVTHaVBmNVMFBAGBlAO0oJgIJQANclIzqUE5TMoEZ1KEZGUzKGZCoEdpQjByhBQzCgmO0qAMDKmCoEwEADtKRrlBKQcqgZhDMFCTGiE40KCc9igTATPJAcDBUATHBBOdEB7EB1wqGZTM8E46JOdEAzIQz2ITkYQzyQDMLEsBqNeQN5oIBngdf2C13NwaPRhtMvdUdugAgcCf5LB1W9J9i2ox/t1iD+jSoLImAgmFVFTaEf2a2+e70KWPvv71C3HdWcf+1BZE5TMqrv38mLe2I7a7h/2Kd++4ULeeXTH0oLOZQzIVYVL6fatqHwrn0qTcXDQN6ze4/8Apvaf3IQa76hcVbi3fSdFNhmo3fI3hII+OJ8RxU0roX9NtWwuaLqbXFr/AGS6Y4aiCsjeQ327e4b2dHvf/GVqp1bOgD0NrUYc/dtnDjPLmUGt1YWF8595XDKVw8MoNLiRvH4cceKs0doW1d7WUarKjt4sO6ZgiZnloVTo7Uo3lMitZ1HNNUsYOjLg7EzkY4ra+vRtqAeywrRSO+AymBHAkZ5EoOgJjggnOirNvJaCLeuQf9hbKNwKtRzNx7XAB0OEY/yEG3MpmeCcdEnOioGZCGexCcjCGeSAZgqc9ignGhTPJQBMBBMcEBwMFAexUBOdEzKA64TjogZngq15P2iy/wB+f+m9WZzoq13/AB7MwcVv+xygsmexDMFDPJCcaFUTnsUCYCZ5IDgYKgCY4IJzogPYVWq1X1ajregSCP4lQf3ByH+1+2vKaIqvfXqOo0nbrW/xKgOnYO39llbNa4NLGMFJmKY3YI4E/FQ5gG7a0qTTSIipJ0b+8mf3VnQjCgGZCq7Tn7K3/fUv+o1WicjCr38fZTvNJAc0+DggsGYU57FBONEzyQBMBBMcEBwMFAexUBOdEzKA64TjogZnghmQk50QnIwgGexDMFDPJCcaFBOexQJgJnkgOBgqAJjgmcoD2JJzhUJx9EJwfJCRGqEiNUEz3+CgHA8kkc0BEDKgA4+iA6+SAjmgIzlAnP0Sc/RJE6pInVAJyPJCe/wQkSMoSOaCtefxbV0aVT+rHBRf0q1WnT+zuLXNfLoduy2DIU37g2lTdOlamPFwH81Zkc0HPqbXZR6fftrkNogy8U5DoE4j+at29yy4a8s3vZdB3mkHQHQ9hC2yN3gsWNYxznNABeZd28EGFxTdVpFrHFr5kEEj9lnSBYxrXEucAJPNZAjOUkTqqOb0lawuaj67a1wys8w9mejbkgFvxOR/JdHeBg/yUyJ1QkSMqDXVFUlpp1A0A+0CyZ/VTXD30XNpP6N5GHbsx8FmSOaEiNVRwrcNbZ0aLrt1KvTY9k9CQckSY5yJntVuzqXJuW0xVZVtm0wDvU3h5PEycFa2Clc7ZvW1DUD2tbTaRTIAaQCfaiNf5pQuL20M1rR76REuNKHOc8nPHQcO/hGYOjas6G2p0i7e6NobPcsGn/8AI1P9039ysLUNbcOdSaWU6rekLS3dh0mTHAnj3LJjgdpVeyk393eSoszn6JOfokidUkTqoBOR5IT3+CEiRlCRzQCcHyUz3+CgkRqkjmgA4HkgOPogIgZQEc1QB18knP0QEZykidVAnP0Ve7P9daf77/sco+2sNy2l0b91ziwVMQXASRrPA8Fpur20NW1i6omKs/xB/gcqL5Pf4ITg+S0G+tPxVD5gWX2u2cDu3FI9zwg3T3+CjeAbJ5KrX2nZW7ZqXVKZgNDwSTyAVCptChcVG061zbhuvR9KN1o5uP8AeP8Asj481BZqXdS7qChZGGnL68SAOzv5/vmMreoyhs8/ZWOcQ9zGB+pdvEZ+IJnkpoXuz6LIF7bFx+841Gy48yq1hcW7Xuq3V7aufvv6JjajfZaXE+Mf51kOlQpNoggSXOO850Zcea2Tn6KuNoWhdAuGE8gZWX2yhP3/APlKDcTkeS0Xx/8AB1e5Sbyhj2/+UrVd3Vu63c0ue6Yw1pnVBbJx9FM9/gq772gB95x/9rHH9gsft9GYDLj/APnfH7ILIOB5IDj6Kv8AbGQIpVz/APrI/dPtmDFvXP8AwjzQWAdfJJz9Fy7TadX+j693d0w1lMktIEb7eECT2d8rpU6gqMa+I3gDB1CDKc/RCcjySROqEiRlUCe/wQnB8kJHNCRGqCZ7/BQDgeSSOaAiBlQAcfRJ1QEc0kZyqB0UnQqCBGiECDhQSg0CQOQUACBhFSNEHFQAOQQAZwEQJAknAAUMe2o1r2ODmuEggyCqW2KvRWUBhear2Uw0CZl2R4So2ZdW1Uvs6NTpKlq1oqGMGZEz3tKDoHUIVWfe2jLsWrq9NtcgEMJgmVQu9qGzrMZUqUnuNTdfTp03EsG45w7zIHwQdaoxlRsPaHAEOAI4gyD4rJUtlVnXOy6FWo+nUqOZ7bmaTxWqx2nTubqvSMCKu5SECSAxricd/HmgtfZW7xcKtYE/+oY8FItsfxq351q2m6rS2bXq27g2oxhcDuzpwVenta3pteazj/GLGjdyPZDjI7Mygtm1fJ3Luuzu3T+4KgWtaf7fcHvbT9Kr3W07elWYwvLBIc924SN2Dx5+zpqrNrfWt4XC3qte5o9pmjm94OQgGhXn2bt//Exp/knQ3Eibs/CmFtrVGUKT6tTDKbS5xjgFRp7WpValuBTew1XBpY9vtAlm+BA0xM9yCy63rEf22sO5rPSoNrXj+33H5afpU311TsqIq1R7G9BPLBP8lqbtWxqNqFlUHo6gpOG6ZDid0DxQbBZvBLjeVy46mGAnwap+xtc0b9au7HvS39oVgwBJhcqltu0uWU3W9emz+sDHiqwh2TAgdp46Iro0aFOjJZvZ/wATy791sGpK4ljtH7ZttzS/cp0xUpsZOKhDyC7/AJf1711rmtTtqD6zx7LYHDU4H6ojdxTiqGyLtt3aMO9vvDGlx1GRIzxxCzv7+jYdEarSekJHsiSAAST8IQXDqEK5hu6r9itueic2o9o3QIn2jiMxx4rY24FpasFZ7Xig1ra75ktdjJ8ZKC+dCiq0r2hWrCgARW6MVCwty0HnyPYrUDkEAaBBooAEDCADkEEjinFQAM4CQJ0CDTTtaNOu6q0HecScuJAnWBwWb6VN7mFzQSx28O+CP5rOBOgWFbfDP6ljS8mPaMAdqCarqdJhfULWtGpOAFVe51yIo0GtZ72qz9m6n4x8VsZaNDxUrONaoDILhhvcOH79qsECDhBRt206Ny6m21qvfgOrOLDI15yBnSFeAEDCrts2NvXXQe7feN1wMQQNBp/me5WABAwitN3c07O2dXqg9Gz7xHATr8FlbvNWkHuZub2QJnHBZljXAbzQYMiQqL6V+y7cLb7MLTot1rS0hzH8+0diIv8AFTxVWgy6FVhrupkFh3wwQA7Efz/RWYE6BAOoWq6uadrTFSrIYXASOHLxOO8hbSBIwtdehSr09ytTa9sgwRhBkxxfSDnNLCRJaeCzVOjTvKdzUFR1J9t/5eCHtxz4q5A5BAGgQaKABAwgA5BBBY14hzQQCCJHFS1oaN1oAaBAA4IAM4CQJ0CCeKHUKIE6BCBIwgkodCoIHIIQIOEEoNAkDkFAAgYRUjRRzQAcgkDOEQgR9UIwfNCMalCMHKCYCgDA81MdpUAYGUADCAaoB2lANclBhWoU69N1OoCWmNCQfFaLfZtpa1S+3oikTqGEgeGnE+JVqM6lIzqUFW42baXNTpKtFrqhj2+IjQg8E/o20FTpOgb0hiXnUwIye5WiMjKEdpQaqNrQtmFtCkymIA9kRgCB+i0/0XZtqioyj0b+dNxZ36H/ADAVsjByVrbcUn1Oja8l+cAaRqgzLGuZuuEtIgg6Faqtnb181aYdme/v5rcBgZWDatM1DTFVpeMlsiQO5BX/AKLsoLegbG810SdQZH6qbXZllZum2t2UiBEtwYVik9lVu/TqB7ToWkELVQrGpvveWtYXlrOZgx+pCDc5jXgtcN5rgQQcghVxs+1bVp1GUg17HbwI1mHDPPDj4rKvUqUq9AgjonEsfjIJ0Pjj4hbyMjKDGpTZUYW1GhzTqHZC1GxtRT3RQptAGIER3eAW8jtKEYOUANDWgDQY1WmpZW1cDpaFN+hlzQVvjtKjRsk4QVqWzrOmS5ltSaSQTDeIJI/UnxW6rRp16T6VZgqU3iHNcJBCwbWmu2mAS1zS4PBBGDotwGuSg00bS3oOmjQp0yBujcaBjktN/ai4rWh6Nr2sqy/eEwIP84VyM6rQ+uTWoMoua7pJJM43REx4hBo2latdasY3fDG1G+xTJE+0OXBbRs61a2o0UoZVcHPZJ3XEcSPh8VsfULbllIj2Xgw6eI4eGfgVsMSBvZ5INT7Sg92+6mC6AJ7tFugIRjVYU61OqSKdQOgAmNM6IMwMDzQDCxc9rA3feGzgSQJWQGNSqNdaq2i3ecHGTAa3JJ7FNGoytTbVpmWPaHNPYVz9ptqOubRri8UDVz0Y9qdx2OwdyOZc29Zz6FtvOLw2ZEdHGIzwjTmVB0ab6dVofTeHtM5aZCyIyPNV7dpZdXDATukh47JEEfpPxVgjIygEBCMHzQjtKEYOUEwFAGB5qY7SoAwMoAGEA1QDtKAa5KBGfqkZSM6lIzqUAjI80ICEZGUI7SgEYPmpgKCMHKmO0oIAwPNAMIBgZQDtKABqkZ+qAa5KRnUoEZQjI80jOpQjIyqBAQjB80I7ShGDlQTAUAYHmpjtKgDAygAYSBlAO0pGuVQMxqhmEMwhmCoJgqBMBTnsUCYCAJQTlBMcEE50QMymZTMpmeCAZkIQUMyEM9iAZgqjXtK9F9e5sXMNxViW1B7JjTTP/wB9yvGYU57EFBlxehtQ1bcNDGsIgzOfa08VuqWnSuql1aoGVaZplgiOOdJnKsCYCCYQVNn0/s/S2weXim6Q52SQ7Oe2Z/Ra6WznsuaVV1ZruikBu5EjMTnUSc9pwrjKTGPe5jGtc8y4jiYWeZQU2U7y4ovp3baLHTLHUySJBkGD3BXM4ymZ4IZkKgQUMwhnsQzBUEwVourf7VaVKBeWB7Y3gMhb89igTAQULLZVO06M061YuZvZcR7UmXYiBJziFfE5QTHBBOUGq6oG4tqtHpCzpGFu83USqGzaAIDqdbfNJ9Qb5bh7Xw/ThqPBdTM8Fqt7ena0xTot3WjhJP7qivaWDqDGtrVulc2pvh8QdI5nhju56rdXtGVq9Oq8zuaN3GkHxE+BW8zIQz2KAZjVaLW0ZaNe2k526528Qc5Op+K3mYKnPYgwcxtRoD2tcORErIAxwQTAQTHBBWu6NaoWOZWZSZTO+ZplxmO8dqrUrt9WlTqms9gfvTvU2jd3ZknPYrd7SNezr0dwP6RhZul0TIjVU7dlWvsekKdGlT3mscGB5ILZBIJjjlBctOlNFr6j3OLxvAPaAWyBgwt5mQtdFlRm8ajt5znE64A4AfCFsMyEAgoZhDPYhmCgmCoEwFOexQJgIAlBOUExwQTnRAzKZlMymZ4IBmQhBQzIQz2IBmFMFQZgqc9iCBMBBKCYCCY4IAnKZlBOdEzKBmUMyEzPBDMhUCChmEM9iGYKgmCoEwFOexQJgIAlIOUExwTOVQJMaITjQpOPohOD5IGeSA4GCpnv8FAOB5KAD2IDrhAcfRAdfJA46JOdEnP0Sc/RAJyMIZ5ITkeSE9/ggE40KZ5ITg+Sme/wQQDgYKA9iA4HkgOPoqAOuE46IDr5JOfogTnRCcjCTn6ITkeSAZ5ITjQoT3+CE4PkgZ5IDgYKme/wUA4HkoAPYgOuEBx9EB18kDjok50Sc/RJz9EAnIwhnkhOR5IT3+CATjQpnkhOD5KZ7/BBAOBgoD2IDgeSA4+ioc8KGNDAGtbutAAAGgCkHXySc/RAnOiE5GEnP0QnI8kAzyQnGhQnv8EJwfJAzyQHAwVM9/goBwPJQAexAdcIDj6IDr5IHHRJzok5+iTn6IBORhDPJCcjyQnv8EAnGhTPJCcHyUz3+CCAcDBQHsQHA8kBx9FQB1wnHRAdfJJz9ECc6ITkYSc/RCcjyQDPJCcaFCe/wQnB8kDPJAcDBUz3+CgHA8lAB7Ek5wgOPok6qgSI1QkRqh0UnQqCJHNARAypQaBFQCOaAjOVI0QcUREidUkTqp4pxQQSJGUJHNSdQhQQSI1SRzUnQoggEQMoCOakaBBoggEZykidVI4pxQRInVCRIyp4odQggkc0JEaqSh0KCJHNARAypQaBFQCOaAjOVI0QcUREidUkTqp4pxQQSJGUJHNSdQhQQSI1SRzUnQoggEQMoCOakaBBoggEZykidVI4pxQRInVCRIyp4odQggkc0JEaqSh0KCJHNARAypQaBFQCOaAjOVI0QcUREidUkTqp4pxQQSJGUJHNSdQhQQSI1SRzUnQoggEQMoCOakaBBoggEZykidVI4pxQRInVCRIyp4odQggkc0JEaqSh0KCJHNARAypQaBFQCOaSM5UjRRzRAgRohAg4SBH1QjB80EwOQUACBhTAUAYHmgADkEAGcBAMIBqgQJ0CQJ0CRn6pGUAgSMIQOQQjI80ICAQIOFMDkFBGD5qYCCABAwgA5BAMDzQDCAAM4CQJ0CAapGfqgQJ0CECRhIyhGR5oBA5BCBBwhAQjB80EwOQUACBhTAUAYHmgADkEAGcBAMIBqgQJ0CQJ0CRn6pGUAgSMIQOQQjI80ICAQIOFMDkFBGD5qYCCABAwgA5BAMDzQDCAAM4CQJ0CAapGfqgQJ0CECRhIyhGR5oBA5BCBBwhAQjB80EwOQUACBhTAUAYHmgADkEAGcBAMIBqgQJ0CQJ0CRn6pGUAgSMIQOQQjI80ICAQIOFMDkFBGD5qYCCABAwgA5BAMDzQDCAAM4CQJ0CAapGfqgQJ0CECRhIyhGR5oBA5BCBBwhAQjB80EwOQUACBhTAUAYHmgADkEgZwgGEgZQCMalCMHK+YdY22Pw9j+R/qTrG2wf9XsfyP9SD6hHaVAGBlfMOsfbH4ex/I/1J1jbYj+z2P5H+pFfTwMalANclfMOsbbH4ex/I/1J1jbY/D2P5H+pEfT4zqkZ1K+YdY22Pw9j+R/qTrG2xP9nsfyP9SD6eRkZQjtK+YdY22Pw9j+R/qQ/wD+jbYP+r2P5H+pB9PIwcqY7Svl/WNtg/6vY/kf6k6xtsfh7H8j/Ug+ngYGUA7SvmHWNtiP7PY/kf6k6xtsfh7H8j/Ug+nga5KRnUr5h1jbY/D2P5H+pOsbbH4ex/I/1IPp8Z1KEZGV8w6xtsT/AGex/I/1J1jbY/D2P5H+pB9PI7ShGDlfMOsbbH4ex/I/1J1jbYP+r2P5H+pB9QjtKgDAyvmHWNtj8PY/kf6k6xtsR/Z7H8j/AFIr6eB2lANclfMOsbbH4ex/I/1J1jbY/D2P5H+pEfT4zqUjOpXzDrG2x+HsfyP9SdY22J/s9j+R/qQfTyMjKEdpXzDrG2x+HsfyP9SdY22Pw9j+R/qQfTyMHKmO0r5f1jbYP+r2P5H+pOsfbH4ex/I/1Ir6eBgZQDGpXzDrG2xH9nsfyP8AUnWNtj8PY/kf6kR9PA1yUjOq+YdY22Pw9j+R/qTrG2x+Hsvlu9SD6fGdShGRlfMOsbbHuLL5bvUnWNtj3Fl8t3qQfTyO0oRg5XzDrG2x7iy+W71J1jbY9xZfLd6kH1CO0qAMDK+YdY22Pw9l8t3qTrG2x7iy+W71Ir6eBjUoBrkr5h1jbY9xZfLd6k6xtse4svlu9SI+nxnVIzqV8w6xtse4svlu9SdY22PcWXy3epB9PIyMoR2lfMOsbbHuLL5bvUnWNtj3Fl8t3qQfTyMHKmO0r5f1jbY9xZfLd6k6xtsfh7L5bvUivp4GBlAMalfMOsbbHuLL5bvUnWNtj3Fl8t3qRH08DXJSM6r5h1jbY9xZfLd6k6xtse4svlu9SD6fGdShGRlfMOsbbHuLL5bvUnWNtj3Fl8t3qQfTyO0oRg5XzDrG2x7iy+W71J1jbY9xZfLd6kH1CO0qAMDK+YdY22Pw9l8t3qTrG2x7iy+W71Ir6eBjUpGuV8w6xtse4svlu9SjrF2x7iy+W71IjyKIiKIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiD/2Q==\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/h0XUxE47l00\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f852a19ffd0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('h0XUxE47l00')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans notre cas on veut toujours résoudre notre système matriciel non linéaire.\n",
    "\n",
    "On va donc définir la méthode itérative par analogie avec la méthode 1D. En 1D on a :\n",
    "\n",
    "$$\n",
    "x^{k+1} = x^k - \\frac{f(x^k)}{f'(x^k)}  \\quad \\textrm{ou bien} \\\\\n",
    "f'(x^k) \\; (x^{k+1} - x^k) = - f(x^k)\n",
    "$$\n",
    "                               \n",
    "Cela donne en $n$ dimensions :\n",
    "\n",
    "$$\n",
    "J_{\\bf f}({\\bf x}^k) \\; ({\\bf x}^{k+1} - {\\bf x}^k) =  - {\\bf f}({\\bf x}^k)\n",
    "$$\n",
    "\n",
    "avec $J_{\\bf f}$ la matrice Jacobienne de ${\\bf f}$ :\n",
    "\n",
    "$$\n",
    "J_{\\bf f}\\left({\\bf x}\\right)=\n",
    "\\begin{pmatrix} \n",
    "\\dfrac{\\partial f_1}{\\partial x_1} & \\cdots & \\dfrac{\\partial f_1}{\\partial x_n} \\\\\n",
    "\\vdots & \\ddots & \\vdots \\\\\n",
    "\\dfrac{\\partial f_n}{\\partial x_1} & \\cdots & \\dfrac{\\partial f_n}{\\partial x_n}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "et pour notre système non linéaire\n",
    "\n",
    "$$\n",
    "{\\bf f}({\\bf x}) = A({\\bf x})\\, {\\bf x} - {\\bf b}\n",
    "$$\n",
    "\n",
    "On notera que ${\\bf f}$ est une fonction vectorielle à savoir qu'elle retourne\n",
    "un vecteur."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test\n",
    "\n",
    "On reprend le même système matriciel :\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "x_0 - 2 x_1 &  x_1 \\\\\n",
    "x_0 & 2 x_0 + x_1 \\\\\n",
    "\\end{bmatrix}\n",
    "\\;\n",
    "\\begin{bmatrix}\n",
    "x_0 \\\\\n",
    "x_1 \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 \\\\\n",
    "9 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "La matrice Jacobienne de la fonction ${\\bf f}$ définie ci-dessus est : \n",
    "\n",
    "$$\n",
    "J_{\\bf f}({\\bf x}) = \n",
    "\\begin{bmatrix}\n",
    "2 x_0 - 2 x_1 &  2 x_1  - 2 x_0\\\\\n",
    "2 x_0 + 2 x_1 & 2 x_0 + 2 x_1 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^0 =  [0.80714286 1.40214286]  erreur² =  0.1738391124362817\n",
      "x^1 =  [0.82661178 1.47589539]  erreur² =  0.14013193945500102\n",
      "x^2 =  [0.84449678 1.53832416]  erreur² =  0.11305725653981773\n",
      "x^3 =  [0.86066582 1.59186587]  erreur² =  0.09127299702681221\n",
      "x^4 =  [0.87518409 1.63820488]  erreur² =  0.07372379277598218\n",
      "x^5 =  [0.88818095 1.67857971]  erreur² =  0.0595731611254003\n",
      "x^6 =  [0.89980142 1.71393945]  erreur² =  0.04815469077192336\n",
      "x^7 =  [0.91018754 1.74503332]  erreur² =  0.038935573039892035\n",
      "x^8 =  [0.91947123 1.77246602]  erreur² =  0.03148872349291875\n",
      "x^9 =  [0.92777209 1.79673415]  erreur² =  0.025471148220392294\n",
      "x^10 =  [0.9351973  1.81825117]  erreur² =  0.020606983474307096\n",
      "x^11 =  [0.94184231 1.83736516]  erreur² =  0.01667410290306686\n",
      "x^12 =  [0.94779187 1.85437189]  erreur² =  0.013493486464829903\n",
      "x^13 =  [0.9531212  1.86952455]  erreur² =  0.010920748632859209\n",
      "x^14 =  [0.95789703 1.88304134]  erreur² =  0.00883936731257032\n",
      "x^15 =  [0.96217855 1.89511126]  erreur² =  0.007155259635295223\n",
      "x^16 =  [0.96601835 1.90589887]  erreur² =  0.005792428863449946\n",
      "x^17 =  [0.96946319 1.91554793]  erreur² =  0.0046894658536456735\n",
      "x^18 =  [0.97255463 1.92418452]  erreur² =  0.003796734034536249\n",
      "x^19 =  [0.97532973 1.93191956]  erreur² =  0.0030741021940233564\n",
      "x^20 =  [0.97782148 1.93885083]  erreur² =  0.0024891170269789866\n",
      "x^21 =  [0.98005934 1.94506477]  erreur² =  0.002015529173331623\n",
      "x^22 =  [0.9820696  1.95063793]  erreur² =  0.0016321037099287408\n",
      "x^23 =  [0.98387574 1.95563824]  erreur² =  0.0013216597493956443\n",
      "x^24 =  [0.98549877 1.96012605]  erreur² =  0.0010702947077985679\n",
      "x^25 =  [0.98695747 1.96415505]  erreur² =  0.0008667575172292292\n",
      "x^26 =  [0.98826867 1.96777309]  erreur² =  0.0007019420351895113\n",
      "x^27 =  [0.98944742 1.97102285]  erreur² =  0.0005684774962019743\n",
      "x^28 =  [0.99050723 1.9739424 ]  erreur² =  0.0004603973426083823\n",
      "x^29 =  [0.99146019 1.9765658 ]  erreur² =  0.0003728713824135844\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return A(x) @ x - b\n",
    "\n",
    "def J_f(x):\n",
    "    return 2 * np.array([[x[0] - x[1], x[1] - x[0]],\n",
    "                         [x[0] + x[1], x[0] + x[1]]])\n",
    "\n",
    "x = np.array([1, 1.1])\n",
    "for i in range(30):\n",
    "    delta = lin.solve(J_f(x), -f(x))\n",
    "    x = x + delta\n",
    "    print(f'x^{i} = ', x, end='  ')\n",
    "    print(f'erreur² = ', np.square(A(x) @ x - b).sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On converge alors que la méthode du point fixe oscille \n",
    "sans convreger pour cette valeur initiale.\n",
    "\n",
    "On converge néanmoins moins vite qu'avec la méthode du point fixe sur-relaxée.\n",
    "\n",
    "**Remarques** sur la méthode de Newton-Raphson\n",
    "\n",
    "* Le coût de construction de la matrice Jacobienne peut être très lourd\n",
    "* Une technique pour aller plus vite consiste à dire que d'une itération à une\n",
    "  autre on n'a pas beaucoup bougé et donc la matrice Jacobienne n'a pas du trop \n",
    "  changer. Aussi on recalcule cette matrice que toutes les 3 ou plus itérations\n",
    "* Il s'agit d'une matrice pleine qu'il rend plus compliqué la résolution du système\n",
    "  lors de l'iteration (une méthode de gradiant ne marchera pas en général)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}