{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# ANU - TP5 Optimisation\n", "\n", "B. Boutin, Y. Le Hénaff et M. Bouchereau" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## I - Gradient à pas fixe\n", "\n", "Cette première partie propose d'implémenter la méthode du gradient à pas fixe afin de résoudre des systèmes linéaires de la forme\n", "\n", "\\begin{equation}\n", " Ax=b,\n", "\\end{equation}\n", "\n", "où $A\\in\\mathcal{M}_d(\\mathbb{R})$ est une matrice symétrique, via l'optimisation de fonctionnelle quadratique suivante:\n", "\n", "\\begin{equation}\n", " J(x) = \\frac12\\langle Ax | x \\rangle - \\langle b | x \\rangle.\n", "\\end{equation}\n", "\n", "On utilise pour cela la méthode du gradient à pas fixe:\n", "\n", "\\begin{equation}\n", " x_{k+1} = x_k - \\rho\\nabla J(x_k)\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1 - Etude générale\n", "\n", "On se propose d'implémenter la méthode du gradient à pas fixe afin de résoudre un tel système.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q1 -** Ecrire une fonction permettant de résoudre un système linéaire de cette forme via le gradient à pas fixe. On veillera à ce que la méthode s'arrête lorsque le nombre d'itérations ne dépasse pas un seuil maximal et/ou lorsque que la suite $(\\nabla J(x_k))_{k\\in\\mathbb{N}} = (Ax_k-b)_{k\\in\\mathbb{N}}$ ne varie presque plus. Cette fonction devra retourner les itérations $x_0,x_1,\\cdots,x_N$ et la solution approchée $x_N$ où $N$ est le nombre d'itérations.\n", "\n", "Tester cette fonction dans le cas suivant:\n", "\n", "$$A = \\begin{bmatrix}\n", " 1 & 0 \\\\\n", " 0 & \\frac{1}{20}\n", " \\end{bmatrix} , \\hspace{10mm} b = \\begin{bmatrix}\n", " 1 \\\\\n", " \\frac{1}{10} \n", " \\end{bmatrix}$$\n", " \n", "en prenant $\\rho = 1$, $\\verb|Max_iter| = 1000$ et $\\verb|tol| = 10^{-8}$. Représenter l'erreur $\\|Ax_k-b\\|$ en échelle logarithmique ($\\log\\|Ax_k-b\\|$ en fonction de $k$, les itérations). Procéder de même en prenant $\\rho = 0.5$ et $\\rho = 1.5$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution - rho = 0.5: [1. 1.9999998]\n", "Solution - rho = 1.0: [1. 1.99999981]\n", "Solution - rho = 1.5: [1. 1.99999981]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnmElEQVR4nO3dfbRU9X3v8feXpwCCGB7k2hzJwWBMEAQRTaw1C2tuFJIYb1ZuIj02ZmlKfEDT9ubmarmJ1ZY0fVhdzWoEyvWpuSGwsoyNhljTSD03bUqCYCB6MFZaBI8P5akiiAic871/zB6YM8zDntl7z977zOe11qwzs2dmz+8nx/me/fv+vr+fuTsiIiLVDEm7ASIikm0KFCIiUpMChYiI1KRAISIiNSlQiIhITcPSbkASJk6c6J2dnU2998033+SUU06Jt0Etpj6kL+/tB/UhK1rZh02bNu1x90nlxwdloOjs7GTjxo1Nvbe7u5t58+bF26AWUx/Sl/f2g/qQFa3sg5ntqHRcQ08iIlKTAoWIiNSkQCEiIjUNyhyFiEjcjh49Sm9vL4cPH27p544bN47nnnsu1nOOHDmSjo4Ohg8fHur1ChQiIiH09vYyduxYOjs7MbOWfe6BAwcYO3ZsbOdzd/bu3Utvby9Tp04N9R4NPRWtAjqBTcHPVWk2RkSy5vDhw0yYMKGlQSIJZsaECRMaujLSFQUUgsIi4FDweEfwGKArlRaJSAblPUgUNdoPXVEALOFEkCg6FBwXEWlzChQAOxs8LiKSgscff5xzzjmHadOm8fWvf/2k57u7uxk3bhyzZ89m9uzZ3H333bF8roaeAKZQGG6qdFxEJAP6+vq45ZZb+PGPf0xHRwcXXnghV111FdOnTx/wuksvvZS1a9fG+tm6ogBYCowuOzY6OC4i0oRVz6yi8686GXLXEDr/qpNVz0SbIbNhwwamTZvGWWedxYgRI7jmmmt45JFHYmptbQoUUEhYrwTeHTx+d/BYiWwRacKqZ1ax6AeL2LF/B46zY/8OFv1gUaRg8fLLL3PmmWcef9zR0cHLL7980uvWr1/PrFmzmD9/Pj09PU1/XikFiqIu4EXgguCngoSINGnJuiUcOjpwhsyho4dYsq75GTLuftKx8tlLc+bMYceOHWzZsoVbb72Vq6++uunPK6VAUWoV8AyF/yqdqJZCRJqyc3/lmTDVjofR0dHBSy+9dPxxb28vv/ZrvzbgNaeeeipjxowBYMGCBRw9epQ9e/Y0/ZlFChRFxVqKI4BzopZCwUJEGjRlXOWZMNWOh3HhhRfywgsvsH37do4cOcKaNWu46qqrBrzmtddeO37lsWHDBvr7+5kwYULTn1mkQFGkWgoRicnSy5cyevjAGTKjh49m6eXNz5AZNmwY3/zmN7niiit4//vfz6c//WnOPfdcVqxYwYoVKwB46KGHmDFjBrNmzeK2225jzZo1sRQJanpskWopRCQmXTMLSc4l65awc/9OpoybwtLLlx4/3qwFCxawYMGCAcduvPHG4/cXL17M4sWLI31GJQoURaqlEJEYdc3sihwYskJDT0WqpRARqUhXFEXFwL8PMApXEkvRNFkRaXu6ohARkZoUKIo0PVZEpCIFiqIvH9T0WBGRChQoALavglfLM9kBTY8VkYy4/vrrOf3005kxY0bF592d2267jWnTpnHeeefx9NNPx/K5ChQAW5bAhCoRQdNjRSQjPve5z/H4449Xff7v//7veeGFF3jhhRdYuXIlN910Uyyfq0ABcGgnfPoPYMSbA4+PeFPTY0WkOdtXwfc74TtDCj+3R094fuhDH2L8+PFVn3/kkUf47Gc/i5nxwQ9+kNdff51XX3018ucqUACMngKXrIbP/w4MOwL0w8QX4eY7ND1WRBq3fRVsWASHdgBe+LlhUSzBopawS5E3SoEC+OfxC3izf+CxfodfjZmdSntEJOe2LIG+stkxfYcKxxMUZinyZmQ+UJjZWWZ2n5k9lNRnXPv0Yzzw2EL67/0/cGwEMIQhezuZ8qfXaHqsiDTuUJWcZ7XjMQmzFHkzEg0UZna/me0ys2fLjl9pZs+b2TYzu73WOdz93939hiTbuXP/Tj722NcYcuSUAcdHHx2t6bEi0rjRVWbBVDsek6uuuopvfetbuDs/+9nPGDduHGeccUbk8ya9hMeDwDeBbxUPmNlQ4B7gvwK9wFNm9igwFPiTsvdf7+67Em4jU8ZNYcr+Kv+Amh4rIo2atbSQkygdfho6unA8goULF9Ld3c2ePXvo6Ojgrrvu4ujRo0BhFdkFCxbw2GOPMW3aNEaPHs0DDzwQ6fOKrNKYVpzMrBNY6+4zgscXA3/o7lcEj+8AcPfyIFF+nofc/VM1nl9EoZaayZMnX7BmzZrQbdz31j7GPD+GEX0jONhxkDG9Y44/1z+8nyHnZX6EboCDBw8e3+Uqr/Leh7y3H9SHcuPGjWPatGmhXz+s97u841d3YW/14qM6ePt9d3Ks49MNf25fXx9Dhw5t+H31bNu2jf379w84dtlll21y97knvdjdE71R2FT02ZLHnwLuLXn828A3a7x/ArAC+DfgjjCfecEFF3ijVlz7294/4qA/+RdPHj9TP32+7pIVDZ8rbU8++WTaTYgs733Ie/vd1YdyW7duje1cjXjjjTcSOW+l/gAbvcJ3ahp/KldKwVe9rHH3ve5+o7u/x+tcdUSx8CP/F/vQwMs0YwiXPXUt93/p/qQ+VkQk89IIFL3AmSWPO4BXUmjHAGOHAJs/dtJxO3IKCx/4Tewu4+Yf3tz6homIpCyNQPEUcLaZTTWzEcA1wKMptONkeyontEfum8JfT4TlG5cz/I+Gs+oZzZkVkfaR9PTY1cB64Bwz6zWzG9z9GLAY+BHwHPBdd+9Jsh1h9ANMrDzFySbu5JbT4K8nwrH+Y1z78LWM+uNRChgi0hYSDRTuvtDdz3D34e7e4e73Bccfc/f3BnmHTKym9ML4y/FZays84zB7LWYcDxYAh/sOc+3D1/Lhb324lc0UEWm5fM37TND75j/Bsc1XVXjG4GefKdwrCxYA67avU/5CRFqi3jLj3d3djBs3jtmzZzN79mzuvvvuWD5XgaLE8H1Viu4OToSfLgROBIuFZVOzl29croAhIomqt8w4wKWXXsrmzZvZvHkzX/3qV2P5XAWKUlWr6w3uW3HikcG9kyu/UglvEQEK68R1UviW7SSWdePqLTOeFAWKUrWyJW+PPX5VATDKYPOZlV+qhLdIm1tFYZ2IYJVxdgSPW/B1sH79embNmsX8+fPp6YlnnpACRamae0+cfFVx3juqBws4kfBWwBBpM0uAslXGOUTii4zOmTOHHTt2sGXLFm699VauvvrqWM6rQFGu1jKJZVcVYYIFaIaUSNuptphowouMnnrqqcfXtlqwYAFHjx5lz549kc+rQFGu5pe+wd/cO/BIECweD7GS77rt65S/EGkH1fKdya4yzmuvvXZ886INGzbQ39/PhAkTIp9XgaLceKDWYpN9o+BrA2cdmMFHThk4bbYa5S9E2sBSYHTZsdHUzoOGsHDhQi6++GKef/55Ojo6uO+++1ixYgUrVhSGxR966CFmzJjBrFmzuO2221izZk0sO9wlvR9FPq0Arq32pEHPRwpDUJesPnE0mDYLcGuIK73icNQDv3iAJz77RMQGi0imFPOdSygMN02hECRq5kHrW716dc3nFy9ezOLFi6N9SAW6oqiki9pXFWWJ7eNHKxTk1aOCPZFBqgt4kcL6QC8SOUikSYGimpPjwEBlie2iZoIFqGBPRLJLgaKaLuDyWi84ObF9/BmDW955cvV2GCrYE8muYqI47xrthwJFLU8AI2s8XyGxXWTAt84YyoghIxr+WCW8RbJn5MiR7N27N/fBwt3Zu3cvI0fW+nIbSMnseu6l4cR20TD6eHv2Oazq/AOu//71HOk/0tBHFxPen3/k89z7iXvpmpnjQU6RnOvo6KC3t5fdu3e39HMPHz7c0Jd6GCNHjqSjoyP06xUo6ukCbgQOVntBkNiuECgAOLCVrhe/RtdX3mbVM6u47uHr6KOvoSZohpRI+oYPH87UqVNb/rnd3d2cf/75Lf/cUhp6CqPJxPZxB7bCD86la2YXx+48xk1zb2qqGZohJSJpUKAII0xiu8J02QEObIUnCkt4LPvoMvxObzpgaIaUiLSSAkVY9RLbb4+FB/669jl2rYPtJ5LTxYBx+dSaUaiq5RuXc+495zb1XhGRsBQoGlF5NmzA4Ilbag9BAay/7qRDT3z2Cb79yW83NUNq656turoQkUQpUDQiVMX2yjon6YMfnHwV0DWzi7e/8nbTAUPDUSKSFAWKRtVNbI+BTV+q/ZoguV1JacAYytCGm7d843KefvVp1V+ISGwUKBpV96oC+Js/r3+ekuR2xY+JMEPKcRXsiUhsFCiaUe+q4k3ge0/WP8+udbCh9lBRlIS3NkwSkTgoUDSjC6j3h/7D8+oPQQFsW143WMCJhHczw1GqvxCRKBQomrWMcENQY6fXP1fIYBG1YE8JbxFphgJFFGGGoN7oiTVYQDwFe1qhVkTCUqCIIkxi+wvAx3tgSIhFvRoIFhAtf6EVakUkLAWKqMJcVdwMfKBmtd4J25YPqN4OI0rBXjHhrYAhItUoUEQVJrG9HPiXLpgWcqho/fWNN6Ok/qIZxYCh/IWIlFOgiEOYxPYXgIuWhQwWR6oW5NXTNbOLC864IFL+QlcXIlJKgSIuYYagVhE+WNSo3g4jSsJbw1EiUkqBIi5hE9vQsmABKtgTkegUKOIU9qoCGgsWNZb6CCtKwlsFeyLtTYEiTmES218ouX/RMjg9xF/6IZb6CCOOBQcVMETajwJF3OoltkuvKgA+/ETsBXn1FCu8o2yYpII9kfahQJGEekNQXyh7/PH4q7fDiDIcpYI9kfahQJGErjrPl19VQGrBIuqGSZohJTL4KVAkZUKd58uvKqCxpT4arN6uJ2r+QjOkRAYvBYqkfKPO88WlPcqFXeqjiertMKKuUKsZUiKDT+YDhZm938xWmNlDZtbct1cawi7tUX5hMDXsUh/NV2+HEccKtcPuHqbhKJFBINFAYWb3m9kuM3u27PiVZva8mW0zs9trncPdn3P3G4FPA3OTbG/swi7tUa7FNRa1RCnY6/M+DUeJDAJJX1E8CFxZesDMhgL3APOB6cBCM5tuZjPNbG3Z7fTgPVcB/wysS7i98WukCK9U2GARU41FPSrYE2lfiQYKd/8JsK/s8EXANnf/d3c/AqwBPuHuz7j7x8puu4LzPOruv079+UTZ08jSHuXCBosEktuVRJ0hpYI9kXwyd0/2A8w6gbXuPiN4/CngSnf/fPD4t4EPuPviKu+fB3wSeAfwS3e/p8rrFgGLACZPnnzBmjVrmmrvwYMHGTOm3jd7g/YB2+u8ZiowvspzB/4Vjh6ocwKD8XOAhPpQwb639rH99Xodq8wwOk/rZPyoyp1uVR+Skvf2g/qQFa3sw2WXXbbJ3U8a4h/Wkk8fyCocqxqt3L0b6K53UndfCawEmDt3rs+bN6+pxnV3d9Pse2u6mULyupoRwNvVnpwHq4eDH6v9GQemw8d7kutDFTf/8GaWb6zVuepGDh3JvZ+4l66ZAy8WW92HuOW9/aA+ZEUW+pDGrKde4MySxx3AKym0o7XqJbaPALVyvh98sP5nxLDabDOKCe/pE0MUDJZR/YVI9qURKJ4CzjazqWY2ArgGeDSFdrRevcT2OiontiH8tNkDW2F/T4MNi0fPLT1NF+wp4S2SXUlPj10NrAfOMbNeM7vB3Y8Bi4EfAc8B33X3dL7ZWi1KYhvCJ7f7Dic+bbaaqAV7xYT3zv07Y26ZiDQr6VlPC939DHcf7u4d7n5fcPwxd3+vu7/H3Zcm2YbMaXa6bFGLlyZvVtSCvd2HdmuFWpGMyHxl9qDTBdT7nq+3OkcKS5M3K0rBnlaoFckGBYo0PAHUWvuvXmIbUltttllRCvaKCe9z72l9ol5EFCjSU2/tv1qJ7aKcBYuoK9Ru3bNV60eJpECBIi1RE9tFKS5N3qwoCe/i+lEajhJpHQWKNEVNbBelvDR5s6IkvLVhkkjrKFCkKY7ENmRmafJmRUl4q2BPJHkKFGmLI7ENjS1NnsFgAVqhViSr6gYKMxtqZt9uRWPaVhyJbRgUwaKY8J562lSs4rJgtWmFWpH41Q0U7t4HTAqW25AkxJXYhkKwGDmp/utasOlRFONHjaf/zv6mhqOgEDBUsCcSj7BDTy8CPzWzr5jZ7xdvCbar/cSV2AYYPSUX1dthRBmOUsGeSDzCBopXgLXB68eW3CQuYRLbYa8qIFfV2/VE3TBJM6REogkVKNz9rkq3pBvXduoltt+ksK9FWDkryKsnasGeZkiJNCdUoDCzJ83sH8tvSTeuLdVLbC8n/BAUDLpgAdFXqNUMKZHGhB16+hLwP4PbV4DNwMaE2tTe4kxsF+WwejuMqCvUKuEtEk7YoadNJbefuvvvAx9IuG3tK87EdlHY6u2fNxqF0hfHCrUajhKpLuzQ0/iS20QzuwL4Lwm3rX0lcVURtnq7/83cDEGVU8GeSDLCDj1tojDUtInCjnX/A7ghqUYJ4a4qGv1OC1uQl6N8RbmoM6RUsCdysrBDT1Pd/azg59nu/hF3/+ekG9fWuoB63+mNJrahLYIFRJ8hpfyFyAlhh55Gm9n/NrOVweOzzexjyTZNWEb9IagvNnHei5bBsHonJnfJ7UqKM6S0w55I88IOPT1AYXm6Xw8e9wJ/nEiLZKB6Q1B7mzzvhfVOHMjY0uTNimOHPSW8pV2FDRTvcfc/A44CuPtb0MSKbdK4MIntZv7YzfnS5M2IOhylhLe0q7CB4oiZjQIcwMzeA7ydWKtkoHp//Df7R/8gWG22GVEL9pTwlnYTNlDcCTwOnGlmqygsfP3lxFolA9W7qgi7Z0UlbRosQAV7ImGF2Y9iCPBO4JPA54DVwFx37060ZTJQvauKsHtWVNJIsMjw0uTNiqNgTwlvGczC7EfRDyx2973u/kN3X+vue1rQNimVRBFeqYuWDZqlyZulhLdIZWGHnn5sZl8yszNLq7QTbZmcLImlPUoNoqXJmxVHwnv4Hw1n31v7EmidSDrCBorrgVuAn1Cozi5Waksrxb1nRSWDcLXZZkRJeB/rP8b217drOEoGjbA5ituDquzS21ktaJ+UC7Nnxc6In6FgcVyUhLc2TJLBImyO4pYWtEXCqrcQ7G6iDUHBoF2avFlREt7KX0jeKUeRR0kntovCLk0+SKq3wygmvFWwJ+1EOYq8SjqxDW1ZvR2GCvak3TSyeqxyFFnSqquKNi7IqyfKcBSoYE/yo2agMLMvl9z/72XPfS2pRklISexZUYmCRU1R6i9UsCd5UO+K4pqS+3eUPXdlzG2RRiW1Z0UlbV69XU/UDZM0Q0qyrF6gsCr3Kz2WNITZsyKubbBVvV1XMWBMPW1qUwlvzZCSLKoXKLzK/UqPJS2tSGwXqXo7lPGjxkdKeK/bvo5hdw/T1YVkQr1AMcvM3jCzA8B5wf3i45ktaJ+E0arEdlEjBXltUGNRS5SCvT7v03CUZELNQOHuQ939VHcf6+7DgvvFx8Nb1UgJoZVXFRA+WKy/LsYPzS8V7Emeha2jkKwLk9iOuy4uVPV2X1smt6uJMkNKBXuSFgWKwaReYjvKBkfVhKnebuPkdiVRZ0ipYE9aLfOBwszmmdk/mdkKM5uXdnsyL8kNjioJW73d5sntSkoDhjUxiVAFe9IqiQYKM7vfzHaZ2bNlx680s+fNbJuZ3V7nNA4cpLBmam9SbR00uqj/rxpnYhvCT5tVsKioa2YX/Xf2M31iiJxPGRXsSSskfUXxIGWFeWY2FLgHmA9MBxaa2XQzm2lma8tupwP/5O7zgf8F3JVweweHd9d5Pu7ENmjabAx6bunRDnuSSeaebDmEmXUCa919RvD4YuAP3f2K4PEdAO7+J3XOMwL4jrt/qsrzi4BFAJMnT75gzZo1TbX34MGDjBlTb65pth08eJAxr4yBAzVeZMCcBD78P5+GML9TY6bCiOoLEOf93yFq+/e9tY/tr29v+v2TRk9iyrgpTb8f8v9vAOpDoy677LJN7j63/Piwlnz6QO8CXip53At8oNqLzeyTwBXAacA3q73O3VcCKwHmzp3r8+bNa6px3d3dNPverOju7mbeonkwCjhc44WXU9gIKU7bX4b119Z/3eER8FtvV3067/8OcbX/5h/ezPKNy5t+/01zb2LZR5c19d68/xuA+hCXNJLZlbJ2Vf8EdfeH3f0L7v4Zd+9OrlmDUL0JSXEntkFLk8csSsEeKOEt8UgjUPQCZ5Y87gBeSaEdg1+rK7aLtNps7KIU7BUT3ufeo//W0pw0AsVTwNlmNjXIO1wDPJpCO9pDqyu2ixQsEhGlYG/rnq2qv5CmJD09djWwHjjHzHrN7AZ3PwYsBn4EPAd81917kmxHW+uikIuoJamdTLU0eSJK6y+aWaFWw1HSqEQDhbsvdPcz3H24u3e4+33B8cfc/b3u/h53X5pkG4RCwrrWShtJVGwXaWnyxETZklX1F9KIzFdmS0zSSGwXqcYiUVES3towScJQoGgXaSW2ixpZmlzBoilaoVaSokDRTtJKbBc1EiwO7UywIYNbMeHdzPpRWqFWKlGgaCdhEttfTLgNoZYmBw7v1pVFBMX1o6LUX2x6dZMChgAKFO2nXmJ7bwvaEGZpctAOeTGIMhwFmiElBQoU7aje93TSw9Shq7eB9UnN3W0vUeovNENKFCjaUVed55OcAVUUtsZCS33EJuqGSZoh1b4UKNrVhDrPJzkDqkjV26mIWrCnGVLtR4GiXX2jzvNJz4AqUrBITZSCPSjMkFL+oj0oULSrNJf2KNdIsNBMqNhFKdhT/qI9KFC0szSX9ijXyHaqkggV7Ek1ChTtLs2lPcqFXepDQ1CJUsGelFOgaHdpL+1RLkz1tvIVieua2cWcM+Y0PUNq+cblChiDiAKFpL+0R7mP91B5I8QSChYtUZwhpR322psChWQrsV00prP+a7SPRcss++gyFey1MQUKKchSYhtgxHjtY5ExcRXsKeGdPwoUckKWEtugfSwyKmrBnhLe+aNAISdkLbEN2sciw6IW7CnhnR8KFDJQ1hLboGCRcVEK9kAJ7zxQoJCBwiS2W31VAeH3sdDS5KmJUrBXTHgrf5FNChRysnqJ7TeBNP5wD7uPhZYmT1WUJc2Vv8gmBQqprN538nJaPwQVeh8LLU2etqgJbw1HZYsChVSWxcQ2aLXZnImS8Fb9RXYoUEh1WUxsg4JFDkVJeGvDpPQpUEh1Wb2qgMaChaq3M6MYMKZPDDGLrYwK9tKjQCG1hbmqSCvvGHZpclVvZ07PLT2RNkxSwru1FCikti6g3v/PaSS2i1S9nVtx1F8oYLSGAoXUt4zsDkGBCvJyLkr9BWiGVCsoUEg4WU1sFylY5F6U+gvNkEqWAoWEk+XEdpGqt3MvrhVqFTDipUAh4WX9qgJUvT1IRC3YKwaMnft3JtC69qNAIeGFSWyn/f2r6u1BJeoKtbsP7dbVRQwUKKQx9RLbrd7gqBLVWAw6KthLlwKFNK7eEFSrNziqJGywUI1FrkSZIaWCveYpUEjj8pDYhvDBQsnt3NEKta2lQCHNyUNiG8JXbyu5nTtxrFCrgBGOAoU0J8wGR1n57v3wE2DD6rxIye28Kia8VbCXHAUKaV69DY6OkN46UOU++GD912i12VxTwV5yFCgkmjAbHGVB2GmzCha5FlfBnhLeAylQSDRhEttZ+QNN02bbRjFgTD1talP5CyW8B8p8oDCzS81shZnda2b/knZ7pIJ6ie0szIAq0tLkbWX8qPGRCvaU8C5INFCY2f1mtsvMni07fqWZPW9m28zs9lrncPd/cvcbgbXA3ybZXmlSF7VzFWnuWVGJliZvO3EsaX7uPe07JJn0FcWDwJWlB8xsKHAPMB+YDiw0s+lmNtPM1pbdTi95628BqxNurzQrTK4iK0NQoNVm21SUgr2te7a27dWFuXuyH2DWCax19xnB44uBP3T3K4LHdwC4+5/UOMcU4Cvu/js1XrMIWAQwefLkC9asWdNUew8ePMiYMfUG3bMttT78Auiv8fwQ4Pxwp2pZH/b3QN/h+q8bOQlGTwl9Wv0eZUOtPux7ax8vvv4iTnPfgZNGT2LKuPC/E81q5b/DZZddtsnd55Yfrze5PAnvAl4qedwLfKDOe24AHqj1AndfCawEmDt3rs+bN6+pxnV3d9Pse7MitT68DFxb5zXfpjBUVUfr+jAP1oyC/jrB4jDw/m8XZk+FoN+jbAjTh1XPrOK6h6+jj76Gzz9syDAevPpBumaG+71oRhb+HdJIZluFYzVDurvf6e5KZGddXpb2KKelydtalBVq26X+Io1A0QucWfK4A3glhXZIEvKytEcpLU0uxLPg4GANGGkEiqeAs81sqpmNAK4BHk2hHZKEPOxZUUkjNRYKFoNascI7yoZJg61gL+npsauB9cA5ZtZrZje4+zFgMfAj4Dngu+7ek2Q7pMXysGdFJQoWEoi6YdJgK9hLNFC4+0J3P8Pdh7t7h7vfFxx/zN3f6+7vcfelSbZBUpKHPSsqUfW2lIij/mIwBIzMV2ZLTuU1sQ2q3paTRMlfQP5XqFWgkOTkMbFdpOptqaBdV6hVoJDk5GnPikpUvS0VtOMKtQoUkqwwe1Zk+f8XBQupIuoOe+u2r8vNcJQChSSvXj1bVhPbRR/vgSG1ol1Ae2+3pXYo2FOgkOTlObFdFLZ6++dZ74gkJcoMqawX7ClQSGvkObEN4au3+9/UEFSbi6PCO2v5CwUKaY28J7YhfI3FtuVwaGfy7ZFMi1LhnbWCPQUKaZ28J7YhfLA4vFtXFhK5wnv5xuVsenVT6gFDgUJaK++JbSgEi2Eh9gdQclsCeS/YU6CQ1hoMiW2AC+slXQJamlxK5LVgT4FCWi9MYntfKxoSgZYmlyblsWBPgUJaL0xie0crGhKRVpuVCOIo2GtVwluBQtJRL7HdD+QhF6xgIRHFkfAedvewRIejFCgkPfUS28vJfmIbtDS5xCJKwV6f9yU6HKVAIekZLIlt0NLkEpsoM6SSGo5SoJB05b1iu1QjS5Nr2qzUUZwhZVjD712+cXmswUKBQtIV5qrii61oSEzCrja7MU+dkrR0zexizhlzmkp4r9y0MrZ2KFBI+updVextSSvi8/EeGFpntdmjeeuUpKmZhHef98X2+QoUkr4uoLkJH9k17txwS5OLNKCYv5g+sf5V61BrfMptNQoUkg3LqD4ENaGVDYlRraXJR+S1U5IFPbf01B2OWnTBotg+T4FCsmMFUF6oOgL4RgptiUO16u0hI+CCvHZKsqLacJRh3DT3JpZ9dFlsnzUstjOJRNUV/FwS/Hw3sLTkeB5dtAwmXQJblhSWHh89BWYtLQQRkRgs++iyWINCJQoUki1dwa0beDHVlsRnapcCg+Sahp5ERKQmBQoREalJgUJERGpSoBARkZoUKEREpCZz97TbEDsz203zW99MBPbE2Jw0qA/py3v7QX3Iilb24d3uPqn84KAMFFGY2UZ3n5t2O6JQH9KX9/aD+pAVWeiDhp5ERKQmBQoREalJgeJk8S3inh71IX15bz+oD1mReh+UoxARkZp0RSEiIjUpUIiISE0KFAEzu9LMnjezbWZ2e9rtqcbM7jezXWb2bMmx8Wb2YzN7Ifj5zpLn7gj69LyZXZFOqwcyszPN7Ekze87Meszsi8Hx3PTDzEaa2QYz2xL04a7geG76AGBmQ83sF2a2Nnict/a/aGbPmNlmM9sYHMtbH04zs4fM7FfB/xMXZ64P7t72N2Ao8G/AWRS2ytkCTE+7XVXa+iFgDvBsybE/A24P7t8O/Glwf3rQl3cAU4M+Ds1AH84A5gT3xwL/GrQ1N/0ADBgT3B8O/Bz4YJ76ELTr94HvAGtz+rv0IjCx7Fje+vC3wOeD+yOA07LWB11RFFwEbHP3f3f3I8Aa4BMpt6kid/8JsK/s8Cco/LIR/Ly65Pgad3/b3bcD2yj0NVXu/qq7Px3cPwA8B7yLHPXDCw4GD4cHNydHfTCzDuCjQOmerblpfw256YOZnUrhj7/7ANz9iLu/Tsb6oEBR8C7gpZLHvcGxvJjs7q9C4UsYOD04nvl+mVkncD6Fv8hz1Y9g2GYzsAv4sbvnrQ9/BXwZ6C85lqf2QyE4/4OZbTKz4ibReerDWcBu4IFgCPBeMzuFjPVBgaLAKhwbDPOGM90vMxsDfA/4XXd/o9ZLKxxLvR/u3ufus4EO4CIzm1Hj5Znqg5l9DNjl7pvCvqXCsdT/DYBL3H0OMB+4xcw+VOO1WezDMApDycvd/XzgTQpDTdWk0gcFioJe4MySxx3AKym1pRn/YWZnAAQ/dwXHM9svMxtOIUiscveHg8O56wdAMFTQDVxJfvpwCXCVmb1IYaj1N83s2+Sn/QC4+yvBz13A31EYhslTH3qB3uBqFOAhCoEjU31QoCh4CjjbzKaa2QjgGuDRlNvUiEeB64L71wGPlBy/xszeYWZTgbOBDSm0bwAzMwpjss+5+1+WPJWbfpjZJDM7Lbg/Cvgw8Cty0gd3v8PdO9y9k8Lv+z+6+7XkpP0AZnaKmY0t3gc+AjxLjvrg7q8BL5nZOcGhy4GtZK0PaWf8s3IDFlCYffNvwJK021OjnauBV4GjFP66uAGYAKwDXgh+ji95/ZKgT88D89Nuf9Cm36BwufxLYHNwW5CnfgDnAb8I+vAs8NXgeG76UNKueZyY9ZSb9lMY398S3HqK/9/mqQ9Bm2YDG4Pfpe8D78xaH7SEh4iI1KShJxERqUmBQkREalKgEBGRmhQoRESkJgUKkSaYWZeZTUm7HSKtoEAhUsbMDgY/O83styo8fwMwyd13NnHuPyh7/C9NN1SkRTQ9VqSMmR109zFmNg/4krt/rIH3DnX3vnrnjqGZIi2jKwqR6r4OXBrsdfB7wSKAf25mT5nZL83sCwBmNs8K+2t8B3gmOPb9YKG6nuJidWb2dWBUcL5VwbHi1YsF53422F/hMyXn7i7Zr2BVUNmOmX3dzLYGbfmLlv/XkbYxLO0GiGTY7ZRcUQRf+Pvd/UIzewfwUzP7h+C1FwEzvLD0M8D17r4vWN7jKTP7nrvfbmaLvbCQYLlPUqjQnQVMDN7zk+C584FzKazp81PgEjPbCvw34H3u7sXlRESSoCsKkfA+Anw2WFr85xSWWTg7eG5DSZAAuM3MtgA/o7CI29nU9hvAai+sSPsfwP8DLiw5d6+791NY7qQTeAM4DNxrZp8EDkXsm0hVChQi4Rlwq7vPDm5T3b14RfHm8RcVchsfBi5291kU1oQaGeLc1bxdcr8PGObuxyhcxXyPwqY2jzfQD5GGKFCIVHeAwlatRT8CbgqWSMfM3husWlpuHPCf7n7IzN5HYYvUoqPF95f5CfCZIA8yicKuZ1VXBQ328hjn7o8Bv0th2EokEcpRiFT3S+BYMIT0IPANCsM+TwcJ5d2c2KKy1OPAjWb2SworfP6s5LmVwC/N7Gl37yo5/nfAxRRWQnXgy+7+WhBoKhkLPGJmIylcjfxeUz0UCUHTY0VEpCYNPYmISE0KFCIiUpMChYiI1KRAISIiNSlQiIhITQoUIiJSkwKFiIjU9P8BkZ8aMYLUM/kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def grad_pas_fixe_syst_lin(A , b , rho = 1 , Max_iter = 1000 , tol = 1e-8):\n", " \"\"\"Réalise la méthode du gradient à pas fixe pour résoudre un système linéaire de la forme Ax = b, où A est une matrice *symétrique*\n", " Entrées:\n", " - A: Array de taille (d,d), symétrique\n", " - b: Array de taille (d,)\n", " - rho: Flottant - Pas de la méthode du gradient à pas fixe. Defaut: 1\n", " - Max_iter: Int - Nomnbre maximal d'itérations et critère d'arrêt. Defaut: 1000\n", " - tol: Tolérance, variation du gradient maximale et critère d'arrêt. Defaut: 1e-8\"\"\"\n", " \n", " d = A.shape[0] # Dimension du problème\n", " \n", " x_k = np.random.uniform(size=(d,))\n", " k = 0\n", " \n", " X = [x_k] # Historique des itérations x_0,x_1,...\n", " \n", " while (k < Max_iter) and (np.linalg.norm(A@x_k-b) > tol):\n", " k+=1\n", " #print(k,np.linalg.norm(A@x_k - b))\n", " x_k = x_k - rho*(A@x_k-b)\n", " X.append(x_k)\n", " \n", " if (np.linalg.norm(A@x_k - b) > tol):\n", " raise RuntimeError(f\"La méthode n'a pas convergé après {Max_iter} itérations.\")\n", " else:\n", " return X,x_k\n", " \n", "A = np.array([[1,0],[0,1/20]])\n", "b = np.array([1,1/10])\n", "\n", "X_1 , sol_1 = grad_pas_fixe_syst_lin(A , b , rho = 0.5 , Max_iter = 1000 , tol = 1e-8)\n", "X_2 , sol_2 = grad_pas_fixe_syst_lin(A , b , rho = 1.0 , Max_iter = 1000 , tol = 1e-8)\n", "X_3 , sol_3 = grad_pas_fixe_syst_lin(A , b , rho = 1.5 , Max_iter = 1000 , tol = 1e-8)\n", "\n", "print(\"Solution - rho = 0.5:\" , sol_1)\n", "print(\"Solution - rho = 1.0:\" , sol_2)\n", "print(\"Solution - rho = 1.5:\" , sol_3)\n", "\n", "Err_1 = np.array([np.linalg.norm(A@x-b) for x in X_1])\n", "Err_2 = np.array([np.linalg.norm(A@x-b) for x in X_2])\n", "Err_3 = np.array([np.linalg.norm(A@x-b) for x in X_3])\n", "\n", "\n", "plt.figure()\n", "plt.scatter(np.arange(0,np.size(Err_1),1),Err_1,label=\"0.5\",color=\"green\")\n", "plt.scatter(np.arange(0,np.size(Err_2),1),Err_2,label=\"1.0\",color=\"orange\")\n", "plt.scatter(np.arange(0,np.size(Err_3),1),Err_3,label=\"1.5\",color=\"magenta\")\n", "plt.grid()\n", "plt.xlabel(\"Itérations\")\n", "plt.ylabel(\"Erreur\")\n", "plt.yscale(\"log\")\n", "plt.legend()\n", "plt.show()\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2 - Application: Laplacien discret 1D\n", "\n", "On se propose d'utiliser cette méthode du gradient à pas fixe résolvant un système linéaire afin de résoudre de manière approchée le problème aux limites suivant:\n", "\n", "$$-u''(x) = f(x) , \\hspace{4mm} x \\in [0,1], \\hspace{4mm} u(0) = u(1) = 0.$$\n", "\n", "Pour cela, on discrétise l'intervalle $[0,1]$ en $(x_j)_{0 \\leqslant j \\leqslant J}$ tels que, pour tout $0 \\leqslant j \\leqslant J$, $x_j=\\frac{j}{J}$ ($x_0=0$, $x_1 = \\frac{1}{J}$,..., $x_J=1$). On note $u_j$ la solution approchée de $u$ au point $x_j$: $u_j \\approx u(x_j)$. On résout donc le système:\n", "\n", "$$J^2\\begin{bmatrix}\n", "2 & -1 & 0 & \\cdots & 0 \\\\\n", "-1 & 2 & -1 & \\ddots & \\vdots \\\\\n", "0 & \\ddots & \\ddots & \\ddots & 0 \\\\\n", "\\vdots & \\ddots & -1 & 2 & -1 \\\\\n", "0 & \\cdots & 0 & -1 & 2\n", "\\end{bmatrix} \\begin{bmatrix}\n", "u_1 \\\\\n", "u_2 \\\\\n", "\\vdots \\\\\n", "\\vdots \\\\\n", "u_{J-1}\n", "\\end{bmatrix} = \\begin{bmatrix}\n", "f(x_1) \\\\\n", "f(x_2) \\\\\n", "\\vdots \\\\\n", "\\vdots \\\\\n", "f(x_{J-1})\n", "\\end{bmatrix},$$\n", "\n", "et, par hypothèse, on a automatiquement $u_0 = u_J = 0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q2-** Résoudre cette équation à l'aide de la méthode du gradient à pas fixe. On prendra $f(x) = \\pi^2\\sin(\\pi x)$, $J=100$, $\\rho = 10^{-5}$, $\\verb|Max_iter| = 10^6$ ,$\\verb|tol| = 10^{-5}$. Tracer la solution exacte, donnée par $u_{ex}(x) = \\sin(\\pi x)$. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxvklEQVR4nO3deZzN5f//8cdrmMlOxlIZMZUlstTYd2WPfHyKUAoJZa9kl2whWWJkC6kkSoV8kmRSCVEiDZpkGcvYCoMx2/X7Y+Y7v0kzzIxz3tc5c1732+3cbvOec53zfl7OeL/Oe7suMcaglFJKAfjZDqCUUspzaFFQSimVQouCUkqpFFoUlFJKpdCioJRSKkVO2wFuRpEiRUzp0qWz9NpLly6RN29e1wbycNpn36B99g030+edO3eeMcYUTes5ry4KpUuXZseOHVl6bVhYGI0aNXJtIA+nffYN2mffcDN9FpHD6T2nh4+UUkql0KKglFIqhRYFpZRSKbz6nIJSSrlKXFwckZGRxMTE2I6SIQULFiQ8PPy6bXLlykVQUBD+/v4Zfl8tCkopBURGRpI/f35Kly6NiNiOc0MXL14kf/786T5vjOHs2bNERkYSHByc4fd15PCRiCwSkVMi8ms6z4uIvCkiESKyW0QecCKXUkr9n5iYGAIDA72iIGSEiBAYGJjpPR+nziksAVpc5/mWQJnkR0/gLQcyKaXUP2SXgvB/stIfRw4fGWM2i0jp6zRpCyw1SeN4bxWRQiJyuzHmhBP5lMoKk5AAx44hJ06wd9s2tu7YweXoaGKuXOHq1av4+/vTp2FD8t16K3/Ex3M6IIDb772XoBo1yJErl+34SqVJnJpPIbkorDXG3JfGc2uBScaY75KXNwJDjDH/ujNNRHqStDdB8eLFQ5YvX56lPNHR0eTLly9Lr/VW2uesS4iJ4eQ33/D7li38duAAEX/9xe9Xr7IXuBOYAgxJ43VngEBgKDA5+XcBwD3+/txTsCATWrUivkoVLpQvT2KePDedE/RzzqqCBQtyzz33uCiR+yUkJJAjR44btouIiOD8+fP/+F3jxo13GmOqpdXeU040p7WPk2a1MsbMB+YDVKtWzWT1jj69A9I33EyfD4WFUeCHHyi8eTPvb9zIk3FxABT386PSrbdS9957kQ4doHJlnsmdm8fz5iVfsWLkKliQgHz5iLtyhdwAV67Qe9cuGuzYwfE//uBAeDjhhw5x6vx5Ki5digDdgD/y56dlzZq06t6dyh06IBn4D+/qPnsrV/Q5PDz8uiduPc2NTjT/n1y5cnH//fdn+H09pShEAiVTLQcBxy1lUT7s2PbtfDh2LMu+/pqdV64wC+hbtixNnn6aZXnyUKdTJ+6sWfNfx2oDkx+p+QcEJP1QsCClW7SgdIs0TqudOwfbtnHfG2+wa8sWhn/1FcO/+oq7nnqKvg0aMGj6dKhUCbLZsW6Vvlq1arF8+XJKly7NsWPHaNu2bZaH88kKTykKq4G+IrIcqAmc1/MJyjHx8cR++in/fe451p05gwFC8ubl9UceofULL0DDhhQHOrlj3YULQ8uWvNiyJS8Cx3fvZt2MGaxYs4ZjmzZBlSokVq3KiurVaTtxIrmLFHFHCnWtgQNh1y7XvmfVqjBjxnWbGGM4cuQIpUqVAmD37t1UqlTJtTluwKlLUj8AfgDKiUikiDwjIr1FpHdyk3XAQSACWAA870Qu5dsuRUWxvmdPuOsuAtq3p9ClS4xo0ID9X33FjuhoXvrsM0o3bOhopjsqV6bHokV8efo0r588CbNn8+3583RasIASxYoxrGZNonbvdjSTck5ERATBwcEpe6KZKQpdu3Z1SQanrj667pes5KuO+jiRRanzhw8zq1s3ZoSF8bcxHKlThztmzuS9Nm0gp6fsPIMUKwZ9+tDguefYFBrK7ClTmLx9OzOqVKFHxYqMf/99ClapYjtm9nSDb/TusmfPnn8UgR07dtCrVy+OHj3K2LFjKViwIA0bNmTTpk1cvHiR4sWLM378eC5fvuyyocN17CPlM2L++otpbdpwV3AwozZtolaxYnwzbx53fP89tGvnUQUhNfHzo1G/fnx09Cj7NmzgifLl+WLvXnJXqwYDB5IYFWU7onKRc+fOkTt3biDpxPfnn39OpUqV2LdvHwEBAfTv35/w8HA6d+7MtGnT2LdvHwA//fQTDzzgmnt+tSio7M8YeO89TlaowNC1awkJDGTHsmWsPXmSuj172k6XKWWbNGFheDh7Dx4koFs3rrz5JlXuuIPXW7Yk9uJF2/HUTWrevDkbN26kQ4cOrFy5ksDAQIoXL07Tpk3p168fffv2ZcOGDVSqVInY2FjyJF/G/OOPP1K9enWXZPDMr0ZKuUjUd9/x6tNP88qRI5SuVo29r79OmSeftB3rpgUEB8P8+Zzv0oVS7dvz8hdfsCgwkFmvvkrO2rVtx1NZVLJkSXanOmc0evRoAIYMGUJCQgJ33nknTZs2pWfPnvj7+zNs2DAA9u7dS//+/V0TwhjjtY+QkBCTVZs2bcrya72VL/X5yrlzZmjNmiYnmFtFTOSUKcYkJNiO5TZrx441d/n7G8C0K1rURP/xh+1IjnLF3/Zvv/1280EcdOHChQy1S6tfwA6TznZVDx+pbOeHefO4/7bbmLRtG48GBXFg/35KDB4Mftn3z/3hUaPYe/o0rz70ENFnziSdb3j//aRDZ0plQvb9X6J8z9WrXBkwgHa9e3M5IYH1kybR+913KVKmjO1kjshVsCCjv/qKsW+/jV+5cpx48km6lSrFud9/tx1NeREtCipbOPz11yTWqkXuN99k9X//y54jR2g2JK3RiLK/mOBg+O47tj79NO8fPUrl8uUJe+MN27GUl9CioLzeez17UvGhh5hx4AB89hk1Pv6YAnfcYTuWXTly0G7JEn547z3y5MjBgy+9xKi6dUm4etV2MuXhtCgor3X1/Hmer1CBLgsWEFKwIB3CwuCRR2zH8ighTzzBT0eP8nSZMozfsoUJ5cvDqVO2YykPpkVBeaXIbdtoWKIEb4WHM7hmTTaePEmQi67Tzm7yFS/O4gMHeL9nT/qfOAEPPID54QfbsZSH0qKgvM+WLfzRvDm/X77MRy+/zJStW8mpk9bcUOd58yj0ww9cyZmTRnXrsuy552xHUh5Ii4LyKuFTpkDjxjQsUoTDP/7Io5Mn3/hF6v+7/36ubNwIBQrwxNy5jK5Xj8T4eNuplAfRoqC8gklM5NWGDak4ZAhflS8P27aRLyTEdiyvVPjuu9kQGUm3smUZ9/33dAoOJuavv2zHUh5Ci4LyePExMfS6917GbN7MU/fcQ/3vvoPAa6e0UZkRkC8fb4eHM6V1a1ZERtK7fHnQwuARatWqxaFDhwA4duwY1aqlOWum2+jYR8qjXT59mk6VKrE6KooR9eszLiwMycZ3JjtJ/PwYvGYNpQYNIiQ0FOrVgy++gJIlb/xiH5DW9J4dOnTg+eef5/Lly7Rq1epfz3ft2pWuXbty5swZHnvssX88FxYWdsN1Gl+ZZEepLPnrLz6tU4c1UVGEduzI+M2btSC4QYfp07l7/XrM0aMMrViR39evtx3JZ93MJDuuonsKyiOZU6eQ5s3pfOQIVWfOpIKrRoBUaWvcmGMffcTbLVqwpFUrvlqxgvsefdR2Kquu980+T548132+SJEiGdozuFZ6k+zExcXxyiuvcPnyZWJjY3n11VcZNGgQkyZNYvTo0cybNw9/f/9Mry8t+rVLeZwTP/1EvdKl2REeDqtXa0FwSFCzZmxeuxY/oFH79uxcutR2JJ+T3iQ78+fP58qVKxQqVIjo6GiKFi1KUFAQL774Im+++abLCgLonoLyMMd37KBRnTocj4sjesYMaN7cdiSfcm+rVny7cSMPNW3Kg08/zYa4OGo884ztWD6jefPmvPnmm3To0IH77rsvZZKdn3/+mdDQUG655RYAoqOjOXToEDlz5iRfvnwuzaBFQXmMEz/9ROM6dTgRF8eX8+ZRx8tmRcsu7m7UiG+//56W9erxV//+UKkS1KhhO5ZPSG+SnbZt29K1a1dKlixJgwYNWLVqFaNGjWLdunWEhYWleVI8q7QoKI9w+tdfebB2bY7FxfHFnDlaECwrWaMGu/bvJ2eTJtCsGec//ZSCLtzwqMxp06YNbdq0SVlu3bo1Fy9eZPDgwS5fl55TUPadPk3Bxx4jJCGBdbNmUU+HX/AIOYODYdMmlvv7c8+DD7Lno49sR1IO0KKgrLpw9ChnH3qIgMOHeW/jRhr07Ws7kkrtzjupuWoVt/j50ezxx4nYsMF2IuVmWhSUNVfOnuWR++6jyZ49xK9YAQ0b2o6k0hBcvz4bPvuMOGNo0rIlkdu22Y6k3EiLgrIi7vJlOtx7L5svXODl/v3Jmep4qfI89z78MOvffZdzCQk0bdCA6D//tB3JLUw2m9M6K/3RoqAcZxIT6V2lCmtPn2ZO5850mjnTdiSVASFPPMHaWbPolphI3scfh0uXbEdyqVy5cnH27NlsUxiMMZw9e5ZcmRxWXq8+Uo6b1qwZiyIiGN2wIb3ff992HJUJDfr2pUFQEDz6KAfbtOHOdeuyzVwWQUFBREZGcvr0adtRMiQmJuaGG/xcuXIRFBSUqffVoqCc9dZbdNq4kZhq1Rj+9de206is+M9/ODFxIiFDh9KhalXm/vZbthiTyt/fn+DgYNsxMiwsLIz777/f5e/r/Z+k8hq/hoaS0KcPd7Rpw4gffsgWGxJfdfuQITxXpw7z9+9nShqjhSrvpf8rlSN2r1xJ7b59GVK0KHzwAeTUnVRvN+Hbb+lUqhRD16/noxdesB1HuYgWBeV2J3ftonWnThTw82PQunWQN6/tSMoFxM+PRb/8Qu38+ekyfTq/vPuu7UjKBRwrCiLSQkT2i0iEiAxN4/mCIrJGRH4Rkb0i0s2pbMp9rpw9yyN163I2IYG1y5ZRQqfQzFZyFSzIZ9u306dgQcq/+CIcPmw7krpJjhQFEckBhAItgQpAJxGpcE2zPsBvxpgqQCPgDREJcCKfchNj6FO9OjsuX2bZ8OHc//jjthMpNyhavjxTt23jlthY/m7dmstecvWOSptTewo1gAhjzEFjTCywHGh7TRsD5JekKYfyAeeAeIfyKXeYNInef/7Jm23a0HbCBNtplDuVK0fM0qXU+fVXuoeEYBITbSdSWSRO3KghIo8BLYwxPZKXuwA1jTF9U7XJD6wGygP5gceNMZ+n8V49gZ4AxYsXD1m+fHmWMkVHR7t8HHJP52SfE7/4gsZTpnCqcWPCR46E5OkFnaafs7PWDR3K69u28VJICA9PnerYevVzzpzGjRvvNMZUS/NJY4zbH0B7YGGq5S7ArGvaPAZMBwS4B/gTKHC99w0JCTFZtWnTpiy/1ls51ef969aZgmDeDAoy5tIlR9aZHv2cnZWYkGA6lS5tBMzno0c7tl79nDMH2GHS2a46dfgoEiiZajkIOH5Nm27AquTMESQVhfIO5VMuEn3iBP9t146cIjyyYgXkyWM7knKQ+PmxcMcOquTOzRNjx/LHxo22I6lMcqoo/AiUEZHg5JPHHUk6VJTaEeAhABEpDpQDDjqUT7mASUykR82ahF+9yvLJkylVu7btSMqCPIGBrFq/nrr+/uTq0yfbjZGU3TlSFIwx8UBfYD0QDqwwxuwVkd4i0ju52TigjojsATYCQ4wxZ5zIp1xjZrt2fHj0KONbtKCJG2aEUt4juH591q5dS4kDB0js0UNPPHsRx24rNcasA9Zd87u5qX4+DjRzKo9ysW+/pciaNXQuWZKhn//r+gDli5o1I3rkSNqPG0ebhASeX7HCdiKVAXpHs7ppJioKOnbkyXvu4f1ff9UxjVSKPKNHI0WLMmjlSnbqHc9eQf/3qpuSGBfHY5Urs/T0aVi5EgoUsB1JeRC/nDlZ+t13FMuRgw7du3Ne73j2eFoU1E2Z1LIlq06d4nKHDlCliu04ygMVKVuW5aGhHI6Pp3vt2np+wcNpUVBZ9u2sWYzauJFOpUrRa+lS23GUB6vbqxeTWrdm64kTHBs3znYcdR1aFFSWnIuIoPOgQQT7+zNP50ZQGfDiZ5+xt3lzgiZOhJ9/th1HpUP/J6vMM4bV7dtzKiGB5YsWkf/2220nUl5A/Pwo9N57xBUuzBvNmnEpKsp2JJUGLQoq8+bMoeuuXewfOZJqTz5pO43yJkWKsGPECAafOcPABg1sp1Fp0KKgMuW3Tz9l+6BB0LIlpV991XYc5YVq9+3LkDp1WHjgACsHDbIdR11Di4LKsKvnz9O5UyfaJSRwdd480PMIKovGfvUVNfLmpdfMmURu3247jkpF/1erDBvVpAm/xMQwb9QobilZ8sYvUCod/rlz896nn3LVGHq2bAl6marH0KKgMmTTG28wdccOelesSOsxY2zHUdlAmSZNeOe555hw7hxMm2Y7jkrm2NhHynv9fegQTw0ZQtmAAN4IC7MdR2Ujj4WGwsmTMHw4l+rVI2+tWrYj+TzdU1A3lG/YMHobw3tvv02eIkVsx1HZiQjMn8+QgADqNm5M7MWLthP5PC0K6rrMypXkXL6cEWPG6OWnyj2KFKHeCy/wS0wMY1u0sJ3G52lRUOmK2rOHKp06EVauHAwbZjuOysbajB1LtzJleG3LFrYtXGg7jk/ToqDSZBIT6dW0KQcSEig2dSrk1NNPyr2mb9xIUI4cPPX881w+fdp2HJ+lRUGl6b3evfksKooJjzxChdatbcdRPqBgyZIsnjyZqLg4dj3/vO04Pku//ql/OfHzz/RfuJC6BQow8KOPbMdRPuTBF1/k8IEDFJw/H779FurXtx3J5+iegvonY1jSuTMxxrBo5Upy+PvbTqR8TMFp0zClS/Ne+/ZcPqPTtDtNi4L6p/ffZ+i+ffw0ZAhlm+mU2cqCvHnZOWQIXaKiGK1XIzlOi4JKEbV7Nwf79EFq1+beCRNsx1E+rFrv3vSqUIFpO3fyw7x5tuP4FC0KKkXfVq2ofuEC0bNmQY4ctuMoHzflyy8pmSMHz/Tvz9Xz523H8RlaFBQAnw4dykfHjvFC06bkCwmxHUcpCpQowdwxYwiPjWVimza24/gMLQqK84cP0+f116mcKxcvf/qp7ThKpWg5ciQvV6xI7e+/h927bcfxCVoUFC+3aMHJxETenj8f/zx5bMdR6h8mf/MNLQIDoUcPSEiwHSfb06Lg4xK//poc+/bxYvXqVOvSxXYcpf4tMJCEGTMY9eOPzG7f3naabE+Lgi+LicGvd2/m3HUXkzdtsp1GqXT5dezIT8WKMfSTTziyZYvtONmaFgUftuCxx9j6++8wdy6SN6/tOEqlS/z8CP34YwzQt107jM7U5jZaFHzUvrVr6fv558wuXRqaNrUdR6kbKl2vHq+2bs2aU6dYNXiw7TjZlhYFH5QYH0+vJ54grwjT1q61HUepDBu4ciVVc+fmhRkziNMhMNzCsaIgIi1EZL+IRIjI0HTaNBKRXSKyV0S+cSqbr1ny7LNsvnCBKU89RbGKFW3HUSrDcubKxZK332a1MfiPHm07TrbkyCipIpIDCAWaApHAjyKy2hjzW6o2hYA5QAtjzBERKeZENl9z9sABXn7nHeoWKED3t9+2HUepTKvSqRNs3QqzZnG5Y0fyNGhgO1K24tSeQg0gwhhz0BgTCywH2l7TpjOwyhhzBMAYc8qhbD6lwGuvMViEt5YswU+HslDeatw4XsyTh4YtW5IQG2s7Tbbi1HwKJYCjqZYjgZrXtCkL+ItIGJAfmGmMWXrtG4lIT6AnQPHixQkLC8tSoOjo6Cy/1lv5b9+O/5IldOrYkYO33uoT/ffFz9lX+nxbs2bs+OQTxjZrRvWXXvKJPqfmts/ZGOP2B9AeWJhquQsw65o2s4GtQF6gCPA7UPZ67xsSEmKyatOmTVl+rTeKvXTJNM2d26wuUsSYixdtx3GMr33OxvhOnxMTEkyTwoVNATCrFyywHcdxN/M5AztMOttVpw4fRQIlUy0HAcfTaPOFMeaSMeYMsBmo4lC+bG9Wx45suHKF+GeegXz5bMdR6qaJnx+hH3xADLBk5EjbcbINp4rCj0AZEQkWkQCgI7D6mjafAfVFJKeI5CHp8FK4Q/myteM//cQra9bQpFAh/jNxou04SrlM2WbNeLl+fTZGRRG1apXtONmCI0XBGBMP9AXWk7ShX2GM2SsivUWkd3KbcOALYDewnaTDTb86kS+7G/zf/xIL9Bw5EvHTW1NU9jLsk0/4uWhRio8eDXFxtuN4Pce2EMaYdcaYssaYu40xE5J/N9cYMzdVm9eNMRWMMfcZY2Y4lS07+3nBApYdPsyQ+vUpqvMkqGwoT2Ag0QMGkLh3L+GjRtmO4/X0a2N2FhdH1enTWVW0KEM/+cR2GqXc5mydOgy76y5qTp7MiV27bMfxaloUsrGr06cj4eG0W7iQPIGBtuMo5T4i9Jg3j6vAy+3a2U7j1bQoZFNRe/Zw19ChrKhaFXQqQ+UDyjRpwkt16/LeoUN8/9ZbtuN4LS0K2dTwRx/llDFUee01ELEdRylHDF+1ihI5ctD/xRf1Tucs0qKQDW1fvJhFv//OwOrVKdeihe04Sjkmb7FivP7885y7coUjU6fajuOVtChkM4nx8fTv14/b/PwYpSeXlQ/qOGMG4fXqETxtGvz1l+04XifTRUFE8iaPeqo80JbRo9l26RKTe/SgQIkStuMo5Tjx8yPX7NlcOXeOL55+2nYcr3PDoiAifiLSWUQ+F5FTwD7gRPKcB6+LSBn3x1QZcvEi9RYvZnflyjwZGmo7jVL2VKnCuKpVab1mDb+tvnbwBHU9GdlT2ATcDQwDbjPGlDTGFAPqkzSA3SQRedKNGVUGRY0YASdPUmnBAvxyOjUArlKeadCyZeQX4YXu3XVO50zISFFoYowZZ4zZbYxJ+Zc1xpwzxnxsjHkU+NB9EVVG/BkWRvCsWbxduzbUqGE7jlLWFS1fnlfatmX92bOsGzPGdhyvccOiYIyJAxCRGSJpX9v4f22UPYOfeAIBWsyaZTuKUh6jz/vvUy4ggEGTJhEbHW07jlfIzInmaGC1iOQFEJFmIvK9e2KpzAibPp2Pjx9n2EMPUULHN1IqhX+ePEwbMYKicXGcnjLFdhyvkOGiYIwZCXwAhInId8CLwFB3BVMZkxAby6Dhw7kzRw5eXLHCdhylPE6r0aP5rkULSsycCadP247j8TJcFETkIeBZ4BJQFOhvjPnWXcFUxuyZMIEDMTFM7tOH3IUL246jlEeSadM4GR3NB088YTuKx8vM4aMRwChjTCPgMeBDEXnQLalUxly8SNX584kICeHx6dNtp1HKc917L6/ddx9dNmwgfM0a22k8WmYOHz1ojPku+ec9QEtgvLuCqRvb//LLmJMnuT00VCfPUeoGRn74IXmBl555xnYUj5aRm9fSu+LoBPDQ9doo9zm6dSv3z53LhMqVoWZN23GU8nhFy5dn5MMPs+70ab7UaWnTlaGb10Skn4jcmfqXyXMt1xaRdwC9l9xhwzt1IhHookMEK5Vh/ZctIzhnTl4cO5aEq1dtx/FIGSkKLYAE4AMROSEiv4nIn8DvQCdgujFmiRszqmvsfPdd3jt0iEG1a1OqTh3bcZTyGrcUKMDrAwdS9epVoufNsx3HI91wLARjTAwwB5gjIv5AEeCKMeZvN2dTaTCJibzUrx9FRRiml6AqlWmPTpnCo1u2wGuvQffukC+f7UgeJTOXpP4IzAUeBx4QkSJuS6XSdeKddzhw/jyvdOhAgaAg23GU8j4i8MYb7Dp5kmU6iuq/ZOaSlbbASiAA6A0cFpHDbkml0hYXxx2TJ3OgbFl6LlpkO41S3qtWLV4rWZJnV63ixM8/207jUTJzSepxY8wXxpgpxpgOQDVgofuiqWvtHjOG2P37yTt1Kv558tiOo5RXm/jOO8QBozp2tB3Fo2Tm8NE/rj4yxoQDFV2eSKXp4rFjNHntNZ4pVgxat7YdRymvd3fjxvS5/34WHzjA3k8/tR3HY2Tm8NGHIhIpIt+KyBwRmQaUd1cw9U9TO3XitDH0mzQp6ZioUuqmjfzwQ/IDQ3v1sh3FY2Tm8FFtY0wQ0A3YAOwF9CurA07s2sXUb7+lQ8mS1OjWzXYcpbKNwDJlGPvII1Q9dYqEr76yHccjZHp6LmNMBBDhhiwqHWM6dSIOmLhkie0oSmU7/Zcvh3LlYNgwePBB8PEhY3y7914gbvdudu3bx3NVqnD3gzr+oFIulzs3ZuxYVu/YwRfDh9tOY50WBQ/nP3o0P+TPzySdfFwpt0ns3JmRuXLR5403fH6GNi0KHiziww85+9ln+A0ZQu4777zxC5RSWZIjIIDJw4dzMD6e+T5+3k6LgocyiYl069mT+jlzYgYMsB1HqWyvxYgRNCpUiHEff8zF48dtx7HGsaIgIi1EZL+IRIhIutN4ikh1EUkQkcecyuaJ1r7yCt9duED/Dh0QHZtFKbcTPz8mTZvGKWOY1rmz7TjWOFIURCQHEErSxDwVgE4iUiGddpOB9U7k8lQJsbEMe/11yvj788yCBbbjKOUzanbrRt+77+burVshKsp2HCuc2lOoAUQYYw4aY2KB5SSNpXStfsDHwCmHcnmkd597jr1XrzKxXz8dzkIph81at44n4+NhvG9OLJnp+xSyqARwNNVyJPCP6cJEpATQDngQqJ7eG4lIT6AnQPHixQkLC8tSoOjo6Cy/1p38YmPZ+sEHPJArF4EPP+zSjJ7aZ3fSPvsGV/f5zhYtWBsaStDdd1O4alWXva8rue1zNsa4/QG0BxamWu4CzLqmzUqgVvLPS4DHbvS+ISEhJqs2bdqU5de61bRpxoC59PnnLn9rj+2zG2mffYOr+3xs506TG0yX4GCXvq8r3UyfgR0mne2qU4ePIoGSqZaDgGtP71cDlovIIeAxkib1+Y8j6TzEhchI9o8dC02bkqdVK9txlPJZdzzwAP1r1eK9P/9kz8cf247jKKeKwo9AGREJTp7buSPwj7uxjDHBxpjSxpjSwEfA88aYTx3K5xHeeOIJ7vv7b4707287ilI+b8iyZRQARvTpYzuKoxwpCsaYeKAvSVcVhQMrjDF7RaS3iPR2IoOnO7V3L9M2b+Y/QUHcqUNjK2XdrcHBDGnenDVRUXw/Z47tOI5x6kQzxph1wLprfjc3nbZdncjkSSY++SRXgPHz59uOopRK1v/dd/m+ZElyzpsHzz3nE8PW6x3NHuDw99/z1q5ddC1blnItW9qOo5RKlrdoUdbOmEHN3bvhf/+zHccRWhQ8wA8jRpALeGXpUttRlFLX6t6dM6VKMfXZZ0mMj7edxu20KNi2bx8dv/2WyD59KFmz5o3bK6WcFRDAFy1bMvj4cT4ePNh2GrfTomDZ3gEDIE8e8r/yiu0oSql0dJo5k4q33MKo0FDiY2Jsx3ErLQoW/bxsGfd9+SWLGjeGokVtx1FKpSNHQADjBw1if1wcS597znYct9KiYNHIgQO5VYRHZ8+2HUUpdQNtJ0yget68jHn3Xa6eP287jttoUbBky9y5rDt9miEtWlBQJ9BRyuOJnx8TR4+mRkIC57PxFzktChaYxERGDB1KcT8/+r7zju04SqkMajJ4MB81akSxWbPg0iXbcdxCi4IFx1es4Nfz5xnerh159VyCUt5DBMaPZ19UFGuef952GrfQouA0YygxfToHg4LotWiR7TRKqcyqW5eXihal67vvcv7IEdtpXE6LgsMiFy8mYft28o8Zwy0FCtiOo5TKglenTeOcMcx46inbUVxOi4KDEuPjad2nD//Jkwey4R+TUr4i5Mkn+e8dd/DGN99w9vffbcdxKS0KDvp48GB+iYmh49NPg7+/7ThKqZswds4cooGp2ewLnhYFhyTExvLKnDlUuOUWOs6YYTuOUuomVWzblm533UWOHTsgKsp2HJfRouCQDwYMIDw2ljF9+pAjIMB2HKWUCyxct47xiYkwebLtKC6jRcEJ8fF8uHQpVXLl4tFs9MejlK+TcuUwTz3FxtmzOb5zp+04LqFFwQlLl/Lp5cusDQ3FL6dj8xoppRxwvGdPWsTFMbFrV9tRXEKLgpvFRkdzccwYclSvTlC3brbjKKVcrETt2jxz773M//VXDn//ve04N02Lgpst7tWL4KNHOdKnj09M5aeULxqxZAkCTOjRw3aUm6ZFwY2unj/P+A8/pGy+fJTs0sV2HKWUm5SsUYNnK1Vi8b59HAwLsx3npmhRcKOFzz5LZEICY0eNQvz0n1qp7GzY4sWUACJee812lJuiWyo3uXLuHBM//pj6BQrw0Esv2Y6jlHKzEiEh/NG/P802bgQvvstZi4KbrBk0iOOJiYx99VXdS1DKR+QYNox4f3829+1rO0qW6dbKHS5fpsP69eysXp1GAwfaTqOUcspttzG5WjUaf/kl+9ats50mS7QouMHVN9+EqCgemDbNdhSllMOenTuXXMA4L51vQYuCi106dYoyI0bwVvnyUK+e7ThKKYcVq1iRvjVq8MHhw4SvWWM7TqZpUXCxOd26cTQxkSp62EgpnzV46VLyAGO98NyCFgUXij55kin/+x/NAgOp06uX7ThKKUuKlCtHv1q12HvkCFe8bEwkLQouFNqtG2eM4dVJk2xHUUpZNnrlSnblzUtuLxsEU4uCi8SdO8f09etpUaQItbLBre5KqZuTOygIvwEDuLByJcc3brQdJ8McKwoi0kJE9otIhIgMTeP5J0Rkd/Jji4hUcSqbK/jPn8+3xjA9NNR2FKWUh4gfMIDKIgz0ohFUHSkKIpIDCAVaAhWATiJS4ZpmfwINjTGVgXHAfCeyuYK5cAGmTqVMq1aU79DBdhyllIfIWawYXerWZWVkJL+uWmU7ToY4tadQA4gwxhw0xsQCy4G2qRsYY7YYY/5KXtwKBDmU7aZNat+etmfPcnXYMNtRlFIeZtCSJeQHxg4YYDtKhjg140sJ4Giq5Uig5nXaPwP8L60nRKQn0BOgePHihGVxRMLo6Ogsvza1q2fOMHXDBh4oVIgf4uPBg0dIdFWfvYn22Td4ep+frlSJ2Xv2sHzsWG5r0MAl7+m2Phtj3P4A2gMLUy13AWal07YxEA4E3uh9Q0JCTFZt2rQpy69NbWKzZgYw2xYtcsn7uZOr+uxNtM++wdP7fDYiwuQHM75iRZe95830Gdhh0tmuOnX4KBIomWo5CDh+bSMRqQwsBNoaY846lC3LLh4/ztQNG2hVtCg1dFY1pVQ6Ct99N/sGDGDE3r2we7ftONflVFH4ESgjIsEiEgB0BFanbiAidwKrgC7GmAMO5bopc7p355wxvOJl1yErpZx3x+jRkD8/Z0eOtB3luhwpCsaYeKAvsJ6kQ0MrjDF7RaS3iPRObjYaCATmiMguEdnhRLYsi46mx/btLKpUSfcSlFI3VrgwnzVvTok1a/j1k09sp0mXY/cpGGPWGWPKGmPuNsZMSP7dXGPM3OSfexhjbjXGVE1+VHMqW5aEhhL41190W7DAdhKllJeoN3EiAcA4D74SSe9ozoLokydpPmoUP9SsCTWvdxGVUkr9f4FlytCvTh1WHj3K3s8+sx0nTVoUsmBO9+58GReH6HAWSqlMemHJEvIC4/r1sx0lTVoUMunSqVNM/eILmgcG6hhHSqlMCyxThn61a/PJ0aNEbd5sO86/aFHIpLe6d+e0MYyeMMF2FKWUl3rpnXfYlzs3xefOtR3lX7QoZMLlM2d4fd06mhYurPMlKKWyrHCZMgT37w/LlxO3Z4/tOP+gRSETci5axFhjGDdxou0oSikvZ154gQ5+fnR/5BHbUf5Bi0JGXblCwLRp9HrwQWrqXoJS6iZJsWKUuv9+lh06xIH1623HSaFFIYOW9ujBnKgoEkeNsh1FKZVNvLR4MbcA459/3naUFFoUMiDm778Zunw5KwoWxK9RI9txlFLZRPH77qP3Aw/w/sGD/L5hg+04gBaFDFn47LOcSEzkldGjbUdRSmUzLy9eTADwRt++tqMAWhRu6OqFC0xatYr6BQrQaOBA23GUUtnMbZUr80nbtkyKiICDB23H0aJwI4t79eJYYiKjR4xA/PSfSynlei3mzKGQvz+89prtKFoUris2lrIbN9KreHEeeukl22mUUtnVHXfwU9u2VF+4kMPffWc1ihaF61myhAdPn2buO+/oXoJSyq2KvfQSu4FJli951y1dOuIuX2bikCGcuf9+aNbMdhylVDYXVL063StU4O3ffuPo1q3WcmhRSMd7ffsy4u+/2dq6NYjYjqOU8gFD58/HAFN69rSWQYtCGuJjYpjw7rs8kDs3D48ZYzuOUspHlKpbl67lyrFgzx6O79xpJYMWhTR8MGAAf8THM3rgQD2XoJRy1PB58wgVoeiSJVbWr1u8ayTExjJ+8WKq5MrFI+PH246jlPIxwQ0b8szTT+O/cCGcPOn4+rUoXOP8O+9QMS6OUX376l6CUsqO4cMJvXqVKY8/7viqdauXWmIihadPZ1XFijw6ebLtNEopX1WmDFtLleLVzZs5HR7u6Kq1KKSybdIk9oeHw4gRoHsJSimLRsyezRVgWvfujq5Xt3zJEuPj6TluHP8NCMC0b287jlLKx5V/+GE6lCzJ7K1bOffHH46tV4tCsjWjRrE7JoZh3bohOXPajqOUUoycOZNoYEa3bo6tU4sCYBITGTtzJvf4+9NxxgzbcZRSCoD72rXjlfLlafjTT/D3346sU4sC8L+xY/npyhWGP/kkOXPlsh1HKaVSjPngAx66dAlmzXJkfVoUjCFiyRLu9ffnydmzbadRSql/qlqVU82bM2LiRC4cO+b21WlR2LCB/ocPs3vmTPzz5LGdRiml/uVI585MjIkh1IErkXy6KJjERHYMHgxBQeR0+LIvpZTKqGpPPUXLokWZtmEDl06dcuu6fLoohM2YQfXdu1nZrBnccovtOEopla5R48dzxhjmPvOMW9fj00Vh7Lhx3O7nR5s33rAdRSmlrqt2z540KVyY19et4/LZs25bj2NFQURaiMh+EYkQkaFpPC8i8mby87tF5AF35vnjk08I+/tvXn7kEXIVKuTOVSmllEuMHjOGhomJXHjrLbetw5GiICI5gFCgJVAB6CQiFa5p1hIok/zoCbiv18B7S5ZQTISeCxa4czVKKeUy9fv148OGDbntrbfwi411yzqc2lOoAUQYYw4aY2KB5UDba9q0BZaaJFuBQiJyuzvCnFq3jl8vXGBwq1bkKVLEHatQSin3GDWKP48fJ+fKlW55e6fGcygBHE21HAnUzECbEsCJ1I1EpCdJexIUL16csLCwTIcpcOgQ26tU4Y8ePbL0em8VHR3tU/0F7bOv8Kk++/lxW/363FK4sFv67FRRSGuSY5OFNhhj5gPzAapVq2YaNWqU+TSNGhFWoQJNsvJaLxYWFkaW/r28mPbZN/hcnxs3dlufnTp8FAmUTLUcBBzPQhullFJu5FRR+BEoIyLBIhIAdARWX9NmNfBU8lVItYDzxpgT176RUkop93Hk8JExJl5E+gLrgRzAImPMXhHpnfz8XGAd0AqIAC4Dzo0Vq5RSCnDunALGmHUkbfhT/25uqp8N0MepPEoppf7Np+9oVkop9U9aFJRSSqXQoqCUUiqFFgWllFIpJOn8rncSkdPA4Sy+vAhwxoVxvIH22Tdon33DzfS5lDGmaFpPeHVRuBkissMYU812Didpn32D9tk3uKvPevhIKaVUCi0KSimlUvhyUZhvO4AF2mffoH32DW7ps8+eU1BKKfVvvrynoJRS6hpaFJRSSqXI9kVBRFqIyH4RiRCRoWk8LyLyZvLzu0XkARs5XSkDfX4iua+7RWSLiFSxkdOVbtTnVO2qi0iCiDzmZD53yEifRaSRiOwSkb0i8o3TGV0tA3/bBUVkjYj8ktxnrx5tWUQWicgpEfk1neddv/0yxmTbB0nDdP8B3AUEAL8AFa5p0wr4H0kzv9UCttnO7UCf6wC3Jv/c0hf6nKrd1ySN1vuY7dwOfM6FgN+AO5OXi9nO7UCfhwOTk38uCpwDAmxnv4k+NwAeAH5N53mXb7+y+55CDSDCGHPQGBMLLAfaXtOmLbDUJNkKFBKR250O6kI37LMxZosx5q/kxa0kzXLnzTLyOQP0Az4GTjkZzk0y0ufOwCpjzBEAY4y39zsjfTZAfhERIB9JRSHe2ZiuY4zZTFIf0uPy7Vd2LwolgKOpliOTf5fZNt4ks/15hqRvGt7shn0WkRJAO2Au2UNGPueywK0iEiYiO0XkKcfSuUdG+jwbuJekqXz3AAOMMYnOxLPC5dsvxybZsUTS+N211+BmpI03yXB/RKQxSUWhnlsTuV9G+jwDGGKMSUj6Eun1MtLnnEAI8BCQG/hBRLYaYw64O5ybZKTPzYFdwIPA3cAGEfnWGHPBzdlscfn2K7sXhUigZKrlIJK+QWS2jTfJUH9EpDKwEGhpjDnrUDZ3yUifqwHLkwtCEaCViMQbYz51JKHrZfRv+4wx5hJwSUQ2A1UAby0KGelzN2CSSTrgHiEifwLlge3ORHScy7df2f3w0Y9AGREJFpEAoCOw+po2q4Gnks/i1wLOG2NOOB3UhW7YZxG5E1gFdPHib42p3bDPxphgY0xpY0xp4CPgeS8uCJCxv+3PgPoiklNE8gA1gXCHc7pSRvp8hKQ9I0SkOFAOOOhoSme5fPuVrfcUjDHxItIXWE/SlQuLjDF7RaR38vNzSboSpRUQAVwm6ZuG18pgn0cDgcCc5G/O8caLR5jMYJ+zlYz02RgTLiJfALuBRGChMSbNSxu9QQY/53HAEhHZQ9KhlSHGGK8dUltEPgAaAUVEJBJ4BfAH922/dJgLpZRSKbL74SOllFKZoEVBKaVUCi0KSimlUmhRUEoplUKLglJKqRRaFJRSSqXQoqCUUiqFFgWlXEhENolI0+Sfx4vIm7YzKZUZ2fqOZqUseAUYKyLFgPuBRyznUSpT9I5mpVwseYazfEAjY8xF23mUygw9fKSUC4lIJeB24KoWBOWNtCgo5SLJM169T9JsWJdEpLnlSEplmhYFpVwgeWjqVcCLxphwkkbrHGM1lFJZoOcUlFJKpdA9BaWUUim0KCillEqhRUEppVQKLQpKKaVSaFFQSimVQouCUkqpFFoUlFJKpfh/D++O07Nrcq8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "J = 100\n", "\n", "def f(x):\n", " \"\"\"Fonction de l'exemple.\n", " Entrée:\n", " - x: Flottant.\"\"\"\n", " return np.pi**2*np.sin(np.pi*x)\n", "\n", "A = 2*np.diag(np.ones(J-1),0) - np.diag(np.ones(J-2),-1) - np.diag(np.ones(J-2),1)\n", "A = (J)**2*A\n", "b = np.array([f(j/J) for j in range(1,J)])\n", "\n", "U, sol = grad_pas_fixe_syst_lin(A , b , rho = 1e-5 , Max_iter = 10**6 , tol = 1e-5)\n", "\n", "u_sol = np.concatenate((np.array([0.0]),sol,np.array([0.0])),axis=0)\n", "u_ex = np.array([np.sin(np.pi*j/J) for j in range(0,J+1)])\n", "\n", "plt.figure()\n", "plt.plot(np.array([j/J for j in range(0,J+1)]),u_sol , color = \"red\" , label = \"$u_{sol}$\")\n", "plt.plot(np.array([j/J for j in range(0,J+1)]),u_ex , color = \"black\" , linestyle = \"dashed\" , label = \"$u_{ex}$\")\n", "plt.ylabel(\"$u(x)$\")\n", "plt.xlabel(\"$x$\")\n", "plt.legend()\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## II - Régression & moindres carrés\n", "\n", "On peut appliquer l'optimisation afin de résoudre des problèmes aux moindres carrés. ce procédé peut en particulier être utilisé afin d'étudier un problème de régression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1 - Régression polynomiale\n", "\n", "Connaissant une suite de points $(x_1,y_1),(x_2,y_2),\\cdots,(x_n,y_n)$, on cherche à lier les $x_i$'s et les $y_i$'s par une relation polynomiale, autrement dit, s'il existe des réels $a_0,\\cdots,a_d$ tels que, pour tout $1 \\leqslant i \\leqslant n$:\n", "\n", "$$y_i = a_0 + a_1x_i + \\cdots + a_dx_i^d.$$\n", "\n", "Cela conduit donc à un système de la forme:\n", "\n", "$$\\underbrace{\\begin{bmatrix}\n", "1 & x_1 & \\cdots & x_1^d \\\\\n", "1 & x_2 & \\cdots & x_2^d \\\\\n", "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "1 & x_n & \\cdots & x_n^d\n", "\\end{bmatrix}}_{:=X\\in\\mathcal{M}_{n,d+1}(\\mathbb{R})} \\cdot \\underbrace{\\begin{bmatrix}\n", "a_0 \\\\\n", "a_1 \\\\\n", "\\vdots \\\\\n", "a_d\n", "\\end{bmatrix}}_{:=\\alpha\\in\\mathbb{R}^{d+1}} = \\underbrace{\\begin{bmatrix}\n", "y_0 \\\\\n", "y_1 \\\\\n", "\\vdots \\\\\n", "y_n\n", "\\end{bmatrix},}_{:=Y\\in\\mathbb{R}^{n}}$$\n", "\n", "soit alors $X\\alpha = Y$. Lorsque $d \\gg n$ (ce qui est en général vrai puisque nous avons beaucoup de données à traiter), un tel système n'admet en général pas de solution, on dit qu'il est surdéterminé. On recherche donc $\\alpha \\in \\mathbb{R}^{d+1}$ minimisant la quantité:\n", "\n", "$$E(\\alpha) = \\|X\\alpha-Y\\|^2.$$\n", "\n", "On montre facilement que:\n", "\n", "$$\\nabla E(\\alpha) = 2(X^TX\\alpha - X^TY).$$\n", "\n", "on peut donc trouver le $\\alpha$ optimal via la suite:\n", "\n", "$$\\alpha_{k+1} = \\alpha_k - \\rho\\nabla E(\\alpha_k).$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q3 -** En s'inspirant de la fonction créée en question **Q1**, créer une fonction prenant en entrée des listes de données $(x_i)_{1 \\leqslant i \\leqslant n}$ et $(y_i)_{1 \\leqslant i \\leqslant n}$ et réalisant une régression polynomiale entre ces données. pour construire la matrice $X$, on pourra utiliser la fonction $\\verb|np.concatenate|$. Tester sur les séries de données suivantes dans le cas $d=1$:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "x = np.array([0.0 , 1.0 , 2.0 , 3.0 , 4.0])\n", "y = np.array([-2.0 , -1.0 , 0.0 , 1.0 , 2.0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On prendra $\\rho = 10^{-2}$, $\\verb|Max_iter| = 1000$ ,$\\verb|tol| = 10^{-8}$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution: y = 0.9999999988973359 x + -1.999999996856546\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmBUlEQVR4nO3deXhUVbrv8e+bkDAlCAoGEAQ96lEEEQoRB5AIIuCA7YgcwAFEVDx6bbudbmu3Pdnt0IqIyCCoqHFAlLFppIND97WVKDKqoKAgCCJDCCCB5L1/pJQYE5JKVbIrqd/nefZD7dqr9vq5DC+VVbvWNndHRERqv6SgA4iISPVQwRcRSRAq+CIiCUIFX0QkQajgi4gkiDpBBziYpk2betu2bSv12l27dtGwYcPYBooB5YqMckVGuSJTG3Pl5ORscfdmpR5097jdQqGQV1Z2dnalX1uVlCsyyhUZ5YpMbcwFLPIyaqqmdEREEoQKvohIglDBFxFJECr4IiIJQgVfRCRBRF3wzay1mWWb2UozW25mt5TSxsxstJmtNrMlZtY52n5FRGqVRo3ArGjLyTnwuFGjmHURi+vw9wO/dPcPzSwdyDGz+e6+olibfsCx4e1U4MnwnyIiArBzZ2TPV0LU7/DdfaO7fxh+vBNYCRxRotkA4NnwZaLvAY3NrEW0fYuI1DYrSv/KVEyYx3A9fDNrC7wNtHf33GLPzwIecPd3w/sLgDvcfVEp5xgBjADIyMgIZWVlVSpLXl4eaWlplXptVVKuyChXZJQrMvGUa/eifzEhfw6v7/s39xxxO713ND9wMBSq8HkyMzNz3L1LqQfL+kZWpBuQBuQAF5dybDZwZrH9BUCovHPqm7bVR7kio1yRUa6Dm/PZHG/9f3C7D7+lLz7nwT+6w4EtAlT1N23NLAWYBjzv7q+V0mQ90LrYfitgQyz6FhGpyV5c+iL9X+hPWj78axI8+neob3WrpK9YXKVjwCRgpbs/UkazGcDQ8NU63YAd7r4x2r5FRGoid2fL7i0ADDh+AA+e8yAfPZ/GaetLaZyeHrN+Y/EO/wxgCHC2mS0Ob/3NbKSZjQy3mQN8AawGJgA3xqBfEZEaZ+POjVz88sV0m9iN3ft20yClAbeffjt1t+08MIkTCh14nJtb/kkrKOrLMr3og1grp40DN0Xbl4hITeXuTF48mdvm3cbegr3c3/N+UpNTqzVDXK+HLyJSG2zbs43LX72cN794kx5tejDhggkcd9hx1Z5DBV9EpIo1qtuIOkl1ePK8JxkRGkGSBbOqjdbSERGpAiu+XcH5L5zP5l2bSU5KZs6gOYzsMjKwYg8q+CIiMZVfkM8f3v4DnZ7qxHvr3+OTLZ8AUHRBY7A0pSMiEiOLNixi2IxhLNm0hIHtB/JY38c4vOHhQcf6kQq+iEiMPPL/HmHL7i28MfANLvzvC4OO8zMq+CIiUXhr7VtkpGVwfNPjebzf4yQnJdO4XuOgY5VKc/giIpWQuzeXG2bdQM9nevK7t34HwGENDovbYg96hy8iErE5q+Zw/azr2bBzA7d1u437M+8POlKFqOCLiETgxaUvMui1QZzY7ERevexVTm1Vc+7lpIIvIlKOHxY7a9awGQOOH8DDfR5mVNdR1b40QrQ0hy8ichBf537NRS9dxGmTTvtxsbPbTrutxhV7UMEXESmVuzMhZwLtxrZj/ufzuaHLDTWyyBenKR0RkRK27dnGJS9fQvbabHq27cmECyZwzKHHBB0rair4IiIlHFLvEOqn1Gf8+eMZ3nl4XCyLEAua0hERAZZtXka/5/uxKW8TSZbErCtncV3oulpT7CFGBd/MnjazzWa2rIzjPc1sR7E7Yt0bi35FRKKVX5DP7xb+js5PdWbRhkV89t1nQHwsdhZrsZrSmQKMAZ49SJt33P38GPUnIhK1lbkruXn8zSzbvIxBHQbxWN/HaNqgadCxqkxMCr67v21mbWNxLhGR6vLq16+ybc82Zl45k/OPq/3vR63odrMxOFFRwZ/l7u1LOdYTmAasBzYAt7v78jLOMwIYAZCRkRHKysqqVJ68vDzS0tIq9dqqpFyRUa7IKFf5Ptr2EYemHkqbhm3YsH0Dh6QdQsM6DYOO9RPRjFdmZmaOu3cp9aC7x2QD2gLLyjjWCEgLP+4PrKrIOUOhkFdWdnZ2pV9blZQrMsoVGeUq2/Y9233EjBHOb/ErX73S3eMjV2miyQUs8jJqarVcpePuue6eF348B0gxs9o7USYicWXGpzNoN7YdEz+ayO2n3c7ECycGHSkQ1XIdvpk1Bza5u5tZV4quDvquOvoWkcT2/JLnGTx9MB0O78DrV7zOKUecEnSkwMSk4JvZi0BPoKmZrQfuA1IA3H0ccClwg5ntB/YAA8O/eoiIxJy7s3nXZjLSMrj4hIv52+6/ceMpN9b4pRGiFaurdK4s5/gYii7bFBGpUut2rOOG2TewcstKlt6wlAYpDbi1261Bx4oLWlpBRGqFQi9kQs4EfjX/VxR4AX86+0/UTa4bdKy4ooIvIjXe1j1bufili3nry7fodVQvxl8wnqObHB10rLijgi8iNV7jeo1Jr5vOxAsmcm2na2vlsgixoMXTRKRGWrJpCedOPZdv8r4hyZKYeeVMhnUepmJ/ECr4IlKj7N2/l/uy7yM0PsRHGz9i9dbVQUeqMTSlIyI1xnvr32PYjGGs+HYFg08azKPnPsphDQ4LOlaNoYIvIjXG6P+MZufencweNJv+x/YPOk6No4IvInFtwRcLaJHegnbN2jGm/xjqJNWhUd1GQceqkTSHLyJxafv32xk+Yzi9n+vNH97+AwCH1j9UxT4KeocvInHn9U9e58bZN7J512buPONO7j1LN8mLBRV8EYkrPyx21jGjIzOvnEmoZSjoSLWGCr6IBM7d2bRrE83TmnPxCRczes9oRnYZSUpyStDRahXN4YtIoL7a8RXnvXAep086nV35u6ifUp+bT71Zxb4K6B2+iASi0AsZt2gcd7x5B+7On3v9mXp16gUdq1ZTwReRard1z1YuyrqId756h3OOPofxF4ynbeO2Qceq9VTwRaTaNa7XmCb1mzB5wGSu6niV1r+pJjGZwzezp81ss5ktK+O4mdloM1ttZkvMrHMs+hWRONaoEZgVbTk5LG5h9L7K+KZ5GkmWxBsD3+Dqk69Wsa9GsfrQdgrQ9yDH+wHHhrcRwJMx6ldE4tXOnQB8Xwcm7p1LlxGw7HD4PHVXwMESV6xucfi2mbU9SJMBwLPh+9i+Z2aNzayFu2+MRf8iEp/+1RqGXwif7FvA1R/Dw/+AQ/cEnSpxWazuJR4u+LPcvX0px2YBD7j7u+H9BcAd7r6olLYjKPotgIyMjFBWVlal8uTl5ZGWllap11Yl5YqMckUmrnLl5PCH759nWcFabmw1kh7biq1qGYqPL1PF1XgVE02uzMzMHHfvUupBd4/JBrQFlpVxbDZwZrH9BUCovHOGQiGvrOzs7Eq/tiopV2SUKzLxkGve6nm+fPNyd/Ct9fCdqXj2Qw+5w4EtTsTDeJUmmlzAIi+jplbXF6/WA62L7bcCNlRT3yJSDbbt2cY1b1zDuVPP5U/v/AmAJt9DWn7AweRH1VXwZwBDw1frdAN2uObvRWqN11a+Rrux7Xju4+e4+8y7mXjhREhPL71xWc9LlYvJh7Zm9iLQE2hqZuuB+4AUAHcfB8wB+gOrgd3ANbHoV0SCN3XJVIZMH0Kn5p2Y+z9zObn5yUUHcnMPNFq4sGgyRwIVq6t0riznuAM3xaIvEQmeu/NN3je0SG/BJSdcwo7vdzAiNELr38Q5LZ4mIhFZu30t5049lzOePuPHxc5u6nqTin0NoKUVRKRCCr2QJ95/grsW3IWZ8Zfef6F+Sv2gY0kEVPBFpFxb92zlghcv4N/r/k3fY/oy7rxxtGncJuhYEiEVfBEpV+N6jclomMGzFz3L4JMGa/2bGkpz+CJSqg83fsjZz5zNxp0bSbIkXrviNYZ0HKJiX4Op4IvIT+zZt4e73ryLrhO6snLLStZsXxN0JIkRTemIyI/e/epdhs0Yxmfffca1J1/LQ30eokn9JkHHkhhRwReRH41bNI78gnzmD5lP76N7Bx1HYkwFXyTBzV01l9aHtKb94e15vN/jpCSnkJYafytISvQ0hy+SoL7b/R1Dpw+l/wv9eeDdBwBoUr+Jin0tpnf4IgnG3Xl1xauMmjuKrXu28psev+Ge7vcEHUuqgQq+SIKZumQqQ18fSqhFiPlD5nNSxklBR5JqooIvkgDcnQ07N3BEoyO47MTL2LVvF8M7D6dOkkpAItEcvkgtt2bbGvpM7cOZk89kV/4u6tWpx8guI1XsE5D+j4vUUgWFBYx5fwx3//Nuki2Zv57zVy12luBi8g7fzPqa2admttrM7izleE8z22Fmi8PbvbHoV0RK993u7zhz8pncOu9WerbtyfIblzOyy0iSTL/UJ7Ko3+GbWTLwBHAORfeu/cDMZrj7ihJN33H386PtT0TK16R+E1qmt2TqL6YyqMMgrX8jQGymdLoCq939CwAzywIGACULvohUoQ++/oBbFt/C3NBcWqa3ZNrl04KOJHHGPMr7TJrZpUBfdx8e3h8CnOruo4q16QlMo+g3gA3A7e6+vIzzjQBGAGRkZISysrIqlSsvL4+0tPj7AolyRUa5yvd9wfdMWTuFV9a/QuOUxvy+/e9p16hd0LF+Ip7Gq7jamCszMzPH3buUetDdo9qAy4CJxfaHAI+XaNMISAs/7g+sqsi5Q6GQV1Z2dnalX1uVlCsyynVw2Wuy/ZjRxzi/xa+bcZ3PnD8z6EilipfxKqk25gIWeRk1NRaf4KwHWhfbb0XRu/ji/6jkunte+PEcIMXMmsagb5GENumjSRR6IQuGLmD8BeNJqxN/71YlfsRiDv8D4FgzOwr4GhgIDCrewMyaA5vc3c2sK0VXB30Xg75FEs7sz2bTpnGbA4udJaXQMLVh0LGkBoj6Hb677wdGAfOAlcDL7r7czEaa2chws0uBZWb2MTAaGBj+1UNEKmjL7i0Mfm0w5794Pn/911+BolsPqthLRcXki1fhaZo5JZ4bV+zxGGBMLPoSSTTuzkvLX+LmuTez4/sd3HfWfdzd/e6gY0kNpG/aisS555Y8x1WvX8UpLU9h0oWT6JDRIehIUkOp4IvEoUIv5Ovcr2l9SGsuP/Fyvt//PcM6DSM5KTnoaFKD6XvWInFm9dbV9Hq2F90nd/9xsbMRoREq9hI1vcMXiRMFhQU8+t6j/Cb7N6Qkp/Bwn4dpkNIg6FhSi6jgi8SB73Z/R7/n+/HBhg+44LgLePK8Jzmi0RFBx5JaRgVfJA40qd+Eto3bcttpt3HFiVdosTOpEprDFwnI+1+/T/fJ3fk692uSLImXL3uZge0HqthLlVHBF6lmu/ft5vZ/3M5pk05jzbY1rMtdF3QkSRCa0hGpRtlrshk+czhfbPuCkaGRPND7AQ6pd0jQsSRBqOCLVKPJiyeTZEksvGohZ7U9K+g4kmBU8EWq2IxPZ3BU46PokNGhaLGz5BRdbimB0By+SBXZvGszA18dyICsATz47wcBOKTeISr2Ehi9wxeJMXfnhaUvcMvfb2Fn/k5+n/l7fn3Gr4OOJaKCLxJrz378LFe/cTXdWnVj0oWTaNcsvm43KIlLBV8kBgq9kPW56znykCO5ov0V7C/cz9UnX631bySuaA5fJEqrvlvF2c+c/ZPFzoZ11sqWEn9iUvDNrK+ZfWpmq83szlKOm5mNDh9fYmadY9GvSLVq1AjMiracHPYnGw+eYZz0t+NY/M1i7jvrPn0gK3Et6ikdM0sGngDOoeiG5h+Y2Qx3X1GsWT/g2PB2KvBk+E+RmmPnzh8f7vBdnDYMFh0BF62EJyasoGV6ywDDiZQvFu/wuwKr3f0Ld88HsoABJdoMAJ71Iu8Bjc2sRQz6FglEOvX5r23w8svw2kuo2EuNYNHeS9zMLgX6uvvw8P4Q4FR3H1WszSzgAXd/N7y/ALjD3ReVcr4RwAiAjIyMUFZWVqVy5eXlkZaWVqnXViXlikw85Vr+/jSe3DuT++oNpf6RJ5K2fv2Bg6FQcMGKiafxKk65IhNNrszMzBx371LqQXePagMuAyYW2x8CPF6izWzgzGL7C4BQeecOhUJeWdnZ2ZV+bVVSrsjEQ668vXl+y9xb3O7Dj7wV/88RePZDD7nDgS1OxMN4lUa5IhNNLmCRl1FTY3FZ5nqgdbH9VsCGSrQRiTtvfvEm1828jrXb13LTB/DnNyE9HxYGHUykEmIxh/8BcKyZHWVmqcBAYEaJNjOAoeGrdboBO9x9Ywz6FqlSU5dMJTU5lbevfpsx76STnl9Ko/T0as8lUhlRv8N39/1mNgqYByQDT7v7cjMbGT4+DpgD9AdWA7uBa6LtV6SqvP7J6xzd5GhOyjiJ0f1Gk5KUQv2U+pCbe6DRwoVFkzkiNUhMvmnr7nMoKurFnxtX7LEDN8WiL5GqsilvEzfPvZlXVrzCVR2vYspFU2hUt1HQsURiRksrSMJzd6Yumcqt824lLz+PP579R351+q+CjiUScyr4kvB+WOzs9NanM+nCSRzf9PigI4lUCRV8SUiFXsi6Heto07gNA9sPpNALGdpxqNa/kVpNi6dJwvl0y6ecNeUsuk/uTl5+HnXr1OWaTteo2Eutp4IvCWNfwT4eePcBOo7ryPLNy/l95u9pmNIw6Fgi1UZTOpIQtuzeQp/n+vDRNx9xyQmXMKb/GJqnNQ86lki1UsGXWs3dMTMOq38YJzQ7gXu638Ml7S4JOpZIIDSlI7XWv776F90mdWN97nrMjOcvfl7FXhKaCr7UOnn5efzv3P+l++TubMrbxMadWsVDBDSlI7XMPz7/ByNmjuCrHV8xquso/tTrT6Slxt/ytyJBUMGXWuWFpS9QP6U+71zzDmcceUbQcUTiigq+1HjTVkzjmEOPoWPzjozuN5rU5FTq1akXdCyRuKM5fKmxNu7cyCUvX8Klr1zKo/95FIBGdRup2IuUQe/wpcZxd6YsnsJt/7iNPfv28ECvB/jl6b8MOpZI3FPBlxpnyuIpXDvjWrof2Z2JF07kuMOOCzqSSI2ggi81QkFhAety19G2cVsGdRhEclIyg08aTJJpVlKkoqL622Jmh5rZfDNbFf6zSRnt1prZUjNbbGaLoulTEs/Kb1fSY0oPekzuwa78XdStU5ehHYeq2ItEKNq/MXcCC9z9WGBBeL8sme5+srt3ibJPSRD7CvYx9cupnPzUyXyy5RP+ePYfaZDSIOhYIjVWtFM6A4Ce4cfPAAuBO6I8pwjf7vqWc547h483fczlJ17O6L6jyUjLCDqWSI1mHsWNmM1su7s3Lra/zd1/Nq1jZmuAbYADT7n7+IOccwQwAiAjIyOUlZVVqWx5eXmkpcXfNyyV6+B+WOzM3fnzJ3/mlPRTOKfVOUHH+pl4Ga+SlCsytTFXZmZmTpkzKe5+0A14E1hWyjYA2F6i7bYyztEy/OfhwMdAj/L6dXdCoZBXVnZ2dqVfW5WUq2xvrX3Lu4zv4ut2rPvxuXjIVRrlioxyRSaaXMAiL6Omljul4+69yzpmZpvMrIW7bzSzFsDmMs6xIfznZjObDnQF3i6vb0kMuXtzuevNuxi7aCxHNT6Kb/K+oVWjVkHHEql1ov3QdgZwVfjxVcAbJRuYWUMzS//hMdCHot8QRJi7ai7tx7bnyUVPcuupt7L0hqV0aanP9UWqQrQf2j4AvGxmw4CvgMsAzKwlMNHd+wMZwHQz+6G/F9z971H2K7XEKyteIb1uOv8e9m+6teoWdByRWi2qgu/u3wG9Snl+A9A//PgLoGM0/Ujt4e68suIVjjvsOE5ufjKP9X2M1ORU6tapG3Q0kVpP31yRarNh5wYufvlirnj1Ckb/ZzQA6XXTVexFqomWVpAq5+48/dHT/PIfv2RvwV4ePOdBbu12a9CxRBKOCr5UuSmLpzB85nDOanMWEy+cyDGHHhN0JJGEpIIvVaKgsIAvd3zJ0U2OZlCHQaQkpzCowyCtfyMSIP3tk5hbvnk5Zzx9BmdNOevHxc60sqVI8PQ3UGImvyCf+9+6n05PdWL11tX8pfdftNiZSBzRlI7ExLe7vqXXs71YunkpV7a/ksf6Pkazhs2CjiUixajgS1Q8vNhZ0wZN6dSiE388+49c8N8XBB1LREqhKR2ptIVrF9JlQhfW7ViHmfHMRc+o2IvEMRV8idiO73cwctZIMp/JZPv329m8q9Q180QkzmhKRyIy+7PZXD/rejbmbeSXp/2S+zPv1wezIjWECr5EZNrKaTSp34TXrniNrkd0DTqOiERABV8Oyt15aflL/Pdh/02nFp14rO9j1K1Tl9Tk1KCjiUiENIcvZVqfu54BWQO4ctqVjHl/DFC02JmKvUjNpHf48jOFXsjEDyfyq/m/Yl/BPh7u8zC3nHpL0LFEJEpRvcM3s8vMbLmZFZpZmbcpMrO+Zvapma02szuj6VNiqFEjMCvacnLADDdj8hkNuH7W9XRp2YWlNyzlttNuIzkpOei0IhKlaN/hLwMuBp4qq4GZJQNPAOcA64EPzGyGu6+Ism+J1s6dPz4s8AJWHQrHbIUh7++l4fgXueLEKwjfqUxEaoGo3uG7+0p3/7ScZl2B1e7+hbvnA1nAgGj6ldhaejiM2jOGzKthdwqkFsDA9gNV7EVqGXP36E9ithC43d0XlXLsUqCvuw8P7w8BTnX3UWWcawQwAiAjIyOUlZVVqUx5eXmkpaVV6rVVKZ5y5S/6D8/nL+D5fQtIS07jf+tcSGadjkWFPhQKOh4QX+NVnHJFRrkiE02uzMzMHHcvfYrd3Q+6AW9SNHVTchtQrM1CoEsZr7+Mohua/7A/BHi8vH7dnVAo5JWVnZ1d6ddWpXjJtSlvk594I85v8f+5GH/9wd+5w4EtTsTLeJWkXJFRrshEkwtY5GXU1HLn8N29d6X+mTlgPdC62H4rYEOU55RK8vBiZ80aNOOUr+Ev8+G8VbDw9IZBRxORKlYd1+F/ABxrZkeZWSowEJhRDf1KCf9c8086j+/842Jnk/+ZznmrSmmYnl7t2USk6kV7WeYvzGw9cBow28zmhZ9vaWZzANx9PzAKmAesBF529+XRxZZIbP9+O9fNuI5ez/YiLz+PLbu3FB3IzT0wiRMKHXicmxtsYBGpElFdlunu04HppTy/AehfbH8OMCeavqRyZnw6gxtm38A3ed/w69N/zW97/pb6KfWDjiUiAdA3bWu5GZ/OoGmDprwx8A26tCzzu3EikgBU8GsZd+eFpS9wQrMT6NyiM4/2fZTU5FStfyMiWjytNlm3Yx3nv3g+g6cPZuwHYwFIS01TsRcRQO/wa4VCL+SpRU9xx5t3UOAFPHruo4zqWur32kQkgang1wJTFk/hxjk30vvo3ow/fzxHNTkq6EgiEodU8Guo/YX7WbNtDccediyDTxpMWmoal7W7TOvfiEiZNIdfA338zcd0m9iNzGcy2ZW/i9TkVC4/8XIVexE5KBX8GmTv/r385p+/ocuELqzLXcejfR/VDcRFpMI0pVNDbN61mZ5TerJyy0qGdhzKI30e4bAGhwUdS0RqEBX8OFd8sbPTW5/OI+c+Qt9j+gYdS0RqIE3pxLH5n8+n47iOfLXjK8yMiRdOVLEXkUpTwY9D2/ZsY9gbw+gztQ97C/aydc/WoCOJSC2gKZ04M33ldG6ccyPf7vqWu868i3vPupd6deoFHUtEagEV/Dgze9Vsmqc1Z/ag2XRu0TnoOCJSi6jgB8zdeW7Jc7Q/vD2dW3Tmsb6PkZqcSkpyStDRRKSW0Rx+gL7c/iX9nu/HVa9fxbhF4wBomNpQxV5EqkS0d7y6zMyWm1mhmZW52LqZrTWzpWa22MwWRdNnbVDohTzx/hO0f7I97371LqP7jmbc+eOCjiUitVy0UzrLgIuBpyrQNtPdt0TZX60w+aPJjJo7inOOPofxF4ynbeO2QUcSkQQQ7S0OVwJaw6UC9hXsY832NQAM6TiERnUbcWm7SzV2IlJtzN2jP4nZQuB2dy91usbM1gDbAAeecvfxBznXCGAEQEZGRigrK6tSmfLy8khLS6vUa2Nt1c5VPPjZg2zL38aT7Z6k6SFNg470M/E0XsUpV2SUKzK1MVdmZmaOu5c+xe7uB92ANymauim5DSjWZiHQ5SDnaBn+83DgY6BHef26O6FQyCsrOzu70q+NlT379vhdb97lyb9L9owHM3zaimlxkas0yhUZ5YqMckUmmlzAIi+jppY7pePuvSv1z8xPz7Eh/OdmM5sOdAXejva88Wzzrs30mNyDT7/7lGtOvoaH+zxMk/pNWLhpYdDRRCRBVfllmWbW0MzSf3gM9KHoN4RaycNTZM0aNKNHmx7MGzyPpwc8TZP6TQJOJiKJLtrLMn9hZuuB04DZZjYv/HxLM5sTbpYBvGtmHwPvA7Pd/e/R9Buv5q2ex0njTuLL7V9iZoy/YDx9/qtP0LFERIDor9KZDkwv5fkNQP/w4y+AjtH0E++27tnKbfNu45mPn+H4psez/fvttKFN0LFERH5CSytEadqKadw05ya27N7CPd3v4f/2+L9a7ExE4pIKfpTmfT6Plukt+fvgv3Ny85ODjiMiUiYV/Ai5O1MWT6FDRge6tOzC3879G3Xr1KVOkoZSROKbFk+LwJpta+gztQ/XzriWCTkTgKLFzlTsRaQmUKWqgILCAp744AnuWnAXSZbE2P5jub7L9UHHEhGJiAp+BUxZPIVb/n4L/Y7px7jzx3HkIUcGHUlEJGIq+GXYV7CPz7d9zvFNj2dox6EcWv9QLjr+Ii12JiI1lubwS/Hhxg85ZcIpnP3M2ezK30VKcgq/OOEXKvYiUqOp4BezZ98e7nzzTrpO6MrmXZsZe95YGqY2DDqWiEhMaEonbFPeJrpP7s6qrasY1mkYD/V5iMb1GgcdS0QkZhK+4Bd6IUmWxOENDyezbSZPnvckvY7uFXQsEZGYS+gpnTmr5tB+bHvWbl+LmfHUBU+p2ItIrZWQBX/L7i0MmT6E8144DzMjd29u0JFERKpcwk3pvLz8ZUbNGcW277dxb497ubv73dStUzfoWCIiVS7hCv6CLxbQpnEbFly4gA4ZHYKOIyJSbWp9wXd3nv7oaU7KOIlTjjiFv/X9G6nJqVr/RkQSTrR3vHrQzD4xsyVmNt3MGpfRrq+ZfWpmq83szmj6PKhGjcCsaMvJ4YtDjd5XJzF85nAmfTQJgAYpDVTsRSQhRfuh7XygvbufBHwG3FWygZklA08A/YB2wJVm1i7Kfku3cycABQav5L9NhxvggyNg3EwYe97YKulSRKSmiKrgu/s/3H1/ePc9oFUpzboCq939C3fPB7KAAdH0W57JnWBs/gwy18KKJ+D6HEiyhLwgSUTkR+busTmR2UzgJXefWuL5S4G+7j48vD8EONXdR5VxnhHACICMjIxQVlZWxUPk5ACw3wtY2Phbem3POLD+TSgU4X9R1cjLyyMtLS3oGD+jXJFRrsgoV2SiyZWZmZnj7l1KPejuB92AN4FlpWwDirW5h6KbmVspr78MmFhsfwjweHn9ujuhUMgjAj9u2Q899JP9eJGdnR10hFIpV2SUKzLKFZlocgGLvIyaWu6nl+7e+2DHzewq4HygV7izktYDrYvttwI2lNeviIjEVrRX6fQF7gAudPfdZTT7ADjWzI4ys1RgIDAjmn7LlJ4e2fMiIgkk2k8yxwDpwHwzW2xm4wDMrKWZzQHwog91RwHzgJXAy+6+PMp+S5ebe2ASJxQ68DhXSyeIiER1Qbq7H1PG8xuA/sX25wBzoulLRESio2sVRUQShAq+iEiCUMEXEUkQKvgiIgkiZt+0rQpm9i3wZSVf3hTYEsM4saJckVGuyChXZGpjrjbu3qy0A3Fd8KNhZou8rK8XB0i5IqNckVGuyCRaLk3piIgkCBV8EZEEUZsL/vigA5RBuSKjXJFRrsgkVK5aO4cvIiI/VZvf4YuISDEq+CIiCaJGF/zybo5uRUaHjy8xs85xkqunme0IrzC62MzuraZcT5vZZjNbVsbxoMarvFxBjVdrM8s2s5VmttzMbimlTbWPWQVzVfuYmVk9M3vfzD4O5/pdKW2CGK+K5ArkZyzcd7KZfWRms0o5FtvxKuvOKPG+AcnA58DRQCrwMdCuRJv+wFzAgG7Af+IkV09gVgBj1gPoDCwr43i1j1cFcwU1Xi2AzuHH6cBncfIzVpFc1T5m4TFICz9OAf4DdIuD8apIrkB+xsJ93wa8UFr/sR6vmvwOvyI3Rx8APOtF3gMam1mLOMgVCHd/G9h6kCZBjFdFcgXC3Te6+4fhxzspup/DESWaVfuYVTBXtQuPQV54NyW8lbwqJIjxqkiuQJhZK+A8YGIZTWI6XjW54B8BrCu2v56f/9BXpE0QuQBOC/+KOdfMTqziTBUVxHhVVKDjZWZtgU4UvTssLtAxO0guCGDMwtMTi4HNwHx3j4vxqkAuCOZn7FHg10BhGcdjOl41ueBbKc+V/Fe7Im1irSJ9fkjRehcdgceB16s4U0UFMV4VEeh4mVkaMA241d1L3j4tsDErJ1cgY+buBe5+MkX3ru5qZu1LNAlkvCqQq9rHy8zOBza7e87BmpXyXKXHqyYX/IrcHD2IG6iX26e75/7wK6YX3Q0sxcyaVnGuiojLG84HOV5mlkJRUX3e3V8rpUkgY1ZerqB/xtx9O7AQ6FviUKA/Y2XlCmi8zgAuNLO1FE39nm1mU0u0iel41eSCX5Gbo88AhoY/6e4G7HD3jUHnMrPmZmbhx10p+v/wXRXnqoggxqtcQY1XuM9JwEp3f6SMZtU+ZhXJFcSYmVkzM2scflwf6A18UqJZEONVbq4gxsvd73L3Vu7elqI68U93H1yiWUzHK6p72gbJ3feb2Q83R08Gnnb35WY2Mnx8HEX30e0PrAZ2A9fESa5LgRvMbD+wBxjo4Y/kq5KZvUjR1QhNzWw9cB9FH2AFNl4VzBXIeFH0DmwIsDQ8/wtwN3BksWxBjFlFcgUxZi2AZ8wsmaKC+bK7zwr672QFcwX1M/YzVTleWlpBRCRB1OQpHRERiYAKvohIglDBFxFJECr4IiIJQgVfRCRBqOCLiCQIFXwRkQTx/wFp4jfTVrSOOwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def grad_pas_fixe_reg_pol(x , y , d , rho = 1 , Max_iter = 1000 , tol = 1e-8):\n", " \"\"\"Réalise la méthode du gradient à pas fixe pour faire une régression polynomiale\n", " Entrées:\n", " - x: Array de taille (n,), abscisses\n", " - y: Array de taille (d,), ordonnées\n", " - d: Int - Degré de la régression polynomiale\n", " - rho: Flottant - Pas de la méthode du gradient à pas fixe. Defaut: 1\n", " - Max_iter: Int - Nomnbre maximal d'itérations et critère d'arrêt. Defaut: 1000\n", " - tol: Tolérance, variation du gradient maximale et critère d'arrêt. Defaut: 1e-8\"\"\"\n", " \n", " n = np.size(x)\n", " xx = x.copy()\n", " xx = xx.reshape(n,1)\n", " \n", " X = np.ones_like(xx)\n", " \n", " for i in range(1,d+1):\n", " X = np.concatenate((X,xx**i),axis=1)\n", " \n", " alpha_k = np.random.uniform(size=(d+1,))\n", " k = 0\n", " \n", " ALPHA = [alpha_k] # Historique des itérations alpha_0,alpha_1,...\n", " \n", " while (k < Max_iter) and (np.linalg.norm(2*(X.T@(X@alpha_k-y))) > tol):\n", " k+=1\n", " #print(k,np.linalg.norm(2*(X.T@(X@alpha_k-y))))\n", " alpha_k = alpha_k - rho*2*(X.T@X@alpha_k-X.T@y)\n", " ALPHA.append(alpha_k)\n", " \n", " if (np.linalg.norm(2*(X.T@X@alpha_k-X.T@y)) > tol):\n", " raise RuntimeError(f\"La méthode n'a pas convergé après {Max_iter} itérations.\")\n", " else:\n", " return ALPHA,alpha_k\n", " \n", " \n", "ALPHA , alpha = grad_pas_fixe_reg_pol(x , y , d = 1 , rho = 1e-2 , Max_iter = 1000 , tol = 1e-8)\n", "print(\"Solution: y =\" , alpha[1] , \"x +\" , alpha[0])\n", "\n", "plt.figure()\n", "plt.scatter(x,y,marker=\"s\",color=\"red\")\n", "plt.plot(x,y,color=\"green\",linestyle=\"dashed\")\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2 - La loi d'Ohm\n", "\n", "La loi d'Ohm permet de faire le lien entre la tension aux bornes d'une resistance et l'intensité la traversant. Si on note $R$ la valeur de la résistance, $U$ la tension aux bornes de la résistance et $I$ l'intensité la traversant, alors on a:\n", "\n", "$$ U = R \\times I$$.\n", "\n", "On mesure les tensions et intensités suivantes (repectivement en Volts et Ampères):" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "U = np.array([0.085 , 0.187 , 0.4 , 0.52 , 0.71 , 0.785 , 1.02])\n", "I = np.array([0.012 , 0.021 , 0.045 , 0.049 , 0.067 , 0.08 , 0.11])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q4 -** Déterminer la valeur de la resistance $R$ par régérssion linéaire. Tracer $U$ en fonction de $I$ sur un graphe, ainsi que la droite obtenue par régression linéaire. On prendra $\\rho = 10^{-1}$, $\\verb|Max_iter| = 10^5$ ,$\\verb|tol| = 10^{-8}$. On comparera les résultats avec la méthode $\\verb|np.polyfit|$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "- Méthode programmée:\n", " > Solution: U = 9.75385145422322 I + -0.0054969940211336416\n", " > R = 9.75385145422322 Omhs\n", " \n", "- Méthode np.polyfit:\n", " > Solution: U = 9.753852175346076 I + -0.005497033618984676\n", " > R = 9.753852175346076 Omhs\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyYUlEQVR4nO3deZxN9f/A8dfbvg1lDaOYUg3Zt5JilOyUbwplLUNfheRbqq+tSNtXSSFfKaUs3ygSJT9GC2WJlFH4WifZxjKLMMv798e95juNO2PMzJ1zl/fz8biP5p7zOee8P+fmvu/n8znnc0RVMcYYE7wKOB2AMcYYZ1kiMMaYIGeJwBhjgpwlAmOMCXKWCIwxJshZIjDGmCBnicBkm4hcLSIJIlLQ6ViyIiL3iMhBd6wN8vG4K0Skbz4cp5+IfJvDbaNE5GEPy33ys3XHFJaNcsVFZJ2IdMiPuAKNJQI/5/6HcuGVKiJ/pnv/QF4eS1UPqGopVU3Jy/16wavAo+5Yt3jjACIyTkTmpl+mqu1VdY43jucN7oTyHlzeZysi1UVERaSQt2N0x7THfdz3RGRCJkXfBl5V1eXejikQef2DNN6lqqUu/C0i+4CHVXWVcxH5hGuA7U4HYfKPqvZxOgZ/Zi2CACUiBURklIj8V0RiRWShiJR1r7vwi66viBwQkeMi8my6bZuKyCYRiRORIyIyOcN2hdzvq4jIUhE5ISK7RWRgun2Mcx/zfRGJF5HtItI4i3hvFJGv3Pv6TUTuS7fuPRF5S0Q+d+/rBxG51sM+iopIAlAQ+ElE/uteHu7uEjnljqNLdvctIrXTxXVERJ4RkXbAM8D97pbXT+6yad0u7vP/TxHZLyJH3eehTHbOv4d6lXOf5zgR2QBcm2F9pucuuzx8tlEi8ryIfOc+LytFpLy7+Nfu/55y1/8W9zYDRGSHiJwUkS9F5Jp0+1cRGSwiu9zr3xIRca+7TkTWishp97lYkGG760QkEngAeNJ9zM/c66uIyCIROSYie0Vk6OXW3QCqaq8AeQH7gDvdfw8HvgdCgaK4ms7z3OuqAwr8GygO1APOAeHu9euB3u6/SwE3Z9iukPv9WmAaUAyoDxwD7nCvGwecBTrg+mKeBHyfSdwlgYNAf1yt1IbAcaC2e/17wAmgqXv9h8D8LM6DAte5/y4M7Mb1xV0EaA3EAzdcat9ACPAH8IS7jiFAs3T1m5vhuFG4WmQAA9zHDXOfw8XAB9k5/x7qMx9Y6D5PNwG/A99m59x52FdajBmWZ/xso4D/Ate7Y4wCXvRU1r3sbnd9w91x/BNYl+EzWQZcAVzt/n+lnXvdPOBZXD9MiwEtMvks3wMmpFtXANgMjHF/tmHAHqCt0/8W/e1lLYLANQh4VlVjVPUcri+ue+Wv/brjVfVPVf0J+AnXFxJAEnCdiJRX1QRV/T7jzkWkGtACeEpVz6rqVmAW0DtdsW9Vdbm6+p0/SLf/jDoB+1T1XVVNVtUfgUXAvenKLFbVDaqajOvLun42z8PNuL6IX1TV86q6GtcXUs9s7LsTcFhV/+WuY7yq/pDN4z4ATFbVPaqaADwN9Mjm+U8jrsHbvwFjVDVRVX8B0o9DZOfc5dS7qrpTVf/ElYjqZ1F2EDBJVXe4z+MLQP30rQJcn8EpVT0ArEm3vyRc3XlV3Oc5uwPhTYAKqvqc+7Pdgyu59shuBY2LJYLAdQ3wibs75BSwA0gBKqUrczjd32dwfWECPITrl+CvIrJRRDp52H8V4ISqxqdbth+omsX+i4nnAcZrgGYXYnXH+wBwVTZivZQqwEFVTb2MOC/suxquX8U5UcV9nPTHLET2zn96FdzbHcywrwuyc+5y6nLO+TXAlHQxnACE7J3nJ91lN7i77gZkM75rgCoZ6v4Mfz3HJhtssDhwHQQGqOp3GVeISPWsNlTVXUBPESkAdAM+FpFyGYodAsqKSEi6ZHA1rm6LnMS6VlXb5GDbSzkEVBORAumSwdXAzmzG1TOTdZeatvcQri+qC64GkoEjuLrrsuuYe7tqwK/p9pU+Rm+du8x4qvtBYKKqfnjZO1M9DAwEEJEWwCoR+VpVd1/iuAeBvapa83KPaf7KWgSBawYw8ULTXEQqiEjX7GwoIg+KSAX3F+cp9+K/XFaoqgeBdcAkESkmInVxtSQu+4sAV1fN9SLSW0QKu19NRCQ8B/vK6AcgEdcgY2ERaQV0xtXvnp24rhKR4eIaiA4RkWbudUeA6u5k6ck84HERqSEipXB1lSxwd5tkm7tbbTEwTkRKiEgtIP29Ct48d5k5BqTi6pO/YAbwtIjUBhCRMiLSPTs7E5HuInIhOZ7E9YXv6TLWIxmOuQGIE5GnxHUfQUERuUlEmlxmfYKeJYLANQVYCqwUkXhcA8fNst4kTTtgu7iuwJkC9FDVsx7K9cQ1cHgI+AQYq6pfXW6g7hbFXbj6dg/h6kJ4Cdcgd66o6nmgC9Ae1yDqNKCPqv6a5Yb/i6sNrsRxGNgFRLhX/8f931gR+dHD5rNxjYt8DezFNXD+WA6r8SiubpTDuAZM380Qo1fOXWZU9QwwEfjO3SVzs6p+4j7ufBGJA37Bdc6zownwg/v/t6XAMFXd66HcO0At9zE/dSfJzrjGGvbi+nxnAWVyUb2gJKr2YBpjjAlm1iIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyPndfQTly5fX6tWrOx1GphITEylZsqTTYeQ7q3fwCda6+2u9N2/efFxVK3ha53eJoHr16mzatMnpMDIVFRVFq1atnA4j31m9g0+w1t1f6y0i+zNbZ11DxhgT5CwRGGNMkLNEYIwxQc7vxgg8SUpKIiYmhrNnPc2CkL/KlCnDjh07nA4jTxUrVozQ0FAKFy7sdCjGGC8IiEQQExNDSEgI1atXx/3QI8fEx8cTEhLiaAx5SVWJjY0lJiaGGjVqOB2OMcYLAqJr6OzZs5QrV87xJBCIRIRy5cr5RGvLGOMdAZEIAEsCXmTn1pjAFjCJwBhjAlVSUhK//fab1/ZvicCPJScn8+abb3Lu3DmnQzHGeMlPP/1Eo0aNaN26NWfOnPHKMSwR5JFSpbL7CF0YM2YMq1atytXxVJXhw4dTt25dihb1/AyS6tWrc/z4cQCaN28OwL59+/joo49ydWxjTP4pUqQIiYmJTJs2jRIlSnjlGAFx1ZC/ee6553K8bXJyMoUKFUJEePPNN7O93bp164D/JYJevXrlOAZjjHd9+umnrFmzhilTphAeHs7OnTspWLCg144XkInA0zwg9913H3//+985c+YMHTp0uGh9v3796NevH8ePH+fee+/9y7qoqKgcxbF161YGDx7MmTNnuPbaa5k9ezZXXnkl/fr1o1OnThcdp1WrVtSvX58NGzYQFxfH7Nmzadq0KePGjePQoUPs27eP8uXLM2XKFAYPHsyBAwcAeP3117n11luJjY2lZ8+eHDt2jKZNm5L+6XOlSpUiISGBUaNGsWPHDurXr0/fvn0ZOnQoo0aNIioqinPnzjFkyBAGDRqUo/oaY3Ln999/57HHHuOTTz6hTp06xMXFUbp0aa8mAbCuIa/q06cPL730Etu2baNOnTqMHz/+ktskJiaybt06pk2bxoABA9KWb968mSVLlvDRRx8xbNgwHn/8cTZu3MiiRYt4+OGHARg/fjwtWrRgy5YtdOnSJS1RpPfiiy9y2223sXXrVh5//HHeeecdypQpw8aNG9m4cSP//ve/2bvX0+NijTHekpqayrRp0wgPD2fFihW8+OKLbN68mdKlS+fL8QOyRZDVL/gSJUpkub58+fI5bgGkd/r0aU6dOkXLli0B6Nu3L927d7/kdj179gTg9ttvJy4ujlOnTgHQpUsXihcvDsCqVauIjo5O2yYuLo74+Hi+/vprFi9eDEDHjh258sorL3m8lStXsm3bNj7++OO0uHft2mU3jxmTj06cOMHo0aNp1qwZM2bM4Nprr83X4wdkIvBnGa/Zv/A+/fznqamprF+/Pi0xZLX9pagqU6dOpW3btjmI1hiTU2fPnmX27NkMHjyY8uXLs3HjRmrUqOHIfTvWNeQlZcqU4corr+Sbb74B4IMPPkhrHWRlwYIFAHz77beUKVOGMmXKXFTmrrvu+stA8datWwFXK+LDDz8EYMWKFZw8efKibUNCQoiPj09737ZtW6ZPn05SUhIAO3fuJDExMZu1NMbkxOrVq6lTpw5Dhgxh9erVAISFhXlMAqUnlUbGy0Wv0pPyrtvIWgR55MyZM4SGhqKqiAgjRoxgzpw5aYPFYWFhvPvuu5fcz5VXXknz5s3TBos9eeONNxgyZAh169YlOTmZ22+/nRkzZjB27Fh69uxJw4YNadmyJVdfffVF29atW5dChQpRr149+vXrx7Bhw9i3bx8NGzZEValQoQKffvppbk+HMcaD2NhYRo4cyXvvvcd1113HqlWruOOOO7LcJv58/GUtzxFV9atXo0aNNKPo6OiLljklLi4ux9u2bNlSN27cmIfR5J1LneM1a9bkTyA+JljrrRq8dc9NvW+55RYtVKiQPvPMM3rmzJlsbcM4Mn1dDmCTZvK9ai0CY4zxoj179lCpUiVKlizJ5MmTKVmyJHXq1HE6rL/w2hiBiMwWkaMi8ksm60VE3hCR3SKyTUQaeisWfxEVFUXjxo2dDsMYkweSkpJ4+eWXuemmm5g4cSIAN998s88lAfDuYPF7QLss1rcHarpfkcB0L8ZijDH5ZuPGjTRp0oSnnnqKdu3aMWTIEKdDypLXEoGqfg2cyKJIV+B9d/fV98AVIlLZW/EYY0x+mD59Os2aNePYsWMsXryYxYsXU7Vq1RzvL6SI5wddZbY8J0TTTUOQ10SkOrBMVW/ysG4Z8KKqfut+/3/AU6q6yUPZSFytBipVqtRo/vz5f1lfpkwZrrvuuryvQA6kpKR4/XZwJ+zevZvTp09nuj4hIeGyJt4LFMFabwjeumdW76SkJAoXLsz+/ftZsmQJDz300F/u/3FaRETEZlX12Pfs5GCxp7smPGYlVZ0JzARo3LixZpxLaMeOHT7zeMhAe1TlBcWKFaNBgwaZro+KivI4x1OgC9Z6Q/DWPWO9Dx06xNChQylYsGDafUB9+/Z1KLqccfKGshigWrr3ocAhh2LJNRGhd+/eae+Tk5OpUKECnTp1cjAqY4y3pKamMmPGDMLDw1m2bBn169fHmz0s3uRki2Ap8KiIzAeaAadV9Q9vH7T0pNIeb8QIKRJC3NNxOd5vyZIl+eWXX/jzzz8JCQnhq6++ylW/YHZcmJLaGJO/9uzZQ+/evVm3bh2tW7dmxowZ1KxZ0+mwcsybl4/OA9YDN4hIjIg8JCKDRWSwu8hyYA+wG/g38HdvxZKeN+/Sa9++PV9++SUA8+bNS5tADlyzig4YMIAmTZrQoEEDlixZAsD27dtp2rQp9evXp27duuzatYt9+/Zx003/G1Z59dVXGTduHOCaqvqZZ56hZcuWTJkyhc2bN9OyZUsaNWpE27Zt+eMPVy594403qFWrFnXr1qVHjx65rpsx5n9KlSrF8ePHmTNnDqtWrfLrJABebBGoas9LrFfAt6+pukw9evRgzJgxdO/enW3btjFgwIC0uYYmTpxI69atmT17NqdOnaJp06bceeedzJgxg2HDhvHAAw9w/vx5UlJSOHLkSJbHOXXqFGvXriUpKYmWLVuyZMkSKlSowIIFC3j22WeZPXs2L774Inv37qVo0aJpM5gaY3IuKiqKWbNm0b9/fypWrEh0dHTAXBhi/Qp5qG7duhw4cIB58+Zd9PCblStXsnTpUl599VXANfPggQMHuOWWW5g4cSIxMTF069YtW78s7r//fgB+++03fvnlF9q0aQO4rliqXLlyWiwPPPAAd999N3fffXce1tKY4BIbG8s//vEP3n33XcLCwujcuTNAwCQBsESQ59q3b8/IkSOJiooiNjY2bbmqsmjRIm644Ya/lA8PD6dZs2Z8/vnntG3bllmzZnH99deTmpqaVubs2bN/2ebCJWmqSu3atVm/fv1FcXz++ed8/fXXLF26lOeff57t27fbeIIxl0FVmTdvHsOHD+fEiROMGjWK0aNHs2HDBqdDy3M2DXUe6927N2PGjLnoNvK2bdsyderUtKsKtmzZArgGncLCwhg6dChdunRh27ZtVKpUiaNHjxIbG8u5c+dYtmyZx2PdcMMNHDt2LC0RJCUlsX37dlJTUzl48CARERG8/PLLnDp1ioSEBC/W2pjAc/78ecaOHUuNGjXYvHkzkyZN8trD450WdD8RQ4qEZHrVUF6oWrUqw4YNu2j56NGjGT58OHXr1kVVqV69OsuWLWPBggXMnTuXwoULc9VVVzFmzBgKFy7MmDFjaNasGTVq1ODGG2/0eKwiRYrw8ccfM3ToUE6fPk1ycjLDhw/n+uuv58EHH+T06dOoKo8//jhXXHFFntTPmECWlJTEzJkz6devHyVLluT//u//qFq1akB1A3mU2bSkvvoK5GmofZlNQ+1ZsNZbNfDqvmHDBq1Xr54COnv27EzL+Wu9yWIaausaMsYEtfj4eIYPH87NN9+cNj9Q//79nQ4rXwVd15AxxqQ3cOBAFi5cyCOPPMILL7zg8fGwgS5gEoG6HxFp8p766W3zxmTmjz/+oECBAlSqVInx48czdOhQmjdv7nRYjgmIrqFixYoRGxtrX1heoKrExsZSrFgxp0MxJtdSU1N5++23CQ8PZ8SIEYDr6rtgTgIQIC2C0NBQYmJiOHbsmNOhcPbs2YD70ixWrBihoaFOh2FMrkRHRxMZGcl3331HRERE2rQtJkASQeHChalRo4bTYQCu29Czmq7ZGJP/PvnkE+6//35CQkJ499136du3r3UlpxMQXUPGGOPJuXPnAGjRogUDBgzg119/pV+/fpYEMrBEYIwJOCdOnODhhx+mVatWpKSkUKFCBWbMmEGFChWcDs0nWSIwxgQMVWX+/PmEh4fz3nvvcfvtt5OcnOx0WD4vIMYIjDHm6NGj9OvXjxUrVtCkSRNWrlxJvXr1nA7LL1iLwBgTEEJCQjh8+DBTpkxh/fr1lgQugyUCY4zf2rx5M/fccw+JiYkUL16cTZs2pT1I3mSfJQJjjN9JSEhgxIgRNG3alB9++IFdu3YBUKCAfaXlhJ01Y4xfWb58ObVr1+a1115j0KBB7Nixg/r16zsdll+zwWJjjN9QVV566SVKlSrFt99+y6233up0SAHBEoExxqelpqYye/ZsOnbsSOXKlVmwYAFly5alSJEiTocWMKxryBjjs3bs2EHLli0ZOHAgM2fOBOCqq66yJJDHLBEYY3zOuXPnGDduHPXq1WP79u288847jBkzxumwApZ1DRljfM7o0aN55ZVX6NWrF6+99hoVK1Z0OqSAZonAGOMTTp48ycmTJwkLC2PkyJG0bt2adu3aOR1WULCuIWOMo1SVBQsWEB4eTp8+fVBVKlasaEkgH1kiMMY4Zv/+/XTq1IkePXoQGhrK1KlTbYpoB1jXkDHGEevWraNNmzaICK+//jqPPvqoTQ3hEGsRGGPy1Z9//glAw4YN6dOnD9HR0QwbNsySgIMsERhj8kVCQgJPPPEEderUITExkWLFijF9+nSuvvpqp0MLepYIjDFet3z5cm666SYmT57MHXfcQUpKitMhmXS8mghEpJ2I/CYiu0VklIf1ZUTkMxH5SUS2i0h/b8ZjjMlfiYmJ9OzZk44dO1KiRAm++eYb3n77bUqXLu10aCYdryUCESkIvAW0B2oBPUWkVoZiQ4BoVa0HtAL+JSJ277gxAaJ48eLExsby3HPPsWXLFlq0aOF0SMYDb7YImgK7VXWPqp4H5gNdM5RRIERc14uVAk4A9oBRY/zYb7/9RpcuXTh06BAFChTgiy++YPTo0RQtWtTp0EwmvHn5aFXgYLr3MUCzDGXeBJYCh4AQ4H5VTc24IxGJBCIBKlWqRFRUlDfizRMJCQk+HZ+3WL2Dw5bDW0h1/xMNLRrKv+b9C4ACUoDaZWszb948PvzwQ4oWLcq8efNo1KiRk+F6RSB+5t5MBJ7uCtEM79sCW4HWwLXAVyLyjarG/WUj1ZnATIDGjRtrq1at8jzYvBIVFYUvx+ctVu/gEDE+Iu3vV69/lZE7R7re7Ifw78PZsWMHPXv25LXXXqNSpUoOReldgfiZezMRxADV0r0PxfXLP73+wIuqqsBuEdkL3Ahs8GJcxpi8thnOnDnD8uXLad++vdPRmMvkzUSwEagpIjWA34EeQK8MZQ4AdwDfiEgl4AZgjxdjMsbkBYWfvv8JkoBKQHvYPno7JUuWdDoykwNeSwSqmiwijwJfAgWB2aq6XUQGu9fPAJ4H3hORn3F1JT2lqse9FZMxJg+cApbDBzs/gAa4LgEpjiUBP+bVuYZUdTmwPMOyGen+PgTc5c0YjDF5IyUlBdYDq13vOz/Qmc/CPnM0JpM37M5iY0y2TJ8+3dW+rw78HVp2bOlq6wMhRUIcjMzkls0+aozJVGJiIvv27aN27do8/PDDhIaG0rVrV0SEqKgotGfGCwGNP7IWgTHGoy+++IKbbrqJTp06kZSURLFixbj77rvteQEByBKBMeYvjhw5Qq9evWjfvj1FixZlzpw5FC5c2OmwjBdZ15AxJs2uXbto1qwZiYmJjB07lqefftqmhggClgiMMZw5c4YSJUpw3XXX0a9fPwYOHEh4eLjTYZl8Yl1DxgSx8+fPM2HCBGrUqMGhQ4cQESZPnmxJIMhYi8CYIPXdd98RGRlJdHQ09913H4UK2ddBsLIWgTFBJiUlhUceeYQWLVqQkJDAsmXLWLBgARUrVnQ6NOMQ+wlgjI8oPak08efjL1oeUiSEuKfjPGyRMwULFuTs2bM8/vjjPPfcc5QqVSrP9m38k7UIjPERnpJAVssvx8GDB+nWrRs//fQTALNnz2by5MmWBAxgicCYgJaSksKUKVOoVasWX375JTt27ACwm8LMX1giMCZAbd26lVtuuYXhw4fTokULtm/fTo8ePZwOy/ggGyMwJkAtXryY/fv3M2/ePO6//35rBZhMWYvAmADy1VdfsXq1a57oZ555hh07dtCjRw9LAiZLlgiM8RGZTeWcnSmejx07Ru/evbnrrrt46aWXAChWrBhly5bN0xhNYLKuIWN8RE4uEVVV5syZwxNPPEF8fDxjxozh6aef9kJ0JpBZIjDGj3322Wf079+fW2+9lZkzZ1KrVi2nQzJ+yLqGjPEz58+fZ/PmzQB07tyZxYsX8/XXX1sSMDlmicAYP7J+/XoaNmxI69atOXHiBCLCPffcQ4EC9k/Z5Jz932OMHzh9+jRDhgzh1ltvJS4ujg8//NAGgk2esTECY3zciRMnqFOnDocPH2b48OE2P5DJc5YIjPFRiYmJlCxZkrJlyzJ48GDat29P48aNnQ7LBCDrGjLGx6SkpDB16lSqVauWNknc6NGjLQkYr7FEYIwP2bZtG82bN2fo0KE0adKEMmXKOB2SCQKWCIzxEWPHjqVRo0bs3buXuXPn8sUXX1C9enWnwzJBwBKBMT6kd+/e7NixgwceeMDmBzL5xgaLjXHIsWPHeOKJJ+jRowcdOnRg3Lhx9uVvHJGtRCAiBYB6QBXgT2C7qh7xZmDGBCpV5f333+eJJ54gLi4ubRDYkoBxSpaJQESuBZ4C7gR2AceAYsD1InIGeBuYo6qp3g7UmECwe/duBg0axOrVq2nevDkzZ86kdu3aTodlgtylWgQTgOnAIFXV9CtEpCLQC+gNzPFOeMYElrVr17Jp0yamT59OZGSkTQ1hfMKlEkEfVU3ytEJVjwKvZ7WxiLQDpgAFgVmq+qKHMq3c+ykMHFfVlpcK2hh/8v3333Pw4EG6d+/OgAED6NSpE5UqVXI6LGPSXOrnyO8i8m8RaS2X2YEpIgWBt4D2QC2gp4jUylDmCmAa0EVVawPdL+cYxviyuLg4Hn30UZo3b8748eNJSUlBRCwJGJ9zqUQQDmwCRgMHReR1EWmWzX03BXar6h5VPQ/MB7pmKNMLWKyqByCtlWGM3/v000+pVasW06ZN47HHHmP9+vUULFjQ6bCM8UgydP1nXlCkCq5f7D2AisB8VX02i/L3Au1U9WH3+95AM1V9NF2Z13F1CdUGQoApqvq+h31FApEAlSpVajR//vxsxeyEhISEoJwQzOr9P7t372bgwIGEhYUxcuRIwsPDHYrOu+wz9y8RERGbVdXzPCWqmu0XUAroA2wFjlyibHdc4wIX3vcGpmYo8ybwPVASKI/ryqTrs9pvo0aN1JetWbPG6RAcEez1Tk5O1nXr1qUt/+yzz/T8+fMORZU/gv0z9zfAJs3ke/WSlyyISDER6S4ii4H/AncAT+O6pyArMUC1dO9DgUMeynyhqomqehz4Gtf9Csb4jZ9//pkWLVpw2223sXPnTgA6depE4cKFHY7MmOzJMhGIyEfAAeB+4CPgGlXtq6orVDXlEvveCNQUkRoiUgRXl9LSDGWWALeJSCERKQE0A3bkpCLG5Lc///yTWbNm0bBhQ3bv3s2cOXOoWbOm02EZc9kudfnol7juIYi/3B2rarKIPOreR0FgtqpuF5HB7vUzVHWHiHwBbANScXUl/XK5xzImvyUlJdG4cWOio6Pp378/r7zyCuXKlXM6LGNy5FKJIAVIyGyl+87jyqr6raf1qrocWJ5h2YwM718BXslWtMY4LD4+npCQEAoXLsyQIUM4e/YsI0aMcDosY3LlUmME5YCtIjJbRIaIyH0i0kdEnhORtcDLgM05ZAKeqjJ37lzCwsL4/PPPAfj73/9Ow4YNHY7MmNzLMhGo6hSgITAPqIBroLgh8DvQW1X/pqq7vB6lMQ7673//S9u2benduzc1a9a0ZwSYgHPJ2Ufdg8JfuV/GBJW33nqLkSNHUqRIEaZNm8agQYNsfiATcOx5BMZkoUSJEnTo0IE33niDqlWrOh2OMV5hP22MSScuLo6hQ4cyY4brmoZ+/fqxaNEiSwImoFkiMMZtyZIl1KpVizfffJMDBw4A9rAYExwu9WCajNfFKXAc+FZV93otKmPy0e+//87QoUNZvHgxderUYdGiRTRrlt25FY3xf5caIwjxsKw68KyIjFNV3539zQSc0pNKE3/+4nsbQ4qEEPd0XI73++uvv7J8+XImTZrEE088YVNDmKCTZSJQ1fGelotIWWAVrqmljckXnpJAVsuz8ssvv7Bu3ToiIyO544472L9/PxUrVsxtiMb4pRyNEajqCcA6T43fOXv2LP/85z9p0KABY8eOJSHBdeO8JQETzHKUCESkNXAyj2MxxqvWrFlD3bp1mThxIr169eLnn3/2y3nljclrlxos/hnXAHF6ZXFNJ93HW0EZk9cOHz5Mu3btqFatGl999RV33nmn0yEZ4zMuNVjcKcN7BWJVNdFL8RiTZ1SVr7/+mpYtW3LVVVfx+eefc+utt1K8eHGnQzPGp1xqrqH9GV4HLAkYp4QU8XQRm+fle/bsoW3btrRq1YqoqCgA7rzzTksCxnhgU0wYv5GdS0STkpKYPHky48ePp1ChQrz55pvcdttt+RCdMf7LEoEJGKpKu3btWL16NXfffTdTp04lNDTU6bCM8Xk2xYTxewkJCaSmpiIiDBo0iMWLF/PJJ59YEjAmmywRGL/22WefUatWrbRJ4u677z7uueceh6Myxr9YIjB+6dChQ9x777106dKFMmXK2JPCjMkFSwTG7yxYsIDw8HCWLVvGCy+8wI8//sjNN9/sdFjG+C0bLDZ+p0KFCjRp0oTp06dTs2ZNp8Mxxu9ZIjA+7+zZs0ycOBFVZcKECbRu3ZqIiAh7VoAxecS6hoxPi4qKom7dukyYMIFDhw6h6prxxJKAMXnHEoHxSSdOnOChhx4iIiKClJQUvvrqK2bPnm0JwBgvsERgfNKRI0eYP38+o0aN4ueff7ZJ4ozxIhsjMD5j7969/Oc//+HJJ58kPDycAwcOUK5cOafDMibgWYvAOC45OZlXX32V2rVr8/zzz3Pw4EEASwLG5BNLBMZRmzZtokmTJvzjH/+gTZs2REdHU61aNafDMiaoWNeQccyff/5Jhw4dKFSoEB9//DHdunWzwWBjHGCJwOS7qKgobr/9dooXL86nn35K7dq1KVOmjNNhGRO0rGvI5Js//viD++67j4iICObOnQtA8+bNLQkY4zCvJgIRaSciv4nIbhEZlUW5JiKSIiL3ejMe44zU1FRmzpxJeHg4S5cuZcKECfTo0cPpsIwxbl7rGhKRgsBbQBsgBtgoIktVNdpDuZeAL70Vi3FW//79ef/994mIiGDGjBlcf/31AJSeVJr48/EXlQ8pEpKtp5EZY/KGN8cImgK7VXUPgIjMB7oC0RnKPQYsApp4MRaTz86ePYuqUrx4cfr3709ERAR9+/b9y2CwpySQ1XJjjHd4s2uoKnAw3fsY97I0IlIVuAeY4cU4TD5bu3Yt9evXZ+zYsQC0atWKfv362RVBxvgob7YIPP2r1wzvXweeUtWUrL4kRCQSiASoVKkSUVFReRRi3ktISPDp+LwlISGBpUuX8vbbb7N8+XIqV65M+fLlszwXr17/aqbr/OUcBuvnDcFb90CstzcTQQyQ/s6gUOBQhjKNgfnuJFAe6CAiyar6afpCqjoTmAnQuHFjbdWqlZdCzr2oqCh8OT5v+de//sXLL79MbGwsTz75JGPHjqVEiRJZbhMxPiLTddoz428G3xSsnzcEb90Dsd7eTAQbgZoiUgP4HegB9EpfQFVrXPhbRN4DlmVMAsY/lC9fnpo1a/Lll19Sv359p8MxxlwGryUCVU0WkUdxXQ1UEJitqttFZLB7vY0L+LHk5GRef/11tm/fzrvvvss111zDt99+e1n7CCkSkulVQ8aY/OPVO4tVdTmwPMMyjwlAVft5MxaTdzZt2kRkZCRbtmyhc+fOnDt3Lkf7sUtEjfENdmexybaEhARGjBhBs2bNOHz4MB9//DFLliyhaNGiTodmjMkFSwQm286cOcMHH3xAZGQk0dHR/O1vf7NLQo0JAJYITJYOHz7MmDFjSElJoWLFivz2229Mnz6dK664wunQjDF5xBKB8Sg1NZVZs2YRHh7OSy+9xJYtWwAoW7asw5EZY/KaJQJzkV9//ZVWrVoxcOBA6tWrx7Zt22jcuLHTYRljvMSeR2D+QlXp3r07v//+O++88w79+/e3cQBjApwlAgPAd999R/369SlZsiRz586lcuXKVKxY0emwjDH5wLqGgtzJkyeJjIykRYsWTJ48GYB69epZEjAmiFiLIEipKgsXLmTYsGEcP36ckSNHMmLECKfDMsY4wBJBkBo9ejQTJ06kUaNGrFixggYNGjgdkjHGIZYIgkhycjKJiYmUKVOGBx98kPLly/PYY49RsGBBp0MzxjjIEkGQ2LJlCwMHDiQsLIyFCxdy4403cuONNzodljHGB9hgcYBLTExk5MiRNGnShJiYGO69916nQzLG+BhrEQSwH3/8kW7durF//34iIyN58cUXufLKK50OyxjjYywRBLBq1apRrVo1PvjgA2677TanwzHG+CjrGgogF+YH6tChAykpKVSoUIFvvvnGkoAxJkuWCALEr7/+SkREBAMHDiQxMZGTJ086HZIxxk9YIvBz58+f57nnnkubHG7WrFmsWbOG8uXLOx2aMcZP2BiBn0tNTWXu3Ll069aN1157jauuusrpkIwxfsZaBH7o1KlTjBo1isTERIoVK8aGDRuYN2+eJQFjTI5YIvAjqsp//vMfwsPDeeWVV1i9ejWAPS3MGJMrlgj8xIEDB+jSpQv33XcfVapUYePGjXTu3NnpsIwxAcASgZ945JFHWL16NZMnT+aHH36gYcOGTodkjAkQNljsw7Zu3UqlSpWoXLkyU6dOpUCBAlSvXt3psIwxAcZaBD4oMTGRJ598ksaNGzN69GgAwsLCLAkYY7zCWgQ+5ssvv2Tw4MHs27ePgQMH8tJLLzkdkjEmwFmLwIfMmDGDdu3aUbRoUdauXcvMmTNtkjhjjNdZi8BhqsqJEycoV64cf/vb34iNjWXkyJEULVrU6dCMMUHCWgQO2rlzJ61bt6Zjx45pk8Q9++yzlgSMMfnKEoEDzp8/z4QJE6hbty5btmzhoYceQkScDssYE6Ssayif7dmzh86dOxMdHc3999/P66+/blNDGGMcZS2CfKKqAFSpUoWqVauybNky5s+fb0nAGOM4ryYCEWknIr+JyG4RGeVh/QMiss39Wici9bwZjxNUlUWLFtGiRYu0SeJWrlxJx44dnQ7NGGMALyYCESkIvAW0B2oBPUWkVoZie4GWqloXeB6Y6a14nHDw4EG6du3Kvffey5kzZzhy5IjTIRljzEW8OUbQFNitqnsARGQ+0BWIvlBAVdelK/89EOrFePLFlsNbiBgbARuA1YACd8Hu23YTFhbmcHTGGHMxbyaCqsDBdO9jgGZZlH8IWOFphYhEApEAlSpVIioqKo9CzHtVilTh5etfZsbCGRQOL0y3/t0oV7EcgE/HnVsJCQkBXb/MBGu9IXjrHoj19mYi8HQ9pHosKBKBKxG08LReVWfi7jZq3LixtmrVKo9CzDtnzpzhhRdeoECVAjx/7Hm4GygKk05NglOuMtrTY/UDQlRUFL74uXhbsNYbgrfugVhvbw4WxwDV0r0PBQ5lLCQidYFZQFdVjfViPF6zcuVKbrrpJiZOnEj0j+6er2J4ToXGGONjvJkINgI1RaSGiBQBegBL0xcQkauBxUBvVd3pxVi84ujRozz44IO0bduWwoULs2bNGm654xanwzLGmMvitUSgqsnAo8CXwA5goapuF5HBIjLYXWwMUA6YJiJbRWSTt+LxhrFjx7Jw4ULGjBnDTz/9FHDNRWNMcPDqncWquhxYnmHZjHR/Pww87M0Y8tquXbtISUnhxhtv5LnnnuOxxx6jVq3/XRVbQDzn1pAiIfkVojHGXBabYiKbzp8/zyuvvMLzzz/P7bffzsqVK6lQoQIVKlT4S7kGVzVAewTuoLAxJvBYIsiGdevWERkZyfbt2+nevTtTpkxxOiRjjMkzlggu4bPPPqNr166EhoaydOlSOnfu7HRIxhiTp2zSuUwcO3YMgDZt2jB+/Hi2b99uScAYE5AsEWQQExPD3XffTZMmTdImiRs9ejQhITbYa4wJTJYI3FJSUnjzzTepVasWK1euZMiQIRQpUsTpsIwxxutsjACIjY2lY8eO/PDDD7Rt25bp06dTo0YNp8Myxph8EdQtggsPiylbtixXX301H374IStWrLAkYIwJKkGbCFatWkXjxo05dOgQIsLChQvp1auXPTvYGBN0gi4RHDt2jD59+tCmTRvi4+M5evSo0yEZY4yjgioRvP/++4SHhzN//nz++c9/sm3bNurXr+90WMYY46igGixetWoVN9xwAzNnzqR27dpOh2OMMT4hqBLB9OnTKV68OAUKBFVDyBhjshRUiaBkyZJOh2CMMT4n4BNB6UmliT8ff9HykCIhxD0d50BExhjjWwK+j8RTEshquTHGBJuATwTGGGOyZonAGGOCnCUCY4wJcpYIjDEmyAV8IsjsofH2MHljjHEJ+MtH7RJRY4zJWsC3CIwxxmTNEoExxgQ5SwTGGBPkLBEYY0yQs0RgjDFBTi48t9dfiMgxYL/TcWShPHDc6SAcYPUOPsFad3+t9zWqWsHTCr9LBL5ORDapamOn48hvVu/gE6x1D8R6W9eQMcYEOUsExhgT5CwR5L2ZTgfgEKt38AnWugdcvW2MwBhjgpy1CIwxJshZIjDGmCBnieAyiEg7EflNRHaLyCgP60VE3nCv3yYiDd3Lq4nIGhHZISLbRWRY/kefczmtd7r1BUVki4gsy7+ocy839RaRK0TkYxH51f2535K/0edcLuv9uPv/8V9EZJ6IFMvf6HMuG/W+UUTWi8g5ERl5Odv6PFW1VzZeQEHgv0AYUAT4CaiVoUwHYAUgwM3AD+7llYGG7r9DgJ0Zt/XVV27qnW79COAjYJnT9cmvegNzgIfdfxcBrnC6Tt6uN1AV2AsUd79fCPRzuk55WO+KQBNgIjDycrb19Ze1CLKvKbBbVfeo6nlgPtA1Q5muwPvq8j1whYhUVtU/VPVHAFWNB3bg+kfjD3JcbwARCQU6ArPyM+g8kON6i0hp4HbgHQBVPa+qp/Ix9tzI1eeN6xknxUWkEFACOJRfgefSJeutqkdVdSOQdLnb+jpLBNlXFTiY7n0MF3+ZX7KMiFQHGgA/5H2IXpHber8OPAmkeik+b8lNvcOAY8C77i6xWSJS0pvB5qEc11tVfwdeBQ4AfwCnVXWlF2PNS9mptze29QmWCLJPPCzLeO1tlmVEpBSwCBiuqv7y6LQc11tEOgFHVXVz3ofldbn5vAsBDYHpqtoASAT8pd84N5/3lbh+CdcAqgAlReTBPI7PW7JTb29s6xMsEWRfDFAt3ftQLm72ZlpGRArjSgIfqupiL8aZ13JT71uBLiKyD1dzubWIzPVeqHkqN/WOAWJU9UKr72NcicEf5KbedwJ7VfWYqiYBi4HmXow1L2Wn3t7Y1idYIsi+jUBNEakhIkWAHsDSDGWWAn3cV1XcjKtp/IeICK7+4h2qOjl/w861HNdbVZ9W1VBVre7ebrWq+ssvxNzU+zBwUERucJe7A4jOt8hzJ8f1xtUldLOIlHD/P38HrvEwf5CdentjW9/g9Gi1P71wXS2xE9cVAs+6lw0GBrv/FuAt9/qfgcbu5S1wNRW3AVvdrw5O18fb9c6wj1b40VVDua03UB/Y5P7MPwWudLo++VTv8cCvwC/AB0BRp+uTh/W+Ctev/zjglPvv0plt608vm2LCGGOCnHUNGWNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhjTJCzRGBMNohIQibLi4vIWhEpmG7Z4yJyVkTKpFtWR0Tey4dQjblslgiMyZ0BwGJVTUm3rCeum4zuubBAVX8GQkXk6nyOz5hLskRgTO48ACy58EZErgVKAf/ElRDS+wzXXafG+BRLBMbkkHs6gTBV3ZducU9gHvANcIOIVEy3bhNwW/5FaEz2WCIwJufK45pqIL0ewHxVTcU16Vr3dOuO4pqV0xifUsjpAIzxY38CaY9iFJG6QE3gK9ecaxQB9uCalwd32T/zOUZjLslaBMbkkKqeBAqmey5vT2CcqlZ3v6oAVUXkGvf663FNxmaMT7FEYEzurMQ1uyy4uoU+ybD+E/43QBwBfJ5PcRmTbTb7qDG5ICINgBGq2vsS5YoCa4EWqpqcL8EZk03WIjAmF1R1C7Am/Q1lmbgaGGVJwPgiaxEYY0yQsxaBMcYEOUsExhgT5CwRGGNMkLNEYIwxQc4SgTHGBLn/B1hDPDBEjt9AAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Méthode programmée\n", "\n", "ALPHA , alpha = grad_pas_fixe_reg_pol(I , U , d = 1 , rho = 0.1 , Max_iter = 10**5 , tol = 1e-8)\n", "print(\"- Méthode programmée:\")\n", "print(\" > Solution: U =\" , alpha[1] , \"I +\" , alpha[0])\n", "print(\" > R = \",alpha[1],\" Omhs\")\n", "print(\" \")\n", "\n", "\n", "# Méthode np.polyfit\n", "\n", "alpha_bis =np.polyfit(I , U , 1)\n", "print(\"- Méthode np.polyfit:\")\n", "print(\" > Solution: U =\" , alpha_bis[0] , \"I +\" , alpha_bis[1])\n", "print(\" > R = \",alpha_bis[0],\" Omhs\")\n", "\n", "plt.figure()\n", "plt.scatter(I , U , color = \"green\" , marker = \"s\" , label = \"Mesures\")\n", "plt.plot(I , alpha[0]+alpha[1]*I , color = \"black\" , linestyle = \"dashed\" , label=\"Loi prédite\")\n", "plt.xlabel(\"I (A)\")\n", "plt.ylabel(\"U (V)\")\n", "plt.legend()\n", "plt.title(\"Tension en fonction de l'intensité\")\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3 - Trajectoire d'un ballon\n", "\n", "Un ballon, lancé depuis une hauteur $y_0$ avec une vitesse initiale $v_0$ et selon un angle $\\theta$, et soumis uniquement à la force de pesanteur (on ne prend pas en compte les effets dûs aux frottements de l'air), décrit une trajectoire parabolique dont l'expression de l'altitude $y(t)$ au cours du temps est décrite par l'expression suivante:\n", "\n", "$$y(t) = -\\frac12 gt^2 + v_0\\sin(\\theta)t + y_0,$$\n", "\n", "où $g = 9.81 m.s^{-2}$ est l'accélération de la pensanteur. On fait les relevés d'altitude suivants au cours du temps, toutes les 100 millisecondes:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "t = np.array([0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0])\n", "#y = -0.5*9.81*t**2 + 5*np.sin(np.pi/4)*t + 1.8 + 0.05*np.random.uniform(size=(11,))\n", "#print(y)\n", "y = np.array([1.84 , 2.15 , 2.31 , 2.45 , 2.48 , 2.36 , 2.19 , 1.92 , 1.49 , 1.05 , 0.44])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Q5 -** En utilisant la fonction créée en question **Q3**, trouver l'équation de la trajectoire en fonction du temps (on se placera ainsi dans le cas $d=2$). On prendra $\\rho = 10^{-2}$, $\\verb|Max_iter| = 10^5$ ,$\\verb|tol| = 10^{-8}$. En supposant que le ballon a été lancé avec un angle $\\theta = 45°$, quelle est la vitesse initiale du ballon ? Représenter $t$ et $y$ sur un graphe, ainsi que la parabole obtenuer par régression polynomiale. Comparer les résultats obtenus en utilisant la fonction $\\verb|np.polyfit|$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "- Méthode programmée:\n", "Trajectoire: y(t) = -4.924242341444443 t^2 + 3.5415150651934546 t + 1.8327272869625104\n", "Vitesse initiale: v_0 = 5.00845863654522 m/s\n", "Accélération de la pesanteur calculée: g = 9.848484682888886 m/s**2\n", " \n", "- Méthode np.polyfit:\n", "Trajectoire: y(t) = -4.92424242424243 t^2 + 3.5415151515151573 t + 1.8327272727272725\n", "Vitesse initiale: v_0 = 5.008458758622543 m/s\n", "Accélération de la pesanteur calculée: g = 9.84848484848486 m/s**2\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9SklEQVR4nO3de5xN5f7A8c+XGQYjcmlUmFFRwpCREpWhkDqSckrDqZgzR25diBA6+emiy08uhaJ+1dQkujgq6TJDJcmEkCPCDKnjGsa4zcz398feM2eMPTN7LnvW7Nnf9+u1X7PXep611veZYX/3etZazyOqijHGGJNXJacDMMYYUz5ZgjDGGOORJQhjjDEeWYIwxhjjkSUIY4wxHlmCMMYY45ElCOPXRORTEbmnjI85TkRe9eH+d4rIDcXcVkXkkpIeV0QeF5G3irMfU3EEOR2ACTwikpZrsTpwEsh0L/9DVeO93Zeq3lQK8TwOXKKq/b085pMlPaYx/sAShClzqhqa/V5EdgKxqvpF3noiEqSqGWUZW0n5Y8zG5Me6mEy5ISKdRWS3iIwRkT+A10TkXBFZIiL7ROSQ+33DXNskiUhsruWBIrLZXfczEQnPVdZCRD4XkYMi8h93V1EPYBxwp4ikich6d90LRGSxu+42Efl7rv3kdL+ISIS7W2eQiKQCXxUWh4d2DxCRFBE5ICLj85S9LiL/k/d3VMivsqeIbBeR/SLyrIhUcm97sYh85T7OfhGJF5Hahewr+7i9RGSTiPzp/p03z1W2U0RGichPInJYRN4VkRBv9mvKN0sQprxpANQBwoE4XP9GX3MvNwaOAzM9bSgivXF92PcB6gNfA++4y2oCXwBLgQuAS4AvVXUp8CTwrqqGqmpr9+7eAXa7694BPCkiXQuI+3qgOdC9oDg8xHw58DIwwH2sukBDT3WL4DagHdAWuBUYmH044Cn3cZoDjYDHC9uZiDTDFf+DuNrzCfAvEamSq9pfgR5AEyASuLeEbTDlgCUIU95kAZNU9aSqHlfVA6q6SFXTVfUoMAXXh7En/wCeUtXN7m6eJ4E27m/vtwB/qOrzqnpCVY+q6veediIijYBOwBh33XXAq7g+xPPzuKoeU9XjhcSR1x3AElVdoaongQnu30FJPKOqB1U1FZgG9ANQ1W2q+rn7d7sPeIH8f5e53Ql87N72NPAcUA24Jled6aq6R1UPAv8C2pSwDaYcsARhypt9qnoie0FEqovIHHcXzBFgBVBbRCp72DYceNHdDfIncBDXt+YLcX1b/tXLGC4ADroTUrYU937ys8vLODwdK2dbVT0GHPAyTm9iSXEfAxE5T0QSROQ39+/yLaCeF/u7wL2f7Biz3MfI3Z4/cr1PB0Ixfs8ShClv8g4vPBK4FLhKVc8BrnOvFw/b7sJ1F1TtXK9qqrrSXXaxl8fcA9Rxd0tlawz85mXcBcWR1++4kperUSLVcXUzZTuG606vbA0KiCFbo1zvG+NqD7i6lxSIdP8u++P595jXHlxJLztGcR+joN+HqQAsQZjyriau6w5/ikgdYFIBdWcDY0WkBYCI1BKRvu6yJUADEXlQRKqKSE0Rucpd9h8gIvtirqruAlYCT4lIiIhEAoMAb2+/LSiOvBYCt4hIJ3ef/hOc+f9yHa6LznVEpAGu6wCFecR9cb8R8ADwrnt9TSAN1+/yQuARL9uzALhZRLqKSDCupH0S1+/IVGCWIEx5Nw1Xf/d+YBWui8weqeoHwDNAgrsLZSNwk7vsKHAj8Bdc3SFbgWj3pu+5fx4QkR/d7/sBEbi+PX+A67rI594EXFAcHupuAoYCb+M6mziE6+J4tjeB9cBOYBn//bAvyEdAMq7k8jEwz73+n7guXB92r3/fy/ZswXW2MQPX3+EvwF9U9ZQ32xv/JTZhkPFnIrICeFVV33A6FmMqGjuDMH7L3V9/EbDD6ViMqYgsQRi/JCLn4eoqWg5843A4xlRI1sVkjDHGIzuDMMYY41GFGqyvXr16GhERUaxtjx07Ro0aNUo3oHLO2lzxBVp7wdpcVMnJyftVtb6nsgqVICIiIlizZk2xtk1KSqJz586lG1A5Z22u+AKtvWBtLioRScmvzLqYjDHGeGQJwhhjjEeWIIwxxnjkswQhIo1EJNE9acomEXnAQ53O7glG1rlfE3OV9RCRLe7JWh71VZzGGGM88+VF6gxgpKr+6B4VM1lEPlfVn/PU+1pVb8m9wj2U8yxcY+fsBn4QkcUetjXGGOMjPjuDUNXfVfVH9/ujwGYKHk8/t/bANlXd7h4QLAHXzFjGFE98PEREQHKy62e8twOzGhO4yuRJahGJwDXRS0tVPZJrfWdgEa6zhD3AKFXdJCJ3AD1UNdZdbwCu+QCGedh3HK6pKQkLC4tKSEgoVoxpaWmEhgbWHCcB0+aDByElBbKySGvYkNDdu6FSJQgPhzp1nI7OpwLmb5yLtblooqOjk1W1ncdCVfXpC9fMUslAHw9l5wCh7vc9ga3u931xjdCZXW8AMKOwY0VFRWlxJSYmFntbfxUwbQ4PVwVV0MTnnst5r+HhTkfmcwHzN87F2lw0wBrN5zPVp3cxuScXWQTEq+pZY8+r6hFVTXO//wQIFpF6uM4ocs+K1ZD/zopl/FV2N0+lSmXbzZOaWrT1xhjAhxep3dMSzgM2q+oL+dRpAPxHVVVE2uO6JnIA+BNoKiJNcE1reBdwt69iNWUgPh7i4iA93bWckuJaBoiJOaNqeno66enpnDhxgpMnT3LihGuK6hYtWgCwatUq9uzZc0Z5rVq1uPtu1z+RGTNm8Ouvv+aUnahWjYvT0/kf9/5347oYJo0b+7jRxvg3X97F1BFX19AGEVnnXjcO1xy5qOps4A7gfhHJwDWt5F3uU54MERkGfAZUBuara+Yt46eyxo2jkjs5TMB1WnkiPZ0T99zDieHDqVmzJikprif+Y2Ji+PDDD8/YPiIigh07XNM+TJw4kc8/P3Nyt5YtW+YkiEWLFrF27VpCQkIICQmh6jnnUOnkScjMRFW5EhAReoSH02PBAm688UbOPfdcXzbfGL/kswShqt9QyIToqjoTmJlP2SfAJz4IzZSRzMxMli9fTkJCAv9KTWUTUAcIA1oCIUDVzExCYmKoVatWznaDBg2ia9eurg/3qlUJCQk5o3zmzJkcP378vwmgalWqVauWU56UlHR2MPHxMH48mVlZPFW3Lp9ecgkf/PQTr915J5UqVeLpp5/mkUce+W/fayV7htSYCjVYnykfUlJSeP7553nvvff4448/qFGjBr1r1ODosWPUAYa5X4DrTqIZM87Y/pZbbqEgzZo1K3pQMTEQE0NQUhL37t/PvUBGRgarV69m6dKlXH311QCsWbOGm2++me7du3PTTTfRrVs36tWrV/TjGVMB2NckU2Kqyo8//sj69esBOHXqFK+88godO3bkvffeY+/evbw1Zw7h1aufuWH16jBligMRuwQFBXHNNdfwxBNPcO211wJQtWpVbrzxRj799FNiYmI477zzuOqqq/j1118di9MYp9gZhCm2n3/+mYSEBBISEti6dSt9+/ZlwYIFNG3alP379585Pn32hejx4113DzVu7EoOeS5QOy0yMpL4+HgyMzNJTk5m6dKlfPnll1xwwQUAPP/886xZs4abbrqJ7t27ExYW5nDExviOJQhTLL179+ajjz6iUqVKdOnShdGjR9OnT5+cco+Tl7i7efxB5cqVad++Pe3bt2fixJwhwjh+/DhfffUV2Q9ktm3bljvuuIOxY8c6FaoxPmNdTKZQu3fv5oUXXuCGG27g1KlTAPTq1YsZM2bw22+/8fnnnxMbG0udCv5UMsBjjz3G77//TnJyMlOmTKF69er88MMPOeWjRo1i/vz5/Pbbbw5GaUzpsDMI49GBAwdISEjg3Xff5euvvwYgKiqK3377jSZNmjBw4ECHI3ROpUqVaNu2LW3btmXcuHFkZWUBruEOEhIScpJDZGQkNzVuTExyMq3++KPcdqsZkx87gzA5Dh06xO+//w7Av//9b4YNG8aBAweYPHkyv/zyC2vWrKFJkyYOR1n+ZN8SGxoayq5du1i/fj3PPPMMdU6f5vklS/j+999BFc1+ONAGCjR+whJEgDt69Cjx8fH85S9/ISwsjCeffBKADh06sGHDBjZt2sRjjz1G06ZNHY7UP4gIkZGRjB49msT0dA7iGgYA4EXglvR0fh092sEIjfGeJYgA9vzzzxMWFkb//v1Zt24dI0aMYNCgQYDrW3HLli0djtDPpaZSE9dolQBVgOVAiz17mDhxIunZw44YU05ZgggwBw4cyHlfs2ZNBg4cyDfffENKSgrPPfccbdq0cS64iibPWE9DgC3A7dWrM3nyZFq0aMGXX37pSGjGeMMSRIDIyspi9uzZNGnShH+NHAkREcRFRTFzyRI67txpQ0v4wpQprocBc7mgenXi584lKSmJmjVrEhwc7FBwxhTO7mIKANu3byc2NpbExERuaNmSVi+9BO4RUgsaVdWUUAEPB14PrFu3Licxjx49mqCgIMaPH+/5GRJjHGBfGyu4V155hVatWpGcnMwrr7zCsiNHiMhODtnS010fYqb0xcTAzp2QleX6mSsJZycHVeXAgQM89dRTNG/enIULF2ZPlGWMoyxBVHDBwcFcf/31bNy4kdjYWGTXLs8VbfIcx4gI8+bN45tvvqFOnTr07duXbt26sX37dqdDMwHOEkQFk5GRwdSpU3n11VcBuOeee/j4449p1Mg9QV9+k+TY5DmO69ixI2vWrGHmzJls3Lgx5wE8Y5ziswQhIo1EJFFENovIJhF5wEOdGBH5yf1aKSKtc5XtFJENIrJORNb4Ks6KZOPGjXTo0IExY8bkPP0sIrgm93PzcOHU6VFVzX8FBQUxdOhQdu7cySWXXALA4MGDSUhIsG4nU+Z8eQaRAYxU1ebA1cBQEbk8T50dwPWqGglMBubmKY9W1Taq2s6Hcfq906dPM3nyZNq2bUtKSgoLFizg9ddf91w5JgbmznXNwwCun3Pn2gXqcqZq1aoAHDlyhB9++IF+/frRtWtXNm2yiRVN2fFZglDV31X1R/f7o8BmXFMB566zUlUPuRdXAQ19FU9F9v333zNx4kTuuOMONm3aRN++fc88a8gr+8JpVNRZF05N+XLOOeewevVqXn75ZdatW0ebNm0YOXIkR48edTo0EwDK5BqEiEQAVwDfF1BtEPBprmUFlolIsojE+TA8v3Ty5EmWLVsGQKdOnVi3bh1vv/029evXdzgyU9oqV67M4MGD+eWXX7jvvvuIj48nIyPD6bBMABBf92uKSCiuEQamqOr7+dSJBl4COqnqAfe6C1R1j4icB3wODFfVFR62jQPiAMLCwqKyx+kvqrS0NEJDQwuvWA78/PPPTJ06lV27dhEfH0+DBg2KtR9/anNpqQhtzm5DZmYm06ZNo3fv3lx88cUF1g0k1uaiiY6OTs63Gz97knZfvIBg4DPg4QLqRAK/As0KqPM4MKqw40VFRWlxJSYmFnvbsnLs2DEdOXKkVqpUSRs1aqSffvppifbnD20ubRWpzZs2bdK6detq5cqVdcSIEXro0KGz6lSk9nrL2lw0wBrN5zPVl3cxCTAP2KyqL+RTpzHwPjBAVX/Jtb6GiNTMfg90Azb6KlZ/cPr0adq3b8/zzz9PXFwcGzdupEePHk6HZRx0+eWX88svvxAXF8eMGTO49NJL+b//+z/X7bHx8RARAcnJrp82xLgpBl8OtdERGABsEJF17nXjgMYAqjobmAjUBV5yX1TNUNepThjwgXtdEPC2qi71Yazl1smTJ6latSrBwcEMGzaMZs2a0aVLF6fDMuVEnTp1eOmll4iNjWXo0KFMnTqVu4FKQ4a4npAHG07FFJvPEoSqfgMUcCsNqGosEOth/Xag9dlbBJYvv/yS2NhYXnzxRXr16sXgwYOdDsmUU23btuXbb7/lP//5D8EdOnAkPZ0pQKfsJJE9nIolCFME9iR1OXT48GHi4uK44YYbqFKlit2ZZLxSqVIlzj//fEhN5SvgOeBvzzzDt9kVbDgVU0SWIJyS3UdcqdIZfcSfffYZLVq0YN68eYwePZp169bRoUMHR0M1fqZxY3oDPwKh1apxI/Cxe70xRWEJwgnx8a4+4ZQUUP1vH3F8PH/88Qe1a9dm1apVPPPMM1SrVs3paI2/cQ+n0hqYPnQolwO3Ah/07u1sXMbv2HwQThg//r8XEIEPgCPp6dwzfjx/27GDfv36UaVKFefiM/4t1zwU59asSWKjRoyIiKD9I484G5fxO3YG4QR3X/B+4E6gD677gTUlBRGx5GBKLtdwKjVTU3ltxQouvPBCMjMzeeONN2ykWOMVSxBOaNyYPcC1wIfAFOBLQLIH0DPGRxYtWsQ999zDfffdx+nTp50Ox5RzliAckDZhAteLsBvXGCLjgGAbctuUgb59+/LEE0/wxhtv0Lt3b44dO+Z0SKYcswThgNBBg7i/Xz+WNWjAdSI25LYpMyLChAkTmDNnDkuXLuWGG27g4MGDTodlyim7SF2Gtm7dyqFDh2jfvj0P29AHxkFxcXHUq1eP2NhYfv31V+rUqeN0SKYcsgRRRjZv3kzXrl2pUaMGmzdvJijIfvXGWX369OGGG27gnHPOAeDgwYOWKMwZrIupDGzYsIHOnTuTlZXFhx9+aMnBlBvZyeH111+nWbNmfP99QVO2mEBjCcLH1q5dS3R0NMHBwaxYsYIWLVo4HZIxZ7n22mupXbs2Xbp0YenSgBwX03hgCcLHpk+fTmhoKCtWrKBZs2ZOh2OMRxdffDHffvstzZo14y9/+Qvxdo3MYAnCZ9Q9U9+cOXNYuXIlF110kcMRGVOwsLAwkpKS6NSpEwMGDGDTpk1Oh2QcZgnCBxITE7nmmmvYv38/VapU4YILLnA6JGO8UqtWLT799FMWLVpk3aHGEkRp++yzz+jZsydHjx61ieWNXwoJCeG2224D4JtvvuH++++3f8sBypdTjjYSkUQR2Swim0TkAQ91RESmi8g2EflJRNrmKushIlvcZY/6Ks7StGTJEnr16sWll15KYmIiDRo0cDokY0pk5cqVzJ49mz59+pCea4BJExh8eQaRAYxU1ebA1cBQEbk8T52bgKbuVxzwMoCIVAZmucsvB/p52LZc+fTTT7ntttuIjIzkq6++skl+TIUwevRoZs2axZIlS+jWrRuHDh1yOiRThnyWIFT1d1X90f3+KLAZuDBPtVuBN9RlFVBbRM4H2gPbVHW7qp4CEtx1y63WrVtz11138cUXX9jDRqZCGTJkCAkJCaxevZrrrruOAwcOOB2SKSOSfbeNTw8iEgGsAFqq6pFc65cAT7vnr0ZEvgTGABFAD/ec1YjIAOAqVR3mYd9xuM4+CAsLi0pISChWjGlpaYSGhhZ5u7Vr1xIZGUnlypWLdVwnFbfN/izQ2lya7U1OTubLL79k5MiR5frfe6D9jaFkbY6Ojk5W1XYeC1XVpy8gFEgG+ngo+xjolGv5SyAK6Au8mmv9AGBGYceKiorS4kpMTCzyNnPmzFER0alTpxb7uE4qTpv9XaC12VftTU1N1TVr1vhk3yUVaH9j1ZK1GVij+Xym+vQuJhEJBhYB8ar6vocqu4FGuZYbAnsKWF9uzJw5k3/84x/cdNNNDB8+3OlwjClTgwcPpnPnznz++edOh2J8yJd3MQmuidI2q+oL+VRbDPzNfTfT1cBhVf0d+AFoKiJNRKQKcJe7brnw/PPPM3z4cHr37s37779PSEiI0yEZU6ZeffVVLrroIm6++WaK261ryj9fjhrXEVfX0AYRWedeNw5oDKCqs4FPgJ7ANiAduM9dliEiw4DPgMrAfFUtF491/vbbb0yaNIm//vWvvPXWWwQHBzsdkjFl7vzzz2f58uXceuut3H333ezbt8/OpCsgnyUIdV14lkLqKDA0n7JPcCWQcuXCCy9k1apVXHbZZTYqqwlotWvXZunSpdx9993MmzePuLg4qlatCvHxMH68a+71xo1dMyXaZFh+yT7hvKCqPProozRs2JDhw4fTsmVLp0MyplyoVq0a7733HocPH6Zq1aqceO01goYOJej4cVeFlBSIi3O9tyThd2yojUKoKg8++CBTp05ly5YtOYPwGWNcgoKCqFu3LllZWdw5dCh9jx/nRO4K6emuMwrjdyxBFCArK4v777+f6dOn8+CDDzJjxgxc196NMXlVqlSJrseP8yHQAzicuzA11ZGYTMlYgsiHqhIbG8ucOXN49NFHeeGFFyw5GFOIEeHhvA2sBP4CnMouaNzYsZhM8VmCyIeI0Lp1ax5//HGefPJJSw7GeGPKFPpVr87rwNfAIwDVq7suVBu/Yxep8zh16hRbtmyhVatWPPDAWQPQGmMK4r4Qfff48exJSaHr+efDs8/aBWo/ZWcQ8fEQEQHJyZwMD+eOq6+mY8eO/Oc//3E6MmP8U0wM7NzJKFWu2LMHYmL4448/nI7KFENgJ4j4eNcteCkpnDx9mltTU/nX2rU83acPYWFhTkdnTIXwzDPP0KpVK3bu3Ol0KKaIAjtBjB8P6ekcA8bOm8cy4FVgSFKSs3EZU4HcdtttnD59mltvvZW0tDSnwzFFENgJwn3r3TRg/a+/8gYwKNd6Y0zJNWvWjHfffZeNGzdy7733kpWV5XRIxkuBnSDct96NAf73/vvpn2e9MaZ0dO/enWeffZZFixbxP//zP06HY7wU2HcxTZkCcXEEpacTedFFrnV2S54xPvHQQw+xbds2LrvsMqdDMV4K7ASRfetd9jAA4eE2sJgxPiIivPTSSznLp06dokqVKg5GZAoT2F1MkHNLHlFRrp+WHIzxufj4eFq1asX+/fudDsUUwBKEMabMNW3alJSUFPr27cvp06edDsfkw5czys0Xkb0isjGf8kdEZJ37tVFEMkWkjrtsp4hscJet8VWMxhhntG/fnldeeYWkpCQeeughp8Mx+fDlGcTruAZ19EhVn1XVNqraBhgLLFfVg7mqRLvL2/kwRmOMQwYMGMDIkSOZNWsWc+fOdToc44EvZ5RbISIRXlbvB7zjq1iMMeXTM888w6ZNm2wojnJKfDkBjjtBLFHVfKdgE5HqwG7gkuwzCBHZARwCFJijqvl+vRCROCAOICwsLKq4E6inpaURGhparG39lbW54vOH9mZmZlK5cuVS258/tLm0laTN0dHRyfn21Kiqz15ABLCxkDp3Av/Ks+4C98/zgPXAdd4cLyoqSosrMTGx2Nv6K2tzxedP7f3666+1e/fueuzYsRLtx5/aXFpK0mZgjebzmVoe7mK6izzdS6q6x/1zL/AB0N6BuIwxZejIkSMsW7aM++67z6b2LSccTRAiUgu4Hvgo17oaIlIz+z3QDfB4J5QxpuLo2bMnTz/9NAsWLODJJ590OhyDDy9Si8g7QGegnojsBiYBwQCqOttd7TZgmaoey7VpGPCBewa3IOBtVV3qqziNMeXHI488wk8//cRjjz1Gq1at6NWrl9MhBTRf3sXUz4s6r+O6HTb3uu1Aa99EZYwpz0SEV155hS1btrBw4UJLEA4L7LGYjDHlTrVq1Vi2bBm1atVyOpSAVx4uUhtjzBnOPfdcKlWqRGpqKiNGjCAjI8PpkAKSnUEYY8qtFStWMGPGDESEF1980elwAo4lCGNMudW/f3+Sk5OZNm0arVu3ZuDAgU6HFFCsi8kYU649++yz3HjjjQwePJiVK1c6HU5AsQRhjCnXgoKCSEhIoHHjxkyaNMnpcAKKdTEZY8q9OnXqsGzZMurXr+90KAHFziCMMX7hoosuombNmhw7doyZM2facBxlwBKEMcavvP322wwfPpypU6c6HUqFZwnCGONXYmNjueuuuxg7diwff/yx0+FUaJYgjDF+RUSYN28eV1xxBf369WPz5s1Oh1RhFZogRKSdiDwkIs+KyBMi8tfsuaONMcYJ1atX58MPP6RatWoMHDjQrkf4SL53MYnIvcAIYAeQDGwBQoBOwBgR2QhMUNXUMojTGGPO0KhRI/71r39Rv3593KM/m1JW0G2uNYCOqnrcU6GItAGaApYgjDGOaN/eNZeYqvLFF19w4403OhxRxZJvF5OqzsovObjL16nql74JyxhjvPf666/TrVs33njjDadDqVC8uQbRREReEJH3RWRx9suL7eaLyF53V5Sn8s4iclhE1rlfE3OV9RCRLSKyTUQeLVqTjDGBpn///nTp0oW4uDh+/vlnp8OpMLy5i+lDYCcwA3g+16swrwM9Cqnztaq2cb+eABCRysAs4CbgcqCfiFzuxfGMMQEqODiYBX37ckFmJhPHjuVAo0YQH+90WH7Pm6E2Tqjq9KLuWFVXiEhE0UOiPbDNPbMcIpIA3ArY1wJjjGfx8dQdOZL3MzJol5bGw0eO8H9xca6ymBhnY/NjUtjtYSJyN66L0cuAk9nrVfXHQnfuShBLVLWlh7LOwCJgN7AHGKWqm0TkDqCHqsa66w0ArlLVYfkcIw6IAwgLC4tKSEgoLCyP0tLSCA0NLda2/sraXPEFTHs3bIBTpwBI2LiRS0NCuOKSS6BKFWjVyuHgfK8kf+fo6OhkVW3nqcybM4hWwACgC5DlXqfu5ZL4EQhX1TQR6YmrK6sp4Ol+tXyzmKrOBeYCtGvXTjt37lysYJKSkijutv7K2lzxBUx7u3SB7C+7zz1H51GjANcHhwTAMxK++jt7kyBuAy5S1VOleWBVPZLr/Sci8pKI1MN1RtEoV9WGuM4wjDHGs8aNISUlZ1GBscDxmjWxeeiKz5uL1OuB2qV9YBFpIO6nW0SkvTuWA8APQFP33VNVgLuAQu+aMsYEsClToHr1nEUBTgQFMf3oUb7++mvn4vJz3iSIMODfIvJZEW9zfQf4DrhURHaLyCARGSwig91V7gA2ish6YDpwl7pkAMOAz4DNwAJV3VScxhljAkRMDMydC+HhruXwcKbMnk2TJk0YNGgQx4/n+0iXKYA3XUzFmsJJVfsVUj4TmJlP2SfAJ8U5rjEmQMXEuF5JSbBzJzWAVyIiuOGGG/jnP//J008/7XSEfiffM4js7h9VXe7plbuOMcaUR127dmXQoEHMmDGDffv2OR2O3ymoiylRRIaLSOPcK0Wkioh0EZH/A+7xbXjGGFMyzz33HD/88INNV1oMBXUx9QAGAu+ISBPgT1yjuVbG9UzE/6rqOl8HaIwxJVG7dm1q164NwLZt27jkkkucDciPFDRY3wlVfUlVOwLhQFegraqGq+rfLTkYY/zJjBkzaNGihY3VVARezSinqqdV9XdV/dPH8RhjjE/ceeedhIaGEhsbS2ZmptPh+AWbctQYExDOO+88XnzxRb777jtmzvR4A6XJwxKEMSZgxMTE0LNnT8aNG8eOHTucDqfc82Y+iGEicm5ZBGOMMb4kIsyePZvzzjuPbdu2OR1OuefNg3INgB9E5EdgPvCZ2gzhxhg/1ahRI7Zu3UpQkDcff4Gt0DMIVX0M1yir84B7ga0i8qSIXOzj2IwxxieCgoLIyspi1qxZ7NljY4Hmx9u7mBT4w/3KAM4FForIVB/GZowxPpOamsqoUaMYMmQI1inimTfXIEaISDIwFfgWaKWq9wNRwO0+js8YY3wiIiKCJ554go8++oj33nvP6XDKJW/OIOoBfVS1u6q+p6qnAVQ1C7jFp9EZY4wPPfTQQ7Rr145hw4Zx4MABp8Mpd7y5BjFRVVPyKdtc+iEZY0zZCAoKYt68eRw6dIgHH3zQ6XDKHbuMb4wJaJGRkUybNo0WLVo4HUq547MEISLzcXVB7VXVlh7KY4Ax7sU04H5VXe8u2wkcBTKBjPwm1DbGmNIwdOjQnPeqis1k4OLLJ6lfxzUibH52ANeraiQwGZibpzxaVdtYcjDGlAVVZcyYMQwfPtzpUMoNnyUIVV0BHCygfKWqHnIvrgIa+ioWY4wpjIhw+vRpZs2axfLly50Op1wQX97/KyIRwBJPXUx56o0CLlPVWPfyDuAQoMAcVc17dpF72zggDiAsLCwqISGhWLGmpaURGhparG39lbW54gu09kLJ2nz8+HFiY2MREV599VVCQkJKOTrfKEmbo6Ojk/PtqVFVn72ACGBjIXWigc1A3VzrLnD/PA9YD1znzfGioqK0uBITE4u9rb+yNld8gdZe1ZK3+csvv1RAH3nkkdIJqAyUpM3AGs3nM9XR0VxFJBJ4FbhVVXNuQlbVPe6fe4EPgPbORGiMCTRdunTh73//O7NmzQr4eawdSxDuua7fBwao6i+51tcQkZrZ74FuwEZnojTGBKJnn32WNWvWBPw81r68zfUdoDNQT0R2A5OAYABVnQ1MBOoCL7lvKcu+nTUM+MC9Lgh4W1WX+ipOY4zJq1atWtSqVQuArVu30rRpU4cjcobPEoSq9iukPBaI9bB+O9DaV3EZY4y3Zs6cycMPP8zatWsD8kE6m1HOGGPyceedd3LOOecwaNCggJzH2hKEMcbko379+kyfPp3vv/+eGTNmOB1OmbMEYYwxBejXrx+33HIL48ePZ/v27U6HU6YsQRhjTAFEhJdffpmwsDB+/fVXp8MpUzaaqzHGFKJhw4b88ssvATePtZ1BGGOMF4KCgsjMzGTmzJn89ttvTodTJixBGGOMl3bt2sXo0aMDZh5rSxDGGOOliIgIJk+ezOLFi1mwYIHT4ficJQhjjCmCBx54gCuvvJLhw4ezf/9+p8PxKUsQxhhTBEFBQcyfP58///yTBx54wOlwfCqwLskbY0wpaNmyJS+++CLNmzd3OhSfsjMIY4wphvvvv5/OnTsDoG+9BRERUKmS62d8vJOhlRo7gzDGmGJSVR7t1Yu0pUuZlZHhWpmSAnFxrvcxMc4FVwrsDMIYY4pJRMhcsYKXMjJIyl2Qng7jxzsUVemxBGGMMSXwxJEjXIxr7oL03AWpqc4EVIp8liBEZL6I7BURj7PBict0EdkmIj+JSNtcZT1EZIu77FFfxWiMMSVVPTycV4FfgSdzFzRu7ExApciXZxCvAz0KKL8JaOp+xQEvA4hIZWCWu/xyoJ+IXO7DOI0xpvimTKFz9ercBcwEjgFUrw5TpjgbVynwWYJQ1RXAwQKq3Aq8oS6rgNoicj7QHtimqttV9RSQ4K5rjDHlT0wMzJ3L1AsvJBmoER4Oc+f6/QVqcPYupguBXbmWd7vXeVp/VRnGZYwxRRMTQ6NcCeHkyZNUdTCc0uJkghAP67SA9Z53IhKHq4uKsLAwkpKSihVMWlpasbf1V9bmii/Q2gvOt3ny5MmcPn2aJ554osyO6as2O5kgdgONci03BPYAVfJZ75GqzgXmArRr106zH1wpqqSkJIq7rb+yNld8gdZecL7NK1asYNKkSYSEhHD11VeXyTF91WYnb3NdDPzNfTfT1cBhVf0d+AFoKiJNRKQKcJe7rjHGlHsPP/wwYWFhjB492u+HBPflba7vAN8Bl4rIbhEZJCKDRWSwu8onwHZgG/AKMARAVTOAYcBnwGZggapu8lWcxhhTmkJDQ5k0aRJff/01n3zyidPhlIjPuphUtV8h5QoMzafsE1wJxBhj/E5sbCwvvPACU6dO5eabb3Y6nGKzsZiMMaaUBQcHs3DhQsLDw50OpUQsQRhjjA+0bt0agMzMTLKysggODnY4oqKzsZiMMcZHjhw5Qtu2bZk2bZrToRSLJQhjjPGRc845h4YNG/Lkk09y6NAhp8MpMksQxhjjQ0899RSHDx/mqaeecjqUIrMEYYwxPhQZGcmAAQOYPn06u3btKnyDcsQShDHG+NjkyZMBmDlzpsORFI3dxWSMMT7WuHFjEhMTadeundOhFImdQRhjTBno0KEDwcHBnDhxwulQvGYJwhhjysjatWuJiIhgxYoVTofiFUsQxhhTRi677DKCgoIYM2aMXwzkZwnCGGPKSLVq1fjnP//JqlWr+OCDD5wOp1CWIIwxpgzdc889NG/enLFjx5KRkeF0OAWyBGGMMWUoKCiIp556il9++YUlS5Y4HU6B7DZXY4wpY7169eLrr7+mY8eOTodSIDuDMMaYMiYidOrUCREp17e9+jRBiEgPEdkiIttE5FEP5Y+IyDr3a6OIZIpIHXfZThHZ4C5b48s4jTHGCW+99Rbh4eHs37/f6VA88uWUo5WBWcBNwOVAPxG5PHcdVX1WVduoahtgLLBcVQ/mqhLtLvevxw+NMcYLUVFR7N+/nylTpjgdike+PINoD2xT1e2qegpIAG4toH4/4B0fxmOMMeVK8+bNGThwILNmzWLHjh1Oh3MW8dXDGiJyB9BDVWPdywOAq1R1mIe61YHdwCXZZxAisgM4BCgwR1Xn5nOcOCAOICwsLCohIaFY8aalpREaGlqsbf2VtbniC7T2gv+1ed++ffTv35/rrruO8ePHF2sfJWlzdHR0cn69NL68i0k8rMsvG/0F+DZP91JHVd0jIucBn4vIv1X1rOfT3YljLkC7du20c+fOZ5SfPn2a3bt3F3ohqFatWoSEhBRYp6IJtDaHhISQmprK9ddf73QoZSYpKYm8/ycqOn9s848//shzzz3Ha6+9RsOGDYu8va/a7MsEsRtolGu5IbAnn7p3kad7SVX3uH/uFZEPcHVZFXkAk927d1OzZk0iIiIQ8ZSzXI4ePUrNmjWLunu/FkhtVlUOHDhAjRo1nA7FmLOMGTOG/v37Fys5+JIvr0H8ADQVkSYiUgVXElict5KI1AKuBz7Kta6GiNTMfg90AzYWJ4gTJ05Qt27dApODqfhEhLp161K5cmWnQzHmLLVr16ZFixYAHD9+3OFo/stnCUJVM4BhwGfAZmCBqm4SkcEiMjhX1duAZap6LNe6MOAbEVkPrAY+VtWlxY3FkoMB+3dgyr+HH36Y6667jqysLKdDAXz8JLWqfgJ8kmfd7DzLrwOv51m3HWjty9iMMaa8adOmDf/7v//LggULuOuuu5wOx56kNv7p6NGjvPzyy34xZLIx3oqJiSEyMpLx48dz6tQpp8OxBOFLBw4coE2bNrRp04YGDRpw4YUX5iwX9sdfs2YNI0aMKNZxp02bRnp6eqH1br/9dv78889iHaOokpKSuOWWW4q17eOPP85zzz2Xs3zq1CmGDBnC9ddfb91GpkKpXLkyTz/9NNu3b2fOnDlOh2OD9flS3bp1WbduHeD6kAsNDWXUqFE55RkZGQQFef4TtGvXrtjz106bNo3+/ftTvXr1AustWrTorLuYVBVVpVKl8vvdoUqVKrz55ptOh2GMT/To0YPo6GimTZvGkCFDHL2xovx+CvhI586dz3q99NJLAKSnp3ssf/311wHYv3//WWVFde+99/Lwww8THR3NmDFjWL16Nddccw1XXHEF11xzDVu2bAHO/MZ97NgxBg4cyJVXXskVV1zBRx+5bvjKzMxk1KhRtGrVisjISGbMmMH06dPZs2cP0dHRREdHA/DOO+/QqlUrWrZsyZgxY3JiadmyJfv372fnzp00b96cIUOG0LZtW3bt2sWzzz7LlVdeSWRkJJMmTfLYlqVLl9K2bVtat25N165dAfJtT25paWncd999OXEvWrQI4IwHfRYuXMi999571ra//vorPXr0ICoqimuvvZZ///vfgOtho9tvv50rr7ySK6+8km+//bZIfxdjygsRYe7cuXz33XeO33VnZxAO+OWXX/jiiy+oXLkyR44cYcWKFQQFBfHFF18wbty4nA/MbFOmTKFLly7Mnz+fP//8k/bt23PDDTfwxhtvsGPHDtauXUtQUBAHDx6kTp06vPDCCyQmJlKvXj327NnDmDFjSE5O5txzz6Vbt258+OGH9O7d+4xjbNmyhddee42XXnqJZcuWsXXrVlavXo2q0qtXL1asWMF1112XU3/fvn38/e9/Z8WKFTRp0oSDB13POF522WWFtmfy5MnUqlWLDRs2AHDo0CGvf3dxcXHMnj2bpk2b8v333zNkyBC++uorHnjgAR566CE6depEamoq3bt3Z/PmzUX5sxhTblxyySWA64z+xIkTVKtWzZE4Ai5BJCUleVx/9OhRqlevnm85QL169Qos91bfvn1zvhkcPnyYe+65h61btyIinD59+qz6y5YtY/HixTn98CdOnCA1NZUvvviCwYMH53RT1alT56xtf/jhBzp37kz9+vUB10WwFStWnJUgwsPDufrqq3OOt2zZMq644grA9Y1/69atZySIVatWcd1119GkSZMzju1Ne7744gtyD4ly7rnnevFbc8WxcuVK+vbtm7Pu5MmTOfv8+eefc9YfOXIkoB4ENBVPZmYmnTt3plWrVjm9HGUt4BJEeZD7ad4JEyYQHR3NBx98wM6dOz12W6kqixYt4tJLLz1rfWEXab29yyd3TKrK2LFj+cc//lHgfj0d29v2eNo29zpPQ6NkZWVRu3btnOs6ecu+++47x75pGVPaKleuTGRkJHPmzOHBBx+kWbNmZR5DwF2DKG8OHz7MhRdeCJBzrSOv7t27M2PGjJwP+7Vr1wLQrVs3Zs+enTOvbXY3T82aNTl69CgAV111FcuXL2f//v1kZmbyzjvvFDoWUffu3Zk/fz5paWkA/Pbbb+zdu/eMOh06dGD58uU5I1BmH9ub9nTr1o2ZM2fmLGd3MYWFhbF582aysrI8Tuh+zjnn0KRJE9577z3AlWjWr1/vcZ+ekogx/mbChAmEhITw2GOPOXJ8SxAOGz16NGPHjqVjx45kZmaeUZb9jXrChAmcPn2ayMhIWrZsyYQJEwCIjY2lcePGREZG0rp1a95++23A1U9/0003ER0dzfnnn89TTz1FdHQ0rVu3pm3bttx6a0Gjrrs+bO+++246dOhAq1atuOOOO3ISTrb69eszd+5c+vTpQ+vWrbnzzjsLbU+2xx57jEOHDtGyZUtat25NYmIiAE8//TS33HILXbp04fzzz/e4bXx8PPPmzaN169a0aNEi54L99OnTWbNmDZGRkVx++eXMnj3b4/bG+JMGDRowcuRI3nvvPVavXl32AWTf1lgRXlFRUZrXzz//fNY6T44cOeJVvbKycOFC/dvf/ubTY5S3NpeFH3/80ekQylRiYqLTIZS5itbmI0eOaP369bVfv3751ilJm4E1ms9nql2DKIcWL17M+PHjmT9/vtOhGGMcVrNmTT777DOaN29e5se2BFEO9erVi169ejkdhjGmnMi+o/DEiRMEBweX2fMRdg3CGGP8QEpKCs2aNcu51lgWLEEYY4wfaNSoEeeddx6PPfZYoTNklhZLEMYY4wcqVarEM888Q2pqapk9OOfTBCEiPURki4hsE5FHPZR3FpHDIrLO/Zro7bbGGBNounbtSrdu3ZgyZUqZjMTsswQhIpWBWcBNwOVAPxG53EPVr1W1jfv1RBG39QsiwoABA3KWMzIyqF+/frGHvzbGBK6nn36agwcP8tZbb/n8WL48g2gPbFPV7ap6CkgACn5Cq3S2LZn4eIiIgEqVXD/j40u8yxo1arBx48acuWY///zznKeNfSX76WpjTMVyxRVXsHr1aoYOHerzY/nyNtcLgV25lncDV3mo18E99/QeYJSqbirCtohIHBAHrqEa8g6mV6tWrbOeAvYkMzOT4/PmETJ8OJI9aXhKCvr3v3PixAky/vrXQvdRkK5du7Jw4UJ69+7NG2+8QZ8+fVi5ciVHjx7l2LFjPPLII2zatInMzEzGjh3LzTffzObNm7n//vs5ffo0WVlZvPnmmwQHB/PXv/6V77//HnA9QZyWlsa4cePo2bMnV111FatWraJnz5506tSJcePGcezYMerUqcPs2bNp0KABL7/8MvPnz6dy5cpcdtll+Q6JURGpaqkMuOgv0tLSAqq9EDhtXr58OcePH6datWq+a3N+T9CV9AX0BV7NtTwAmJGnzjlAqPt9T2Crt9t6epX4SerwcFU4+xUe7tU+8lOjRg1dv3693n777Xr8+HFt3bq1JiYm6s0336yqqmPHjtU333xTVVUPHTqkTZs21bS0NB02bJi+9dZbqqp68uRJTU9P1x07dmiLFi1y9v3ss8/qpEmTVFX1+uuv1/vvv19VVU+dOqUdOnTQvXv3qqpqQkKC3nfffaqqev755+uJEyf0yJEjeujQoRK1zd/Yk9QVX6C0eenSpVqrVi39+eef/fJJ6t1Ao1zLDXGdJeRQ1SO53n8iIi+JSD1vtvWJ1NSirS+CyMhIdu7cyTvvvEPPnj3PKMtvOO8OHTowZcoUdu/eTZ8+fWjatGmhx8keE2nLli1s3LiRG2+8EXCdIWWPbxQZGUlMTAzdu3enX79+JW6bMabstW3blqysLMaOHcuDDz7ok2P48hrED0BTEWkiIlWAu4DFuSuISANxj0gnIu3d8RzwZlufaNy4aOuLqFevXowaNeqsD2V1D+e9bt061q1bR2pqKs2bN+fuu+9m8eLFVKtWje7du/PVV18RFBREVlZWzrZ574fOHrZbVWnRokXOPjds2MCyZcsA+Pjjjxk6dCjr1q0jKirKrlcY44fq16/PmB49+Oijj9iwcGGpXTPNzWcJQlUzgGHAZ8BmYIGqbhKRwSIy2F3tDmCj+xrEdOAu91mPx219FWuOKVMg7zzO1au71peCgQMHMnHiRFq1anXG+vyG896+fTsXXXQRI0aMoFevXvz000+EhYWxd+9eDhw4wMmTJ1myZInHY1166aXs27eP7777DoDTp0+zadMmsrKy2LVrF9HR0UyePJk///wzZ1hvY4wfiY/nwSVLaAAkrl8PKSkQF1eqScKnYzGp6ifAJ3nWzc71fiYwM+92+W3rczExrp/jx7u6lRo3diWH7PUl1LBhQx544IGz1k+YMIEHH3yQyMhIVJWIiAiWLFnCu+++y1tvvUVwcDANGjRg4sSJBAcHM3HiRK666iqaNGnCZZdd5vFYVapUYeHChYwYMYLDhw+TkZGRM+lI//79OXz4MJmZmTz00EPUrl27VNpnjClD48dT4/hxVgPbbr0VvvkG0tNdn1+l9Jkl2d9aK4J27drpmjVrzli3efNmr0ZBDMTpKQOxzWvXrs0Z+CwQJCUleZzVryILmDZXquS6jQZIeu45Oo8a5VovArm6oQsjIsmq2s7jIUoepTHGmDLn42umYAnCGGP8k4+vmUKAJIiK1I1mis/+HZgKJSYG5s6F8HDXcni4a7mUrj9AACSIkJAQDhw4YB8OAU5VOXDgQL7zZBvjl2JiYOdOiIpy/SzF5AABMKNcw4YN2b17N/v27Suw3okTJwgJCSmjqMqHQGtzSEgIx44dczoMY/xGhU8QwcHBNGnSpNB6SUlJAXV3CwRmm1NSUpwOwRi/UeG7mIwxxhSPJQhjjDEeWYIwxhjjUYV6klpE9gHF7WSuB+wvxXD8gbW54gu09oK1uajCVbW+p4IKlSBKQkTW5Pe4eUVlba74Aq29YG0uTdbFZIwxxiNLEMYYYzyyBPFfc50OwAHW5oov0NoL1uZSY9cgjDHGeGRnEMYYYzyyBGGMMcajgEoQItJDRLaIyDYRedRDuYjIdHf5TyLS1ok4S5MXbY5xt/UnEVkpIq2diLM0FdbmXPWuFJFMEbmjLOPzBW/aLCKdRWSdiGwSkeVlHWNp8+Lfdi0R+ZeIrHe3+T4n4iwtIjJfRPaKyMZ8ykv/80tVA+IFVAZ+BS4CqgDrgcvz1OkJfAoIcDXwvdNxl0GbrwHOdb+/KRDanKveV7jmPb/D6bjL4O9cG/gZaOxePs/puMugzeOAZ9zv6wMHgSpOx16CNl8HtAU25lNe6p9fgXQG0R7YpqrbVfUUkADcmqfOrcAb6rIKqC0i55d1oKWo0Dar6kpVPeReXAU0LOMYS5s3f2eA4cAiYG9ZBucj3rT5buB9VU0FUFV/b7c3bVagpogIEIorQWSUbZilR1VX4GpDfkr98yuQEsSFwK5cy7vd64pax58UtT2DcH0D8WeFtllELgRuA2aXYVy+5M3fuRlwrogkiUiyiPytzKLzDW/aPBNoDuwBNgAPqGpW2YTniFL//Krw80HkIh7W5b3H15s6/sTr9ohINK4E0cmnEfmeN22eBoxR1UzXl0u/502bg4AooCtQDfhORFap6i++Ds5HvGlzd2Ad0AW4GPhcRL5W1SM+js0ppf75FUgJYjfQKNdyQ1zfLIpax5941R4RiQReBW5S1QNlFJuveNPmdkCCOznUA3qKSIaqflgmEZY+b/9t71fVY8AxEVkBtAb8NUF40+b7gKfV1UG/TUR2AJcBq8smxDJX6p9fgdTF9APQVESaiEgV4C5gcZ46i4G/ue8GuBo4rKq/l3WgpajQNotIY+B9YIAff5vMrdA2q2oTVY1Q1QhgITDEj5MDePdv+yPgWhEJEpHqwFXA5jKOszR50+ZUXGdMiEgYcCmwvUyjLFul/vkVMGcQqpohIsOAz3DdATFfVTeJyGB3+Wxcd7T0BLYB6bi+gfgtL9s8EagLvOT+Rp2hfjwSppdtrlC8abOqbhaRpcBPQBbwqqp6vF3SH3j5d54MvC4iG3B1v4xRVb8dBlxE3gE6A/VEZDcwCQgG331+2VAbxhhjPAqkLiZjjDFFYAnCGGOMR5YgjDHGeGQJwhhjjEeWIIwxxnhkCcKYEhCR2iIypIDyaiKyXEQqF1DnCxE51zcRGlN8liCMKZnaQL4JAhiIa5C8zALqvFnIPoxxhCUIY0rmaeBi9zwLz3ooj8H1FDMicr6IrHDX3Sgi17rrLAb6lVG8xnjNHpQzpgREJAJYoqotPZRVAVJVtYF7eSQQoqpT3F1O1VX1qLtsK3B1BRgLy1QgATPUhjEOqAf8mWv5B2C+iAQDH6rqulxle4ELAEsQptywLiZjfOc4EJK94J7w5TrgN+DNPHMyhLjrG1NuWIIwpmSOAjU9Fbhn6qssIiEAIhIO7FXVV4B5uKaPxD3jWQNgZ1kEbIy3LEEYUwLuawbfui86e7pIvYz/TsLUGVgnImuB24EX3eujgFWq6rfTYZqKyS5SG+NDInIF8LCqDiigzovAYlX9suwiM6ZwdgZhjA+p6logsaAH5YCNlhxMeWRnEMYYYzyyMwhjjDEeWYIwxhjjkSUIY4wxHlmCMMYY45ElCGOMMR79P60GbIUlNIgwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Utilisation de la méthode programmée\n", "\n", "ALPHA , alpha = grad_pas_fixe_reg_pol(t , y , d = 2 , rho = 1e-2 , Max_iter = 10**5 , tol = 1e-8)\n", "print(\"- Méthode programmée:\")\n", "print(\"Trajectoire: y(t) = \" , alpha[2] , \"t^2 + \" , alpha[1] , \"t + \" , alpha[0])\n", "print(\"Vitesse initiale: v_0 = \" , alpha[1]/np.sin(np.pi/4) , \"m/s\")\n", "print(\"Accélération de la pesanteur calculée: g = \" , -2*alpha[2] , \"m/s**2\")\n", "print(\" \")\n", "\n", "\n", "# Utilisation de la méthode np.polyfit\n", "\n", "alpha_bis = np.polyfit(t , y , 2)\n", "print(\"- Méthode np.polyfit:\")\n", "print(\"Trajectoire: y(t) = \" , alpha_bis[0] , \"t^2 + \" , alpha_bis[1] , \"t + \" , alpha_bis[2])\n", "print(\"Vitesse initiale: v_0 = \" , alpha_bis[1]/np.sin(np.pi/4) , \"m/s\")\n", "print(\"Accélération de la pesanteur calculée: g = \" , -2*alpha_bis[0] , \"m/s**2\")\n", "\n", "plt.figure()\n", "plt.scatter(t , y , color = \"red\" , label = \"Mesures\")\n", "plt.plot(t , alpha[0] + alpha[1]*t + alpha[2]*t**2 , linestyle = \"dashed\" , color = \"black\" , label = \"Trajectoire calculée\")\n", "plt.xlabel(\"t (s)\")\n", "plt.ylabel(\"y (m)\")\n", "plt.legend()\n", "plt.title(\"Trajectoire du ballon\")\n", "plt.grid()\n", "plt.show()" ] }, { "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }