{ "cells": [ { "cell_type": "markdown", "id": "c64ac0da", "metadata": {}, "source": [ "# Nuage d'ions" ] }, { "cell_type": "markdown", "id": "f77944b4", "metadata": {}, "source": [ "On considère un nuage de $N$ ions de masse $m$ et charge $q$. On note $\\vec{r_i}$ et $\\vec{v_i}$ la position et la vitesse du ième ion. Les ions sont dans un piège électrostatique. De plus, ils interagissent entre eux par la force de Coulomb. \n", "\n", "La force électrostatique dérive d'un potentiel électrostatique :\n", "\n", "$$V(\\vec{r}) = \\frac12 k_x r_x^2 +\\frac12 k_y r_y^2 + \\frac12 k_z r_z^2$$\n", "\n", "La force de Coulomb s'écrira sous la forme:\n", "\n", "$$\\vec{f_i}(\\vec{r_i}, \\vec{r_j}) = \\kappa \\frac{q^2}{|\\vec{r_i} - \\vec{r_j}|^3}(\\vec{r_i} - \\vec{r_j})$$\n", "\n", "Pour simplifier et avoir une meilleure représentation graphique, on suprimera totalement la dimension $z$. \n", "\n", "On utilisera des unités adimensionnées avec des constantes de l'ordre de 1. Par défaut on prendra : $m=9$, $q=1$, $k_x=1$, $k_y=1.3$, $\\kappa=1$.\n", "\n", "L'état des $N$ particules à un instant donné est représenté par 4 tableaux numpy de taille $N$: ``r_x``, ``r_y``, ``v_x``, ``v_y``. Toutes les fonctions seront écrites avec ces variables. On ne regroupera les 4 tableaux que dans la fonction ``f(t, y)`` qui sera utilisée par ``solve_ivp``.\n", "\n", "Le tableau ``y`` est défini en rassemblant les 4 tableaux numpy. On définit alors les fonctions ``join`` et ``split`` (voir ci-dessous)\n", "\n", "Les paramètres seront des constantes globales.\n", "\n", "1. Écrire la fonction ``force_piege(r_x, r_y)`` qui renvoie la force dérivant du potentiel (``f_x`` et ``f_y``).\n", "\n", "2. Écrire la fonction ``force_coulomb(r_x, r_y)`` qui renvoie $f_x$ et $f_y$, la force de Coulomb. On ne cherchera pas à éviter les boucles for dans cette fonction.\n", " \n", "3. Écrire la fonction ``f(t, y)`` qui définit la dynamique du problème.\n", "\n", "4. On condidère un nuage avec une distribution initiale donnée par une loi de Maxwell-Boltzmann (avec $k_B=1$ et $T=1$) pour les particules sans interaction. Calculer jusqu'au temps $T=20$ l'évolution de la position des particules.\n", "\n", "5. Vérifier que l'énergie totale est conservée. On utilisera les fonctions ci-dessous.\n", "\n", "6. Faire une animation (voir exemple de code ci dessous)\n", "\n", "7. On rajoute une force de dissipation, selon l'axe $x$ : \n", "\n", " $$ \\vec F_{i,x} = -\\alpha v_{i, x} $$\n", " \n", " Simuler l'experience en prenant $\\alpha = 0.1$ et $N=20$. Que se passe-t-il au temps long ?\n", " \n", "8. Utiliser le decorateur ``jit`` de numba pour accélérer le code. Faire une comparaison de vitesse.\n", "\n", "9. On considère un nuage contenant deux types de particules (masse $m_1$ et $m_2$) et de même charge $q=1$. Seuls les particules de type $1$ sont soumises à la force dissipative. Simuler l'expérience en prenant les mêmes paramètres qu'à la question 7. On prendra par exemple $m_1=9$, $N_1=20$ et $m_2=2$ et $N_2=5$" ] }, { "cell_type": "markdown", "id": "501f2f54", "metadata": { "tags": [ "sol" ] }, "source": [ "--------------------------------------------------------------------------------------\n", "\n", "## Solution" ] }, { "cell_type": "code", "execution_count": 1, "id": "11f9d33a", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "import numpy as np\n", "from numpy import pi\n", "\n", "from scipy.integrate import solve_ivp" ] }, { "cell_type": "code", "execution_count": 2, "id": "7ee13f2b", "metadata": { "tags": [ "keep" ] }, "outputs": [], "source": [ "# Paramètres\n", "m = 9; q = 1\n", "k_x, k_y = 1, 1.3\n", "kappa=1\n", "k_B = 1\n", "T = 1" ] }, { "cell_type": "markdown", "id": "06ac87bd", "metadata": {}, "source": [ "**Quelques fonctions pour vous aider :**\n", "\n", "* Fonction pour passer de 4 tableaux à un seul et réciproquement" ] }, { "cell_type": "code", "execution_count": 3, "id": "1ce15680", "metadata": { "tags": [ "keep" ] }, "outputs": [], "source": [ "\n", "def join(r_x, r_y, v_x, v_y):\n", " return np.concatenate((r_x, r_y, v_x, v_y))\n", "\n", "def split(y):\n", " N = len(y)//4\n", " return y[:N], y[N:2*N], y[2*N:3*N], y[3*N:4*N]" ] }, { "cell_type": "markdown", "id": "f5e4a5f5", "metadata": {}, "source": [ "* Calcul de l'énergie" ] }, { "cell_type": "code", "execution_count": 4, "id": "70f02398", "metadata": { "tags": [ "keep" ] }, "outputs": [], "source": [ "def energie_cinetique(r_x, r_y, v_x, v_y):\n", " return np.sum(.5*m*v_x**2 + .5*m*v_y**2)\n", "\n", "def energie_piege(r_x, r_y):\n", " return np.sum(k_x*r_x**2/2 + k_y*r_y**2/2)\n", "\n", "def energie_coulomb(r_x, r_y):\n", " N = len(r_x)\n", " total = 0\n", " for i in range(N-1):\n", " for j in range(i+1, N):\n", " d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2\n", " total += kappa*q**2/np.sqrt(d2)\n", " return total\n", "\n", "def energie_totale(r_x, r_y, v_x, v_y):\n", " return (energie_coulomb(r_x, r_y) +\n", " energie_piege(r_x, r_y) +\n", " energie_cinetique(r_x, r_y, v_x, v_y))\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "b13b87d4", "metadata": {}, "outputs": [], "source": [ "# Question 1\n", "def force_piege(r_x, r_y):\n", " return -k_x*r_x, -k_y*r_y\n", "\n", "# Question 2\n", "def force_coulomb(r_x, r_y):\n", " N = len(r_x)\n", " f_x = np.zeros(N)\n", " f_y = np.zeros(N)\n", " for i in range(N):\n", " tot_x, tot_y = 0, 0\n", " for j in range(N):\n", " if i!=j:\n", " d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2\n", " coef = kappa*q**2/(np.sqrt(d2)**3)\n", " tot_x += coef*(r_x[i]-r_x[j])\n", " tot_y += coef*(r_y[i]-r_y[j])\n", " f_x[i] = tot_x\n", " f_y[i] = tot_y\n", " return f_x, f_y\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "84deaffe", "metadata": {}, "outputs": [], "source": [ "# Question 3\n", "def force_totale(r_x, r_y, v_x, v_y):\n", " piege = force_piege(r_x, r_y)\n", " coulomb = force_coulomb(r_x, r_y)\n", " return piege[0]+coulomb[0], piege[1]+coulomb[1]\n", "\n", "def f(t, y):\n", " r_x, r_y, v_x, v_y = split(y)\n", " f_x, f_y = force_totale(r_x, r_y, v_x, v_y)\n", " return join(v_x, v_y, f_x/m, f_y/m)" ] }, { "cell_type": "code", "execution_count": 7, "id": "27d25e42", "metadata": {}, "outputs": [], "source": [ "# Question 4\n", "N = 50\n", "\n", "r_x_0 = np.random.normal(size=N)\n", "r_y_0 = np.random.normal(size=N)\n", "v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "\n", "t_max = 20\n", "t_eval = np.linspace(0, t_max, 51)\n", "res = solve_ivp(f, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0), \n", " t_eval=t_eval)\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "64035f5e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD6CAYAAABNu5eFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjAUlEQVR4nO3de3wc9Xnv8c+zq7tk627hmywkjMGEYoO42MSUhEsIacIloUCgJA0p5TTJSZqmJT05p3F6evJK0ubSlhTqUkouwKGF0tASQqAJuAFjkI0B27LBMr7I2JZs2ZJt3bVP/9iRUWTJkuzdnZX2+35ZL+3+Zmbn8ezqu7O/mfmtuTsiIpIZImEXICIiqaPQFxHJIAp9EZEMotAXEckgCn0RkQyi0BcRySBjhr6Z3W9mLWa2fkjbDWa2wcxiZlY/pL3GzLrMbF3wc+8Ij/fE0McSEZHUyRrHPA8AdwM/HNK2Hrge+PsR5m9y90UjPZCZXQ8cnkiBFRUVXlNTM5FFREQy3po1a/a5e+Xw9jFD391XmlnNsLZGADMbdwFmVgR8EbgD+OfxLldTU0NDQ8O41yMiImBm20dqT0af/qlm9qqZPW9my4a0/1/g20BnEtYpIiLjMJ7unYnYDVS7+34zOw/4NzM7C6gF6tz9D4d/ahiJmd1B/BMB1dXVCS5RRCRzJXRP39173H1/cHsN0AScDiwB6s1sG/Ar4HQze+44j7PC3evdvb6y8pguKREROUEJDX0zqzSzaHC7FpgPbHX3e9x9lrvXAO8F3nT3SxO5bhERGduY3Ttm9jBwKVBhZs3AV4E24G+BSuBJM1vn7h8ALgH+3Mz6gBhwp7u3Jat4ERGZmPGcvXPzKJMeH2Hex4DHxni8bcB7xlOciIgklq7IFRHJIFM29H/00nb+663WsMsQEUkrUzL0+wZiPLR6B5+4/2X+7rkt6NvBRETipmToZ0cjPHrnEj549ky+9bPN3PnjNRzq7gu7LBGR0E3J0AcozM3i7psX878/dCbPNrZwzfdfYEvLobDLEhEJ1ZQNfYiPDfTpZbX8+PYL6ejq45q7X+CpN3aHXZaISGimdOgPWlJXzr9/7r2cfso0/seDa7nr0dfZtKcj7LJERFLO0v0gZ319vSdqlM2e/gG++dRmfrx6O739Mc6bV8otF1Zz9dkzycuOJmQdIiLpwMzWuHv9Me2ZFPqDDhzp5bG1zTy0egdb9x2hpCCbj507h09eXMOc0oKErktEJAwK/RG4O6u27ufB1Tt4ev0eqqbn8fwfX0pWNCN6vURkChst9DM63cyMpXUVfP/j53L3xxez62AXv9ysC7pEZOrK6NAf6rIzq5gxLZeHVo/4ZTMiIlOCQj+QHY1w0/lzee7NVna26cu9RGRqUugPceMF1RjwyCs7wy5FRCQpFPpDzC7J59IFM3ikYSd9A7GwyxERSTiF/jC3XFhN66Eent24N+xSREQSTqE/zKULZjCrOI+HXt4RdikiIgmn0B8mGjFuPL+a/3prH9v3Hwm7HBGRhFLoj+DG8+cSjRgPv6wDuiIytSj0R3BKcR6XnTGDf2nYSW+/DuiKyNSh0B/Fxy+sZv+RXp7esCfsUkREEkahP4pL5lcypzSfh1brgK6ITB0K/VFEIsbNF1Szaut+mloPh12OiEhCjBn6Zna/mbWY2fohbTeY2QYzi5lZ/ZD2GjPrMrN1wc+9QXuBmT1pZpuC5b6RnP9OYt1QP4esiPGw9vZFZIoYz57+A8BVw9rWA9cDK0eYv8ndFwU/dw5p/yt3PwNYDFxsZh88kYJTaca0PK48q4pH1zbT3TcQdjkiIidtzNB395VA27C2RnffPN6VuHunu/8yuN0LrAXmTLDWUNx60TwOdvbx76+9E3YpIiInLRl9+qea2atm9ryZLRs+0cxKgA8D/znaA5jZHWbWYGYNra3hjm+/pLac+TOK+MGqbaT7F86IiIwl0aG/G6h298XAF4GHzGz64EQzywIeBv7G3beO9iDuvsLd6929vrKyMsElToyZcdvSGtbv6uDVnQdDrUVE5GQlNPTdvcfd9we31wBNwOlDZlkBvOXu30vkepPt+sWzmZabxQ9f3BZ2KSIiJyWhoW9mlWYWDW7XAvOBrcH9vwCKgS8kcp2pUJibxUfPm8OTb+ym9VBP2OWIiJyw8Zyy+TCwClhgZs1mdruZXWdmzcAS4EkzezqY/RLgdTNbBzwK3OnubWY2B/gKsBBYG5zO+elk/IeS5daL5tE34Dzyik7fnMxWvtmqgfQko2WNNYO73zzKpMdHmPcx4LER2psBm3B1aeS0GUUsm1/Bj1/awZ2/WUdWVNe1TTYb3+ngk//0MnNKC3jq88sozB3z5S8y5Si5JuC2JTXs6ejmGX3ByqTj7iz/9w0U5mSx80An33hqU9gliYRCoT8B7z9jBrNL8vnBqm1hlyIT9OQbu3n57Ta+fPUZfOriU/nRS9t5Ycu+sMsSSTmF/gREI8atF83jpa1tbN5zKOxyZJy6egf4+pONLJw5nZvOr+aPP7CA2spC/uTR1znU3Rd2eSIppdCfoBvPn0tOVoQfam9/0rjn+Sbeae9m+UfOIhox8rKjfPuGc9jd3sVf/Edj2OWJpJRCf4LKCnP4yDmzePzVXXRoLzHtNR/o5O+fb+LD58ziglPLjrYvri7l93+zjkcadvLLTS0hViiSWgr9E/CJJTV09g7w2JrmsEuRMXz9p42YwZ9+8Ixjpn3h8vksqJrGXY+9zsHO3hCqE0k9nbN2As6eU8yiuSXc/YstPP/myGMDGfEhHCIWvxcxiJgRjRpZESMrEon/jhrR+EwMDu3j+NHbFiwXsfg8keAx+2NOfyxG/4DTNxDcjjkRM7KDx82ORsiOxtcTjcYfIxosH4nE73f1DXCws5e2I70cONLHgc5eDnT20tsfIz8nSn52lLzs6NHbxfnZVJcXMK+skJryAqrLC5hZnH/0/5As7o7ZxNbxYtM+fvrGHr54xenMKsk/ZnpuVpRv//Y5XPv9F1j+xAa+d9PiRJUrkrYU+ifoj648nW///E0OHBl5D9GBmMfDO+bx0Iq5MxDzeGAHQT14H959o3j3NsHyTiz4TfA7ErGjgZ4djRx983CHvoHY0cfv7Y+/GQzE/OjjDMTeHTguK2KUFuZQWpBNaUEOp80ooqQgh5yo0d0Xo7t/gK7eAbr6BujuG2Dz3kP8Z2MLvQPvfndwTjRC5bRccrIG38jefUPLihiGEfzDjPh9YCDYHoM/8e0So3cgXndPf4yevgF6+mMMuDN/RhGL5pawuLqURXNLOL1q2qhvNv0DMb72xEbmlOZzxyW1oz6P75ldzGfffxrfe/YtSgpyKCvMebfO4LlYWlfO4urS8b0wRNKcQv8ELZtfybL54Q4GdzI8CNxoxCa8Bz0Qc3a3d7Fjfyfb2zrZtv8IrR099AWh3TfgDASfPPoH/OgnFwc8Bk78DSMaMbKzI0SDTz3RSPyTSG52hNysCDlZEXKzouRmRXCgcXcHP9+4l39uiHerFeREOXt2MTXlhVQV53HK9DxOKc6lanoeL2zZx+a9h7j31nPJy44e9//zmfedxstvt/HAKGMr5WZFePTOpZw9p3hC20kkHVm6DxdcX1/vDQ0NYZchacLd2b6/k3U7D/LqjgO8vqudXQe6aD3cw/CX8tK6ch789IXjelPz4FOQu8ffnDzezXbgSB8fvedFYu785LMXM2NaXnL+YyIJZmZr3L3+mHaFvkwF/QMxWg/3sLu9m73t3bQe7uGKhVXMLD62L3+iNrzTzsfuWcXCWdN56PcuJDfr+J8cRNLBaKGv7h2ZErKiEWYW5yck5Ic7a1Yxf3XDOXzmobX8n39bzzc/+hsT7hITSRc6ZVNkHD70GzP53PtP458bmkft+xeZDBT6IuP0h5efzhULq/iLJxv51Vsat0cmJ4W+yDhFIsZ3b1xEXWUhn3loLdv2aVx+mXx0IFdkgnbs7+Qj3/8V2dEIv39JLTdfUJ2QsfkP9/Szp72LPe097G7vou1ILzlZkfjFcUMuksvLil+XMXjRXjTy7u/8YJ6C4GK6SJIvmpP0pbN3RBLojeZ2/t9PN/LS1jZKCrL5xJIaPrm0htLCnF+b71B3H69sa+PFLfvZvPcQ/QODF6LFL8wbcKerd4CWjh4O9fQnvM687AiFOVnUzShi8dwSFs0tYVF1SVIOeI9Xe2cfHd19zC0rCK2GTKDQF0mCtTsOcM9zTTyzcS/52VFuvqCaZadX0LCtjReb9vN6czsDMScnK8IZp0wjLysavwgt+MmKGDlZEaqm53FKcR4zj15klkd5US79AzG6+uJXRXf3xY5eGd0/eIX1kKut+2Men6e3n87egeCnn8M9/WzcfYjGdzqOXkldNT2XRXNLqJqeR152/AK4wd+52VGm52VRGlyhXFqYQ1lBDvk5Ez9VNRZztu47wtrtB1iz/QBrdxzgrZbDANRWFHLFWVVcubCKxXNLE/KpxN052NnH/iM95OdkUVqQTX52NCPPtlLoiyTRW3sPcc/zTfxk3TtHr3Q+Z04xS+sqWFpXzrnzSse8MjjZevoH2PhOB6/tPMi6nQd5vbmdts5eeoLhNsaKgrzs+BtDLObHDA9iBtnRCDnBkCCDt/cf6aW9Kz4abXF+NudWl3DevFIKc7P4xaYWVjXtpz/mVBTlcvmZM1hSV052NBJ/UxscoiO4erxvyPAcvf0x+oI3xJZDPbR0dLOno5u9HT309sd+re6caITigmxKC7IpL8zlSx9YwHnzpv6wGgp9kRTYdbCLt1uPsKi6hKJJ9B287vGB+3r6458oOrr7OHAkGIivs5e2YDC+3v7Y0fGTBgfuGxwjqm8gdnTcp95gOI6i3CiLq0s5t7qU2orCY/bm27v6eG5zC89s3Mtzm1s5PIEurqzguxEqp+UyY1oupxTnUTU9/lNRlEN33wAHOvs42NnHwc5eDnb28erOA/QPOE987r3MHmEQvqlEoS8iaa2nf4C39x3BMKIRjjlInROMx5QTfIo4ke6gptbDXHv3C1SXF/DonUtPqMtqstAVuSKS1nKzopxxyvSkrqOusoi/vnkRt/+ggbsee52/vmlRxvX36zx9Ecko7z+jii9duYAnXnuHFSu3hl1Oyo0Z+mZ2v5m1mNn6IW03mNkGM4uZWf2Q9hoz6zKzdcHPvUOmnWdmb5jZFjP7G8u0t1cRSRt/cGkdHzp7Jt/82aZRvwhpqhrPnv4DwFXD2tYD1wMrR5i/yd0XBT93Dmm/B/g9YH7wM/wxRURSwsz4yxt+g9OrpvG5DLu6eszQd/eVQNuwtkZ33zzelZjZTGC6u7/k8SPHPwSunWCtIiIJU5CTxT/cVk80YvzeDxsmdObQZJaMPv1TzexVM3vezJYFbbOBod8i3hy0jcjM7jCzBjNraG3NrI9eIpI6c8sK+P7Hz2XrviN86p9e4VB3X9glJV2iQ383UO3ui4EvAg+Z2YQPx7v7Cnevd/f6ysrJ+5WEIpL+lp5WwfduXMTaHQe45b7VtI3yvddTRUJD39173H1/cHsN0AScDuwC5gyZdU7QJiISug+fM4sVt53Hpj2HuPHvV7G3ozvskpImoaFvZpVmFg1u1xI/YLvV3XcDHWZ2UXDWzm3ATxK5bhGRk/H+M6p44HfP552DXdxw7yp2tnWGXVJSjOeUzYeBVcACM2s2s9vN7DozawaWAE+a2dPB7JcAr5vZOuBR4E53HzwI/AfAfcAW4p8Ankrsf0VE5OQsravgx5++kPauPm64dxVbgsHhphINwyAiMsymPR3cet/LxNz50e0XcNas4rBLmrDRhmHQFbkiIsOcccp0/uXOJeRlRfj0DxrYd7gn7JISRqEvIjKCUysKWXFbPW1HevnsQ2vpH4iNvdAkoNAXERnFe2YX842Pns1LW9v4xlObwi4nITTKpojIcVy3eA6v7Wznvl+9zdlzirlm0ajXlU4K2tMXERnDVz50JhfUlHHXY6/TuLsj7HJOikJfRGQM2dEId9+ymOL8bH7/R2to75y8wzUo9EVExmHGtDz+7pbz2N3execfeZVYLL1Pdx+NQl9EZJzOm1fKVz98Fs9tbuXelU1hl3NCFPoiIhNwy4XVLK0r57E1zWPPnIYU+iIiE2BmvG/BDJpaj7CnffINzKbQFxGZoCV15QCs2rov5EomTqEvIjJBC2dOpzg/mxe37A+7lAlT6IuITFAkYiypLefFpv2k+6CVwyn0RUROwMWnlbPrYBc7Jtm4+wp9EZETsKSuAoAXmyZXF49CX0TkBNRVFjJjWq5CX0QkE5gZS+vKWdW0b1L16yv0RURO0NLTKth3uJc3906er1VU6IuInKClwfn6LzZNnvP1FfoiIidoTmkB1WUFvDCJztdX6IuInISLTytn9db9k+brFBX6IiInYUldBYd6+tnwzuT4cpVxhb6Z3W9mLWa2fkjbDWa2wcxiZlY/wjLVZnbYzL40pO0Pg2XWm9nDZpaXmP+GiEg4ltQO9utPji6e8e7pPwBcNaxtPXA9sHKUZb4DPDV4x8xmA/8TqHf39wBR4KaJFCsikm4qp+VyelXRpDmYO67Qd/eVQNuwtkZ33zzS/GZ2LfA2sGHYpCwg38yygALgnYkWLCKSbpbWVfDKtjZ6+gfCLmVMCe/TN7Mi4C7ga0Pb3X0X8FfADmA30O7uP0/0+kVEUm1pXTndfTHW7TgYdiljSsaB3OXAd939165WMLNS4BrgVGAWUGhmt470AGZ2h5k1mFlDa2trEkoUEUmcC2vLidjk6NdPRuhfCHzLzLYBXwD+l5l9FrgceNvdW929D/hXYOlID+DuK9y93t3rKysrk1CiiEjiFOdn857ZxZOiXz/hoe/uy9y9xt1rgO8BX3f3u4l361xkZgVmZsBlQGOi1y8iEoaldRW8uuMgnb39YZdyXOM9ZfNhYBWwwMyazex2M7vOzJqBJcCTZvb08R7D3VcDjwJrgTeCda84qepFRNLE0rpy+mPOK9sOhF3KcWWNZyZ3v3mUSY+PsdzyYfe/Cnx1XJWJiEwi9TWlZEeNF7fs4zdPT99uaV2RKyKSAAU5WZwzp4RXtrWNPXOIFPoiIgly2oyitP/6RIW+iEiCVJcXsO9wL4d70vdgrkJfRCRB5pUVArBjf/ru7Sv0RUQSpLqsACCtu3gU+iIiCVJdPhj6R0KuZHQKfRGRBCnOz6Y4P5vt6t4REckM88oL1L0jIpIpqssU+iIiGaO6rIBdB7rS9jtzFfoiIgk0r7yA/pizu7077FJGpNAXEUmg6uBc/XQ9mKvQFxFJoMHTNren6WmbCn0RkQQ6ZXoeOdFI2h7MVeiLiCRQNGLMKctP26EYFPoiIglWXVagPn0RkUwxr6yAnW2duHvYpRxDoS8ikmDV5YUc6unnQGdf2KUcQ6EvIpJgg6Ntbt+ffmfwKPRFRBJsXnn6DrGs0BcRSbC5pUHop+HBXIW+iEiC5edEmTEtV3v6IiKZYl55AdsnY+ib2f1m1mJm64e03WBmG8wsZmb1IyxTbWaHzexLQ9pKzOxRM9tkZo1mtiRx/w0RkfQyt6xg0nbvPABcNaxtPXA9sHKUZb4DPDWs7a+Bn7n7GcA5QOP4yxQRmVzmlRWyp6Ob7r6BsEv5NWOGvruvBNqGtTW6++aR5jeza4G3gQ1D2oqBS4B/DJbvdfeDJ1y1iEiaGzyDp/lAeu3tJ7RP38yKgLuArw2bdCrQCvyTmb1qZveZWeFxHucOM2sws4bW1tZEligikhJzj56rP4VDH1gOfNfdDw9rzwLOBe5x98XAEeDLoz2Iu69w93p3r6+srExwiSIiyZeu5+pnJfjxLgQ+ZmbfAkqAmJl1A48Cze6+OpjvUY4T+iIik115YQ6FOdG029NPaOi7+7LB22a2HDjs7ncH93ea2YLgWMBlwMZErltEJJ2YWfwMnsm2p29mDwOXAhVm1gx8lfiB3b8FKoEnzWydu39gjIf6HPCgmeUAW4HfPZnCRUTS3bzyAppa02v8nTFD391vHmXS42Mst3zY/XXAMef0i4hMVfPKC/nl5lZiMScSsbDLAXRFrohI0swtK6C3P8beQ91hl3KUQl9EJEnmlaXfwGsKfRGRJBk8bTOdxuBR6IuIJMmsknyiEdOevohIJsiORphVkpdWp20q9EVEkqi6LL2GWFboi4gkUXVZITvS6LtyFfoiIkk0r7yAA519dHT3hV0KoNAXEUmq6jQ7bVOhLyKSREdDP0369RX6IiJJVJ1mQywr9EVEkmh6XjalBdlpM8SyQl9EJMmqywvZ0ZYeZ/Ao9EVEkqy2opDNew4Ti3nYpSj0RUSSbdn8CvYd7uGNXe1hl6LQFxFJtvctmEE0YjzbuDfsUhT6IiLJVlqYQ/28Up7ZqNAXEckIVyysYtOeQ+wM+dRNhb6ISApcfmYVQOhdPAp9EZEUqKkoZP6MotC7eBT6IiIpcvnCKla/3UZ7Z3iDryn0RURS5IqFVQzEnOfebAmtBoW+iEiKLJpTQkVRTqhdPOMKfTO738xazGz9kLYbzGyDmcXMrH6EZarN7LCZfWlYe9TMXjWz/zj58kVEJo9IxLjsjCqe39xKb38snBrGOd8DwFXD2tYD1wMrR1nmO8BTI7R/Hmgc53pFRKaUKxZWcainn9Vv7w9l/eMKfXdfCbQNa2t0980jzW9m1wJvAxuGtc8BPgTcdyLFiohMdhefVkFedoRnQ+riSXifvpkVAXcBXxth8veAPwGO+7nGzO4wswYza2htbU10iSIiocnPibJsfiXPNrbgnvoB2JJxIHc58F13Pzy00cx+C2hx9zVjPYC7r3D3enevr6ysTEKJIiLhueLMKnYd7GLj7o6UrzsrCY95IfAxM/sWUALEzKwbmA18xMyuBvKA6Wb2Y3e/NQk1iIikrfedMQMzeHZjC2fNKk7puhO+p+/uy9y9xt1riHfnfN3d73b3P3X3OUH7TcAvFPgikokqp+VybnUpzzTuSfm6x3vK5sPAKmCBmTWb2e1mdp2ZNQNLgCfN7OlkFioiMpVcfmYV63d1sLu9K6XrHVf3jrvfPMqkx8dYbvko7c8Bz41n3SIiU9EVC6v45s828WxjC79z0byUrVdX5IqIhKCuspBTKwpTfnWuQl9EJARmxqULKnlp634GUvjduQp9EZGQLKiaRm9/jF0HUtevr9AXEQlJ3YwiAJr2HR5jzsRR6IuIhKS2ohCAphaFvojIlFdWmENxfjZb9x1J2ToV+iIiITEz6ioL2dqqPX0RkYxQW1lEU6v29EVEMkJtZSGth3o41J2a781V6IuIhKiuMn4Gz9YU7e0r9EVEQlRXGZzBk6J+fYW+iEiIqssKiUZMe/oiIpkgJytCdVkBW1N0gZZCX0QkZLUVhTS1aE9fRCQj1M0o4u39R1Iy8JpCX0QkZLUVhSkbeE2hLyISslQOvKbQFxEJ2eDAa6k4g0ehLyISssGB11Jxrr5CX0QkZKkceE2hLyKSBmori9S9IyKSKeoqi2hJwcBrY4a+md1vZi1mtn5I2w1mtsHMYmZWP8Iy1WZ22My+FNyfa2a/NLONwXKfT+x/Q0RkcqutTM3B3PHs6T8AXDWsbT1wPbBylGW+Azw15H4/8EfuvhC4CPiMmS2cWKkiIlNXqgZeyxprBndfaWY1w9oaIX7wYTgzuxZ4GzgyZP7dwO7g9iEzawRmAxtPvHQRkakjVQOvJbRP38yKgLuArx1nnhpgMbA6kesWEZnMUjXwWqIP5C4HvuvuI1YdvCk8BnzB3TtGexAzu8PMGsysobW1NcElioikp7rK5A+8Nmb3zgRdCHzMzL4FlAAxM+t297vNLJt44D/o7v96vAdx9xXACoD6+vrkj0AkIpIGaiuLWPnWPgZiTjRybPd5IiQ09N192eBtM1sOHA4C34B/BBrd/TuJXKeIyFQxOPDaOwe7mFtWkJR1jOeUzYeBVcACM2s2s9vN7DozawaWAE+a2dNjPMzFwO8A7zezdcHP1SddvYjIFDI48NqWJJ7BM56zd24eZdLjYyy3fMjtXwHJ+awiIjJFDB147X0LkrMOXZErIpImygpzKCnITuoYPAp9EZE0YWbxr05U6IuIZIZkD7ym0BcRSSPJHnhNoS8ikkaSPfCaQl9EJI0MDryWrOEYFPoiImlkcOC1ZA3HoNAXEUkjyR54TaEvIpJm4t+Xqz19EZGMUFtZxNZ9RxiIJX68yUSPsikiIifp/JoyWjq66eztZ1pedkIfW6EvIpJmrlhYxRULq5Ly2OreERHJIAp9EZEMotAXEckgCn0RkQyi0BcRySAKfRGRDKLQFxHJIAp9EZEMYu6Jv8w3kcysFdh+gotXAPsSWE6iqK6JUV0To7omZqrWNc/dK4c3pn3onwwza3D3+rDrGE51TYzqmhjVNTGZVpe6d0REMohCX0Qkg0z10F8RdgGjUF0To7omRnVNTEbVNaX79EVE5NdN9T19EREZYkqEvpldZWabzWyLmX15hOm5ZvZIMH21mdWkoKa5ZvZLM9toZhvM7PMjzHOpmbWb2brg58+SXVew3m1m9kawzoYRppuZ/U2wvV43s3NTUNOCIdthnZl1mNkXhs2Tku1lZvebWYuZrR/SVmZmz5jZW8Hv0lGW/UQwz1tm9okU1PWXZrYpeJ4eN7OSUZY97nOehLqWm9muIc/V1aMse9y/3STU9ciQmraZ2bpRlk3m9hoxG1L2GnP3Sf0DRIEmoBbIAV4DFg6b5w+Ae4PbNwGPpKCumcC5we1pwJsj1HUp8B8hbLNtQMVxpl8NPAUYcBGwOoTndA/x84xTvr2AS4BzgfVD2r4FfDm4/WXgmyMsVwZsDX6XBrdLk1zXlUBWcPubI9U1nuc8CXUtB740juf5uH+7ia5r2PRvA38WwvYaMRtS9RqbCnv6FwBb3H2ru/cC/x+4Ztg81wA/CG4/ClxmZpbMotx9t7uvDW4fAhqB2clcZwJdA/zQ414CSsxsZgrXfxnQ5O4nelHeSXH3lUDbsOahr6EfANeOsOgHgGfcvc3dDwDPAFclsy53/7m79wd3XwLmJGp9J1PXOI3nbzcpdQV//78NPJyo9Y3XcbIhJa+xqRD6s4GdQ+43c2y4Hp0n+ANpB8pTUh0QdCctBlaPMHmJmb1mZk+Z2VkpKsmBn5vZGjO7Y4Tp49mmyXQTo/8xhrG9AKrcfXdwew8w0nfZhb3dPkX8E9pIxnrOk+GzQbfT/aN0VYS5vZYBe939rVGmp2R7DcuGlLzGpkLopzUzKwIeA77g7h3DJq8l3oVxDvC3wL+lqKz3uvu5wAeBz5jZJSla75jMLAf4CPAvI0wOa3v9Go9/zk6r097M7CtAP/DgKLOk+jm/B6gDFgG7iXelpJObOf5eftK31/GyIZmvsakQ+ruAuUPuzwnaRpzHzLKAYmB/sgszs2ziT+qD7v6vw6e7e4e7Hw5u/xTINrOKZNfl7ruC3y3A48Q/Zg81nm2aLB8E1rr73uETwtpegb2DXVzB75YR5gllu5nZJ4HfAm4JwuIY43jOE8rd97r7gLvHgH8YZX1hba8s4HrgkdHmSfb2GiUbUvIamwqh/wow38xODfYSbwKeGDbPE8DgUe6PAb8Y7Y8jUYI+w38EGt39O6PMc8rgsQUzu4D485HUNyMzKzSzaYO3iR8IXD9stieA2yzuIqB9yMfOZBt1DyyM7TXE0NfQJ4CfjDDP08CVZlYadGdcGbQljZldBfwJ8BF37xxlnvE854mua+gxoOtGWd94/naT4XJgk7s3jzQx2dvrONmQmtdYMo5Op/qH+NkmbxI/E+ArQdufE/9DAMgj3l2wBXgZqE1BTe8l/vHsdWBd8HM1cCdwZzDPZ4ENxM9aeAlYmoK6aoP1vRase3B7Da3LgO8H2/MNoD5Fz2Mh8RAvHtKW8u1F/E1nN9BHvM/0duLHgP4TeAt4FigL5q0H7huy7KeC19kW4HdTUNcW4n28g6+xwbPUZgE/Pd5znuS6fhS8dl4nHmYzh9cV3D/mbzeZdQXtDwy+pobMm8rtNVo2pOQ1pityRUQyyFTo3hERkXFS6IuIZBCFvohIBlHoi4hkEIW+iEgGUeiLiGQQhb6ISAZR6IuIZJD/Bs/oYB1kpxocAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Question 5\n", "\n", "energie = []\n", "for i, t in enumerate(t_eval):\n", " r_x, r_y, v_x, v_y = split(res.y[:,i])\n", " energie.append(energie_totale(r_x, r_y, v_x, v_y))\n", "energie = np.array(energie)\n", "plt.plot(t_eval, energie)" ] }, { "cell_type": "code", "execution_count": null, "id": "f8ccd87a", "metadata": {}, "outputs": [], "source": [ "* code pour afficher une animation. res est le tableau provenant de solve_ivp" ] }, { "cell_type": "code", "execution_count": 9, "id": "e4246b9f", "metadata": { "tags": [ "keep" ] }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYjklEQVR4nO3df5BdZZ3n8ffXJDCZSBmQyI8GFFcqLKwLga6MLs6UChLIuhAZR2LVzuCPqoyubA0zu8zCYiGrf4jDOMzO4EhllBqdYRAdJWY1GIK4xTi1IB0Iv2EIDEoahCgkqGQkCd/9456ONzf3dj/d957btzvvV1VXn3vO0+d+c/rmfvo8z3POjcxEkqQSr5ruAiRJM4ehIUkqZmhIkooZGpKkYoaGJKmYoSFJKtaT0IiI6yLiuYh4oGndIRGxISIeq74f3OFnL6jaPBYRF/SiHklSPXp1pvE3wFkt6y4BvpuZxwHfrR7vJSIOAT4B/AawFPhEp3CRJE2/noRGZt4OPN+y+lzgS9Xyl4AVbX50GbAhM5/PzBeADewbPpKkATG3xn0flpnPVMs/Bg5r02YIeKrp8ZZq3T4iYhWwCmDBggWnHn/88T0sVZJmv40bN/4kMxd1s486Q2OPzMyI6Op+JZm5GlgNMDw8nCMjIz2pTZL2FxHxw273UefsqWcj4giA6vtzbdqMAkc3PT6qWidJGkB1hsZaYGw21AXAN9u0WQ+cGREHVwPgZ1brJEkDqFdTbm8A/h+wOCK2RMSHgSuBd0XEY8AZ1WMiYjgivgCQmc8DnwLuqr4+Wa2TJA2gmIm3RndMQ5ImLyI2ZuZwN/vwinBJUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVKzW0IiIxRGxqenrxYi4qKXN2yNie1Oby+usSZI0dXPr3HlmPgqcDBARc4BR4KY2Tf8xM99dZy2SpO71s3vqdODxzPxhH59TktRD/QyNlcANHba9NSLujYibI+LEPtYkSZqEvoRGRBwAnAN8rc3mu4HXZ+ZJwF8CazrsY1VEjETEyNatW2urVZLUWb/ONM4G7s7MZ1s3ZOaLmfnzankdMC8iDm3TbnVmDmfm8KJFi+qvWJK0j36Fxvvp0DUVEYdHRFTLS6uaftqnuiRJk1Dr7CmAiFgAvAv4/aZ1HwHIzGuB9wIfjYhdwA5gZWZm3XVJkiav9tDIzF8Ar21Zd23T8jXANXXXIUnqnleES5KKGRqSpGKGhiSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqZmhIkooZGpKkYoaGJKmYoSFJKmZoSJKKGRqSpGKGhiSpmKEhSSpmaEiSihkakqRihoYkqVjtoRERT0bE/RGxKSJG2myPiPiLiNgcEfdFxCl11yRJmpq5fXqed2TmTzpsOxs4rvr6DeDz1XdJ0oAZhO6pc4EvZ8MdwMKIOGK6i5Ik7asfoZHALRGxMSJWtdk+BDzV9HhLtW4vEbEqIkYiYmTr1q01lSpJGk8/QuNtmXkKjW6oj0XEb01lJ5m5OjOHM3N40aJFva1QklSk9tDIzNHq+3PATcDSliajwNFNj4+q1kmSBkytoRERCyLioLFl4EzggZZma4Hfq2ZRvQXYnpnP1FmXJGlq6p49dRhwU0SMPdffZ+Z3IuIjAJl5LbAOWA5sBl4CPlhzTZKkKao1NDLzCeCkNuuvbVpO4GN11iFJ6o1BmHIrSZohDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVKx2kIjIo6OiO9FxEMR8WBE/EGbNm+PiO0Rsan6uryueiRJ3avzM8J3Af8tM++OiIOAjRGxITMfamn3j5n57hrrkCT1SG1nGpn5TGbeXS3/DHgYGKrr+SRJ9evLmEZEvAFYAtzZZvNbI+LeiLg5Ik4cZx+rImIkIka2bt1aV6mSpHHUHhoR8Wrg68BFmfliy+a7gddn5knAXwJrOu0nM1dn5nBmDi9atKi2eiVJndUaGhExj0ZgXJ+Z32jdnpkvZubPq+V1wLyIOLTOmiRJU1fn7KkAvgg8nJl/1qHN4VU7ImJpVc9P66pJktSdOmdPnQb8LnB/RGyq1v1P4BiAzLwWeC/w0YjYBewAVmZm1liTtI8194xy1fpHeXrbDo5cOJ+Lly1mxRLnbEjt1BYamfl9ICZocw1wTV01SBNZc88ol37jfnbs3A3A6LYdXPqN+wEMDqkNrwjXfu2q9Y/uCYwxO3bu5qr1j05TRdJgq7N7Shp4T2/bMan1Jezu0mzmmYb2a0cunD+p9RMZ6+4a3baD5FfdXWvuGe2iSmlwGBrar128bDHz583Za938eXO4eNniKe3P7i7NdnZPab821m3Uq+6kOrq7pEFiaGi/t2LJUM/GHF4zfx7bduxsu37QOPaiqTA0pB6KDpPMO62fLk411lQZGlIPbXtp37OMduun+6/88cZeDA2Nx9CQCpW80R+5cD6jbcYvmmdjjfdXPvRufGU8jr1oqgwNqY3WgHjH8Yv4+sbRCbtzLl62eK9AgH1nY3X6K/+KtQ/yy12v9KXLqCTcpHacciu1aHetxfV3/KhoKu2KJUN8+rw3M7RwPgEMLZzPp897815v+u3erAG27djZ9jkuunETp115W0+v9ej1VGPtPzzT0IxWx9hAuzOBTnfRbNedM9FsrDkR7J7kfTl7fdbR66nG2n8YGpqx6poBNJl+/QROu/K2Sb3hTjYwxvR6oLqXU421/7B7SjNWXVdfT7Zff3TbDv7wxk18fM39+2xbc88op115G8de8u09XUxDHfa/cP68fbqMWjlQrelmaGjGqmsG0MXLFjNvzuQurEjg+jt+tNe4Q7uxkYu/di/bXnp5n5+fP28OV5xzIp8+783MGeeiDgeqNd0MDc1Yvb7Z4JgVS4ZYcMDke24T9jrLaXcmtPOV5Bcv770ugN8+dWhPd9Er43RfOVCt6WZoaMbq1Qygdl1I29vcCqRE81lO6RlPAl/fOLrnLMWzCQ0yQ0MzVsn01ol0upX5wl9vf6+oiW4H0vyGP5k3/+axmIuXLe74kZfeLVfTrfbZUxFxFvC/gTnAFzLzypbtBwJfBk4Ffgqcn5lP1l2XZoduZwB1Gkz/1127mfeqYOcrv+oqmveqgICdu9t3H7We5Vy8bDEXf+3evfYxnrEzkxVLhrjoxk3jtpms6b5tiWaPWs80ImIO8DngbOAE4P0RcUJLsw8DL2Tmm4Crgc/UWZPUrNObcCYQjRlNY2cxr/61uR0DY07EnnGJvUxiPL35zKTTDKupdF35wVDqpbq7p5YCmzPzicx8GfgKcG5Lm3OBL1XL/wCcHjFo9wTVbDXem/BYQPzLlf+Rf7rknR1vRgiNay+axyWgcRbTKWRaNZ+lfHzN/W3DbKpXbPvBUOqlukNjCHiq6fGWal3bNpm5C9gOvLZ1RxGxKiJGImJk69atNZWr/U27wfRm23bs5A3VAPlEn4nR+kZc2pXUPBbz8TX383d3/GifK9AXHDBn0uM1E9XhNR+aihkzEJ6ZqzNzODOHFy1aNN3laJYYG0wf79oIaHTp/OLlXY1xjXE0vxGXdCUNLZzPP13yzj1hcMOdT7Vt1zpNdzLqmpqs/VPdoTEKHN30+KhqXds2ETEXeA2NAXFpXO2myk7FiiVDfPZ9J03YbufuZOcrWXzx3URnMa3dTWvuGR33FiNTHYfw5oTqpbpD4y7guIg4NiIOAFYCa1varAUuqJbfC9yWOcWb82i/0evB3RVLhpjgJGKP3ZnMe1Xsc9V46xtx65Tgg3993l4D683dTWP/nvFMdRyiF1OTpTG1TrnNzF0RcSGwnsaU2+sy88GI+CQwkplrgS8CfxsRm4HnaQSLNK46PnmucGYs0Liye+H8eSw4cO6401hLpwS3+/e0M5VxCKfbqpdqv04jM9cB61rWXd60/K/A79Rdh2aXOgZ3hzp8MFEn23fsZNMnzpzy8zUrrXuy4xB+Frh6bcYMhEvN6hjcbdf3P16PVS8Hkkv3NdlxiG6n2/Zq3Eizh6GhGamOwd12ff9Xn38yf37+ybUPJE80aD7mqvWPTuqNu5szMi8KVDt+CJNmpLo+eW68MYg6xwXG9vW//s+DvDDORYST7V7q5rPA6xg30sxnaGjGmm2fPDf271lzzyh/9NVNHQfmJ/PGffGyxXuNaUD5WZIXBaodQ0OzTq9nC011MHmqdaxYMsQfdrhh4ZjSN+5uzsi6OUvR7GVoaFapY7bQVLppuq2j0xt28/ZSUz0j6+YsRbOXA+GaVeq4Od9Uumm6rWO8gfF+vXF7UaDa8UxDs0od/fBT6abpdJZQeh1Ic7fS6LYdzIlgdzZuYdIcPnW/gc+2cSN1z9DQrFJHP/x43TSdxi3G3uRbTXRjxGbNb9hepKdBYfeUZpV+Xb/x6fPeDNDxOoZONx4c74aE4/EzMTQoPNPQrNLP6zdOu/K2jm/knW5J0ukT+Sbi9FcNCkNDs06/+uHHeyO/+vyTezrzyOmvGhR2T0lTNN79r3o988jPxNCg8ExDmqKJrmPo5RlPXd1u0mQZGtIU9fuN3OmvGgSGhjQJrVNs33G8n1ev/YuhIRVqd63E393xoz3bvXZC+wMHwqVCJR/J6rUTmu1qOdOIiKuA/wS8DDwOfDAzt7Vp9yTwM2A3sCszh+uoR+qF0msivHZCs1ldZxobgH+Xmf8e+Gfg0nHaviMzTzYwNOhKr4nw2gnNZrWERmbekpm7qod3AEfV8TxSP5V8JKvXTmi268eYxoeAmztsS+CWiNgYEavG20lErIqIkYgY2bp1a8+LlCbS7oK9//yWY7x1uPYrkVO8gVpE3Aoc3mbTZZn5zarNZcAwcF62eaKIGMrM0Yh4HY0urf+ambdP9NzDw8M5MjIypbolaX8VERu7HQqY8kB4Zp4x3vaI+ADwbuD0doFR7WO0+v5cRNwELAUmDA1J0vSopXsqIs4C/hg4JzNf6tBmQUQcNLYMnAk8UEc9kqTeqGtM4xrgIGBDRGyKiGsBIuLIiFhXtTkM+H5E3Av8APh2Zn6npnokST1Qy3UamfmmDuufBpZXy08AJ9Xx/JKkenhFuCSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqZmhIkooZGpKkYoaGJKmYoSFJKmZoSJKKGRqSpGKGhiSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqVltoRMQVETFafUb4pohY3qHdWRHxaERsjohL6qpHktS9Wj4jvMnVmfmnnTZGxBzgc8C7gC3AXRGxNjMfqrkuSdIUTHf31FJgc2Y+kZkvA18Bzp3mmiRJHdQdGhdGxH0RcV1EHNxm+xDwVNPjLdW6fUTEqogYiYiRrVu31lGrJGkCXYVGRNwaEQ+0+ToX+Dzwb4CTgWeAz3bzXJm5OjOHM3N40aJF3exKkjRFXY1pZOYZJe0i4q+Bb7XZNAoc3fT4qGqdJGkA1Tl76oimh+8BHmjT7C7guIg4NiIOAFYCa+uqSZLUnTpnT/1JRJwMJPAk8PsAEXEk8IXMXJ6ZuyLiQmA9MAe4LjMfrLEmSVIXaguNzPzdDuufBpY3PV4HrKurDklS70z3lFtJ0gxiaEiSihkakqRihoYkqZihIUkqZmhIkooZGpKkYoaGJKmYoSFJKmZoSJKKGRqSpGKGhiSpmKEhSSpmaEiSihkakqRihoYkqZihIUkqZmhIkorV8nGvEXEjsLh6uBDYlpknt2n3JPAzYDewKzOH66hHktQbtYRGZp4/thwRnwW2j9P8HZn5kzrqkCT1Vi2hMSYiAngf8M46n0eS1B91j2n8JvBsZj7WYXsCt0TExohYVXMtkqQuTflMIyJuBQ5vs+myzPxmtfx+4IZxdvO2zByNiNcBGyLikcy8vcPzrQJWARxzzDFTLVuS1IXIzHp2HDEXGAVOzcwtBe2vAH6emX86Udvh4eEcGRnpvkhJ2o9ExMZuJxzV2T11BvBIp8CIiAURcdDYMnAm8ECN9UiSulRnaKykpWsqIo6MiHXVw8OA70fEvcAPgG9n5ndqrEeS1KXaZk9l5gfarHsaWF4tPwGcVNfzS5J6zyvCJUnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTM0JAkFTM0JEnFDA1JUjFDQ5JUzNCQJBUzNCRJxQwNSVIxQ0OSVKyr0IiI34mIByPilYgYbtl2aURsjohHI2JZh58/NiLurNrdGBEHdFOPJKle3Z5pPACcB9zevDIiTgBWAicCZwF/FRFz2vz8Z4CrM/NNwAvAh7usR5JUo65CIzMfzsxH22w6F/hKZv4yM/8F2AwsbW4QEQG8E/iHatWXgBXd1CNJqtfcmvY7BNzR9HhLta7Za4FtmblrnDZ7RMQqYFX18JcR8UCPaq3TocBPpruICcyEGsE6e806e2um1Lm42x1MGBoRcStweJtNl2XmN7stoFRmrgZWVzWNZObwBD8y7WZCnTOhRrDOXrPO3ppJdXa7jwlDIzPPmMJ+R4Gjmx4fVa1r9lNgYUTMrc422rWRJA2QuqbcrgVWRsSBEXEscBzwg+YGmZnA94D3VqsuAPp25iJJmrxup9y+JyK2AG8Fvh0R6wEy80Hgq8BDwHeAj2Xm7upn1kXEkdUu/gfwRxGxmcYYxxcLn3p1N3X30UyocybUCNbZa9bZW/tNndH4g1+SpIl5RbgkqZihIUkqNrChMdNuUVI9x6bq68mI2NSh3ZMRcX/Vruvpb5MVEVdExGhTrcs7tDurOr6bI+KSaajzqoh4JCLui4ibImJhh3bTcjwnOj7VJJAbq+13RsQb+lVbUw1HR8T3IuKh6v/SH7Rp8/aI2N70eri833VWdYz7e4yGv6iO530RcUqf61vcdIw2RcSLEXFRS5tpO5YRcV1EPNd8/VpEHBIRGyLiser7wR1+9oKqzWMRccGET5aZA/kF/FsaF6L8X2C4af0JwL3AgcCxwOPAnDY//1VgZbV8LfDRPtb+WeDyDtueBA6dxuN6BfDfJ2gzpzqubwQOqI73CX2u80xgbrX8GeAzg3I8S44P8F+Aa6vllcCN0/C7PgI4pVo+CPjnNnW+HfhWv2ub7O8RWA7cDATwFuDOaax1DvBj4PWDciyB3wJOAR5oWvcnwCXV8iXt/g8BhwBPVN8PrpYPHu+5BvZMI2foLUqq534fcEM/nq8mS4HNmflEZr4MfIXGce+bzLwlf3W3gDtoXMczKEqOz7k0XnfQeB2eXr02+iYzn8nMu6vlnwEPM85dFwbcucCXs+EOGtd4HTFNtZwOPJ6ZP5ym599HZt4OPN+yuvk12Ok9cBmwITOfz8wXgA007hfY0cCGxjiGgKeaHnd9i5Ie+03g2cx8rMP2BG6JiI3VrVGmw4XVKf51HU5ZS45xP32Ixl+Z7UzH8Sw5PnvaVK/D7TRel9Oi6h5bAtzZZvNbI+LeiLg5Ik7sb2V7TPR7HKTX5Eo6/1E4CMdyzGGZ+Uy1/GPgsDZtJn1c67r3VJEYkFuUlCqs9/2Mf5bxtswcjYjXARsi4pHqr4S+1Al8HvgUjf+kn6LRlfahXj5/qZLjGRGXAbuA6zvspvbjOdNFxKuBrwMXZeaLLZvvptHN8vNqfGsNjYtx+21G/B6rsdFzgEvbbB6UY7mPzMyI6Mn1FdMaGjnDblEyUb0RMZfGreJPHWcfo9X35yLiJhpdHT39z1F6XCPir4FvtdlUcoy7VnA8PwC8Gzg9qw7YNvuo/Xi2UXJ8xtpsqV4Xr6HxuuyriJhHIzCuz8xvtG5vDpHMXBcRfxURh2ZmX2++V/B77MtrssDZwN2Z+WzrhkE5lk2ejYgjMvOZqivvuTZtRmmMxYw5isY4ckczsXtqkG9RcgbwSGZuabcxIhZExEFjyzQGe/t6t96WfuD3dHj+u4DjojED7QAap+Nr+1HfmIg4C/hj4JzMfKlDm+k6niXHZy2N1x00Xoe3dQq+ulRjKF8EHs7MP+vQ5vCxsZaIWErjPaGv4Vb4e1wL/F41i+otwPamrpd+6tiTMAjHskXza7DTe+B64MyIOLjqqj6zWtfZdIz0F84GeA+N/rVfAs8C65u2XUZj9sqjwNlN69cBR1bLb6QRJpuBrwEH9qHmvwE+0rLuSGBdU033Vl8P0uiG6fdx/VvgfuC+6kV1RGud1ePlNGbbPD5NdW6m0de6qfq6trXO6Tye7Y4P8EkaIQfwa9XrbnP1OnzjNBzDt9Hohryv6TguBz4y9joFLqyO3b00Jhz8h2mos+3vsaXOAD5XHe/7aZpR2cc6F9AIgdc0rRuIY0kjyJ4Bdlbvmx+mMYb2XeAx4FbgkKrtMPCFpp/9UPU63Qx8cKLn8jYikqRiM7F7SpI0TQwNSVIxQ0OSVMzQkCQVMzQkScUMDUlSMUNDklTs/wPqRZZA/ZAvggAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from IPython.display import HTML\n", "from matplotlib.animation import FuncAnimation\n", "\n", "fig, ax = plt.subplots()\n", "l, = ax.plot(split(res.y[:,0])[0], split(res.y[:,0])[1], 'o')\n", "ax.set_xlim(-10, 10)\n", "ax.set_ylim(-10, 10)\n", "\n", "def animate(i):\n", " r_x, r_y, v_x, v_y = split(res.y[:,i])\n", " l.set_data(r_x, r_y)\n", "\n", "ani = FuncAnimation(fig, animate, frames=len(res.t))\n", "\n", "\n", "HTML(ani.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "id": "9eae45ad", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 22, "id": "617246e3", "metadata": {}, "outputs": [], "source": [ "#Question 7\n", "\n", "alpha = 0.2\n", "def force_friction(r_x, r_y, v_x, v_y):\n", " return -alpha*v_x, 0*v_y\n", "\n", "def force_totale(r_x, r_y, v_x, v_y):\n", " piege = force_piege(r_x, r_y)\n", " coulomb = force_coulomb(r_x, r_y)\n", " friction = force_friction(r_x, r_y, v_x, v_y)\n", " return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]\n", "\n", "\n", "\n", "def f(t, y):\n", " r_x, r_y, v_x, v_y = split(y)\n", " f_x, f_y = force_totale(r_x, r_y, v_x, v_y)\n", " return join(v_x, v_y, f_x/m, f_y/m)" ] }, { "cell_type": "code", "execution_count": null, "id": "d7b4147f", "metadata": {}, "outputs": [], "source": [ "N = 20\n", "\n", "r_x_0 = np.random.normal(size=N)\n", "r_y_0 = np.random.normal(size=N)\n", "v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "\n", "t_max = 1500\n", "t_eval=np.linspace(0, t_max, 400)" ] }, { "cell_type": "code", "execution_count": 24, "id": "3920953b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 10 s, sys: 59 µs, total: 10 s\n", "Wall time: 10 s\n" ] } ], "source": [ "%%time \n", "res = solve_ivp(f, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0), \n", " t_eval=t_eval)" ] }, { "cell_type": "code", "execution_count": 12, "id": "cb10757b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABDdUlEQVR4nO2dd5hcxbG3fzVps3a1WmWttMpCAUUEQgJEEkJCwgYHMB8mmWRzbRzwFcEYg20wvva1MdmXYDAGG0wSCAQWUUQFQAkJ5ZzDKqw29/fHOWfmzAndPWl3dqfe59GjmdM93TVn9lRXV1dXkxACDMMwTPsn0NoCMAzDMC0DK3yGYZgcgRU+wzBMjsAKn2EYJkdghc8wDJMjhFpbABkVFRWiqqqqtcVgGIZpMyxatGiPEKKzV1lWKnwimgFgxoABA7Bw4cLWFodhGKbNQEQb/cqy0qUjhJgthLiqtLS0tUVhGIZpN2SlwmcYhmHST1YqfCKaQUQPV1dXt7YoDMMw7YasVPjs0mEYhkk/WanwGYZhmPTDCp9hGCZHyEqFzz58hmGY9JOVCj8dPvya+kY8/ekm1Dc2p1EyhmGYtktWKvx08Mbynbjx+aW445UVrS0KwzBMVtBuFX6zebDLkx/7bjpjGIbJKbJS4afDh3+0oSmNEjEMw7R9slLhp8OHf7Q+pvAbm9iPzzAMk5UKPx3U2iz82sZmNDUL7D9S34oSMQzDtC7tVuHbXTp1DU146L21GH3Hm1iz63ArSsUwDNN6tF+FXx9z49Q2NmPJZmM94O7XV7aWSAzDMK1Ku1X4tY3xFn4oSACAnQdrW0skhmGYViUrFX46onRqbYu2tQ3N2F9TH33NMAyTi2Slwk9LlI7Nh//xur3Ysv8ogHjLn2EYJpfISoWfDuwK//ZXVmDj3hoAwMGjDfjNqyuw93Bda4nGMAzTKrRfhV/fhEjI/fX21zTgr++vx9VPLgIAfLXzEH70zGc8ADAM0+5ptwq/tqEJZQVh3/KFG/ejtqEJ767ajZc+34ZLH1vQgtIxDMO0PKHWFiBTHG1oQllhGLsOxSz3rh3ysPNg7P2QX7yOiuIIAGDDniMtLiPDMExL0o4t/GaUFUTirgWJXPX2HDaid6xkawzDMO2VdqvwLQvfzvCe/lE/9Wa+nd2H6rBsKx+8wjBM+yMrFX664vA7FsYs/OeumYCTB3WOvn/6yhPi6jc0CTQ0NeMHTy3GOX+Zj+qahqT7ZhiGyUayUuGnEoff3Czw5McbcaiuMc7C71tRhPxwMPr+uKqOGNkrvv2jDU3YbUbr/HvxliSlZxiGyU6yUuGnQiBA+McnmwAARXmxNeni/BDyw7GvGwoGkGcbAABjVtCrYwEAY7MWwzBMe6LdKXwAOH9MTwDApn010Wt5oSDyQvEKPt+h8F/4bCu2Vxu5dmr5LFyGYdoZ7VLhXzi+N04aWIHvHN877rrdwgeAAsf7O19bGU2fXMcnZjEM085ol3H4RXkhPHnF8a7rTove+d5OnWnhCyFAHuGcDMMwbY12aeH7ke906ZjvAx76vK6xGXe+9iXO+OO7EEJACIHbXl6O38/lfPoMw7RN2qWF74fTpWO971gYwV7H8Yeb9h7BQ+8eBADsO1KPT9fvw+MfbgAA3HDWkMwLyzAMk2bavcIfWVmGTkVGPL5r0TZivC8rDLsU/hFbPv0Ne2uwbJuxJ6Bf56JMisswDJMx2r3Cf+kHE6OvXRZ+yFL4EQCxXDrlRRHssw0Am/YdwZE6YwBobOIUDAzDtE1yyofvjLu3Fm0LI/HXy4vic/Bs2FODI3WNAIA6PkCFYZg2SotZ+ETUD8DNAEqFEN9oqX7t+PnwndE6Iccq7p/nrY6+ruP4fIZh2igpWfhE9CgR7SKiZY7rU4loFRGtIaJZACCEWCeEuCKV/lIlEnTG4RuK3n5QyvljemHy4C7R95dM6BP3mTrzTNzFm/bjPyt2ZkpUhmGYtJOqS+dxAFPtF4goCOA+AGcDGArgQiIammI/aYGIcPUp/fDsNRMA2Cx7m1v+D98aiZL82MTnV+cOR/fS/Oj72sYmCCFw3v0f4ntPLER9YzOWb+PsmgzDZD8pKXwhxHsA9jkujwewxrTo6wE8A+Bc3TaJ6CoiWkhEC3fv3p2KeJ7cePYxOK6qHAAQ9ArAB5DnOBqx1HZylhBGZk2Lfy7YhOn3zMfc5TvSLivDMEw6ycSibU8Am23vtwDoSUSdiOhBAKOJ6Ea/DwshHhZCjBNCjOvcubNftbRgbaAVECjJC2HGyB4AYou71njQwXFUon3hdv0eI1/P7+euyqisDMMwqdJii7ZCiL0Armmp/nQgGBpdCGDpr86KXrcs/IA5IpQ6FH6NLUZ/6wFD4W8/cDSjsjIMw6RKJiz8rQAqbe97mde0SccBKHr9GP87TzeMKvxAvMK3wjWvf+bzaN21u434/YZmjs9nGCa7yYTCXwBgIBH1JaIIgAsAvJxIA6kcgJIIVvx9R0fcfZ4jx06RWc9S+B/ZcuWv3W1k16xvbEZjE4dsMgyTvaQalvk0gI8ADCaiLUR0hRCiEcB1AOYC+BLAv4QQyxNst0Us/FMGdcavZg7DLdOPibueZ8bnW4eeh8xwzg75bg+YfXZwpK4J877cidteTujrMgzDtAgp+fCFEBf6XJ8DYE4K7c4GMHvcuHFXJtuGDkSES06scl13+vBDwXjF78fh+kZc8beFAIDrThuAL7cfxKQBFZxemWGYrCCnUivoEnXpmD4da+dt0Ka4zxvTE0O6lcR9zkq/AAC3vLAMFz/yKV78PKHlC4ZhmIyRlQq/pVw6fsQsfON9MGC8tyx9ALjj3OGu6J3qow3R1+9+Zewh+GIzb8piGCY7yEqF31KLtn5YOXYsl87lE6swfUR3XHVyv2idwkgwGp9vpWZYvjWm3I+aRyRur+ZwTYZhsoOsVPitjdOlU1YYwX0XjUG3DrEUC0SEDvmGwu9cnAcAWLBxP4D4JG07Dta1iMwMwzAqslLht7ZLx0q54My84DxApUOBsebducRQ+G8s34GS/BBGVZZF6+zwsPAbmpqjh6UzDMO0FFmp8FvbpWO5aAZ36xB3Pc+RXtmy8CuKjfj8hiaBSQMqUJwXC37aebAOzc0CNfWNWLLlAADgjldW4Iw/votdB2sz9RUYhmFcZKXCb20qivPwxOXjce93Rsdddx6Cbvnw7fn0Jw2scB20sm7PETz83jp8/f4Psf9IPeZ9uQtA/DGKDMMwmSYrFX5ru3QA4ORBnaMWvIXTwrfSKNuzbo6qLHMNDIs37cdHa/eiqVlg2bZq1NQb4Zv1fJgKwzAtSFYq/NZ26fjhPEDFsuzt59wO6lqCgoh1klYAJfkhfLJuHz7bfAAAsHRrddSytyJ5GIZhWoKsVPjZSsCxihs23zfYcuiEg4GohV8QDmJUZRleWbItas0v21odfW359af877s4WNsAhmGYTNJi6ZHbC7+cMTR6gIqVaqGxWeC6UwegsrwAQMzyzwsFMaZ3R7y/eg8AYGyfjlhnZtcEgNqGJqzZdRhf7TyMtbsOY3Tvji35VRiGyTFY4SfIZRP7Rl+P7dMReaEArj65H47v1yl6PbZxCxjduwwAMLBLMfpVFEWVPwAcrW/GviP1AIxoHoZhmEySlQqfiGYAmDFgwIDWFkVKeVEEq359tuu6ZeELAKMrO4IIGN+3HIWRIHYdioVi1tQ3Yr+p8Hcf4hBNhmEyS1b68LN10VaXqMIXQGlhGA/9v7G47rQBKCuMwH5OSm1DE/bXGL77nQfrUF3TgJc+3wrhPJGFYRgmDWSlhd/WyXfE4U8Z1g2A+2zcow1N2F9jWPjbqo9i0u/ewqG6RnTrkI+jDU0Y2avMdTgLwzBMsmSlhd/WKYi6dOIt9TKPs3EtH/7SLdU4ZKZX/mzzAVz62AJc+NePW0BahmFyBVb4GcBatHV6ZsoKPSx8U+Gv2xOL3vnYPEJx5Y5DAICfPfsFXl+2I1PiMgyTI7DCzwD2RVs7ZQUx90yH/BBq65uwz3TpNJnO/WCAogofAL5+/wd4btEWXPP3Ra5+GpuacePzS7HeNlgwDMP4kZUKPxtSK6SCn4VvPzClKC+ED9fuRW1Dc/QwdQA4oV85ahtiG7k+23Qg+nrnwVrsOliLqlmv4pN1e/HZ5gN4+tNN+O9/L8nMF2EYpl2RlQq/vUTpOCm1uXQKwkGsNlMkD+0ey8o5vIf/d/5s0358uNaw/v/+ySYcNVM0OFM+MAzDeMGaIgPEFH68iV9iS5tcZ6ZXGNO7DJdPim3mGtQ1dk7u9WcMjPv8nsP10YXd4rwQakyFXxCJH2DeWL4D877cmdqXYBim3cEKPwPY4/DtBAKEiuII/uu0AdGjD39w6gD0Li8EYPj1u5fFTtXqUVYQ9/m9h+txuNZS+MFo1s1Ch8K/6slFuOJvC9P3hRiGaRdwHH4GsFwsXtunFt5yJgDgL2+tAQD061wcPTClY1EEXW3HKHZyxODvPVIXrVsQDkYPTS/wcSExDMPYYYWfASyLe/Lgzsq6lR0LECBCMEDoWBiv8J2umr2H61FnLujWNTZHd+nmhfwnaos27kNtQzMmDqhI+HswDNO+YIWfAYryQnj/56fGKW8nZw7tinlf7oxm3OxcnIfyokjc8YiFkdjrcJCwYe8RhM36H63biy37DbdQQ7N/KobzH/gIALDhrunJfyGGYdoFrPAzRKXpl/fjr98dF5cz5+pT+qF7abzP3u6br+xYiOXbDkbfL9kSC1mta/A+OYtz8jAMYycrFX5byZaZKkSxA1XsaZd/eNoAdCgIx/nmy4sicbtx7WyvPoolWw7g2F5lcddrfQYChmFyk6yM0mnrcfip8pMpg/G9k/rFWfjWjtzzx/RCxOGz/3DtXsy89wPcOedL7Dkcy6v/2ab90ddrdh1ii59hcpysVPiMgd2Hf/vM4fjaqB646/wRGOWw5C0eem8dHv9gQ/T9d/7vk+jrM/74Ht79ajdqbefobj1wNJrLx4vqow14f/Xu5L8AwzBZRVa6dBgDK0UDABzfrxyTBhqRNlb0zvi+5dhzuA7rdh/BoK7F2FFdi5U7Dnq2BQCXPrYAAPCnb4/C9f/8HICxGHzL9KG4+IQ+rjN7r35yIT5etw9f/HJKXFoIhmHaJmzhZzF2H3/Ylj7B8u2P7l2GnubmrJL8MAZ3K8GqnYeU7VrKHgAamgR++fJyXPLYp3hlyba4el9uN9pqbGrG3sN1ePLjjXHl9Y3NmPGX+ZhvO7aRYZjshRV+GyRoWuId8sPICxnKvzASRO/yImzedzSpNt9fvQfX/eOzuGvNps9/w94juOrJRfjFi8uw2jagbK8+iqVbqzHreU7exjBtAXbptEHqm4zomw4F4eimq+K8EPp0koeCJoq1xmvF8gOxHEAAQKDotd2H6tC5JC+t/TMMk17Ywm+DNFoKPz8UVfiFkfQo/CVbDuB/3/wKQMzC98M60Wv3oToc95v/+NbbeuAoqma9ijlLt/vW2XWoFrsP1fmWMwyTOqzw2yANTYai7VAQRl7YsvCDrmRryTDz3g/w53mrUdvQ5Er+BsRmF3Y5LCbe9Rb+vWhL3LUd1bWYs8RQ9M4yO+N/M086aDQ0NWPJlgMa34BhGD9Y4bdBGkylmxcKRBO1FeaF0Lk45lJxpla2CAfJ87qTIb94HUdtIZwWtWZK5m0HjuIj28lcgGHJ3/DcFwCME7we+2A9TrhzHn4z50sAAHl0vfNgLS5+5BN3gY0bn1+CE+96CzPv/cAzCmlHdS3mLt+Bj9ftjduHYEcIgaZmgbpG93eyc6DGP0wVAA7XNUbPIWCYtgYr/Cznw1mnYf5/nxp3zVL44WAgGslTnBeK86GfPby7Z3t9OhWlJI81CEz63Vv4xYvLXOVWWp8XP9uKX81eEVe2dGs1Xl0Sc+vsPFiL37z6Jd5XRPk8/enmqLtn10G3Qv/Ggx/i6icX4YKHP8b5D3zoKt99qA59b5yD/jfNweBbXvfsY9nWajz+wXqMuv1NvPzFNs86ADD8l3Nxwp3zfMuf+XQTqma9Gs1k6qS5WeDiRz7BO6t2+bbx9qpdeGT+et/ymvpG3PzCUhw2z0bw4pH56/HGcv9zkFfuOIhfv7LCdzOeEAJ/eGMVNkiOz/xo7V6pm+5ofRPufWt11AXpxZyl27F8m//Jdpv31eCFz/xnhoBxz6trvO83YPy2izbu9y2vbWjCS59vlW5MXLBhH3ZU1/qW7zlcF7fR0YtP1++T3osd1caJdpmkxRZtiagIwP0A6gG8I4R4qqX6bst4uWkuntAHizcdQP/OxWhsNv6ACiNBFNkSr/mlTK7sWIA15klbyVDb0Iz/rNgJSb42APC0tHcerMMP/rEY0481Erkd/1t/xemHV7dWEjkA2Li3xpVmYuPeeKW1aW8NupXmx+1YPucv86OvP1i9BzNH9oj7zNsrd+Gyx419DH7KfOa986M5jrYdOOrau7Bgwz68uWIn3l+9Bws27MPKO86OK29uFrj/nTX4nzeMNZQrbAfjWCzfVo2H31uHlz7fhvKiCH46ZXBceVOzwJb9NbjjFWOw9Uqa19DUjIv++gn2HqnHNZP7o6LYvdh+oKYBf3lrDcoKI55yAMCFf/0YALDut9NcezgA4J63VuOBd9aiKC+EC47r7cr+CgDff2oxAOD9n5/qmX/q6/d/iD2H6zCiZxm6lebHJRcEDGU+6/ml+OfCzbjj3OEY3tO9O9/6bf/4rZE459gerp3qv5+7Co/MX48V2w5i6vBuGN27o6uNbz5oBC78fOpgXHtK/7iQaQCY+Zf52FZdi+vPGIhpI7rHHWQEAIs27se3HvoIQ7qVYOaoHvj+ZHfaGMuQGNO7DA/8v7HS5IvJkpKFT0SPEtEuIlrmuD6ViFYR0RoimmVePg/Ac0KIKwHMTKXfXOfro3thw13TUV4UQaPpR3f+Eds3bdnp17k4pb7fX70b33tCfbiKlztIh0m/eyuqoHcdrMVjH8Rbupc8+il++q8vpG3MvPeDONeMc3A6+fdve85OLP65cDP++7n4UNNnF22Oe//84i1odjRsT2j3ybq9roHmmw9+hIffWwfAGDgfeGdtXPk7X+2KKnsAuPKJhfhye7wLa/o98/HS58YM5P531uLu11fGlf/5P1/hlN+/E31/+h/ece2mHnjza9hrXjvv/g9ds4lfzV6O0Xe8CQD43esrUTXr1bjyXQdrMdYsB4B+N81xWeGPf7Ae/1yw2WxvBY791dy48h3VtXjo3dj3P+nut7F2d7whMu/LnVHD4Yw/vuty/e05XIeVO4ww4c82HYgbtO2yWvzkX1/gnnmrXXWsw4geem8dvn6/e4Zo5+7XV3kaTNtM6/9P/1mN8z3asORYueMQ7n59lbSPxZsO4PEPN0jrJEuqLp3HAUy1XyCiIID7AJwNYCiAC4loKIBeAKynhp2gacJaOA0HHAo/EkTItLo+vfl0rPvtNCy9bQoqO6a2sGudqatClrjNqSjtbNl/FH83N3hd+9Ril1sIAP69OKZcrIfVyRGbn90r2uidr/xdKoCh9O00OWT+yb++wPOfbfX9/G2zV8QpXi9+9/rKuCn+0fr4e/bmip3SA+qbmgXudwwazt9n7e4jeNvmPnK6LTbti80GLB6zpeeob3T/jnNX7IwOGBbWIGRx2+wV2Ger41zgv/rJhbjztfjBaqttpnakrtF1attnmw7EvZ9w5zz87Fn/wf/tlbsw3jGL3OFwmRjRY/6urxc/24oxtsENiP8uTc0CP7ZtZATiQ5cB0/W1zL+PhqZmvPR5/N+SKkIuWVJy6Qgh3iOiKsfl8QDWCCHWAQARPQPgXABbYCj9zyEZaIjoKgBXAUDv3r1TES8naDJdOiHHYmx+KIhQkNDYLFAQDiIQIJTkh1GSn1qKhE37arTq1Uos/JqGJtfU3I71QO06JPdn1jY0YcKdb3mWWYPK0i3VeN3jYbMUeGNTM55ZsNlVbmffkXpPf7k9jDTZxHQNTQLm3jk0ebThHGhUeLXRaGvDqXh1EUJE3RheA3aich7wcIvZ23AqTS9U32Wxh09dZmx4ccuLy1y/vV0Zr99zBC84Bn7nb2C5vvx48J21+MObX8VdS1ROXTKxaNsTMUseMBR9TwDPAzifiB4AMNvvw0KIh4UQ44QQ4zp3Vp8YletYh58EHT7UcJDw2KXjMX1EdxTZkrAV53srWq8ImmSZ8r/vuixAO7P+vQTrdvuvI6zZdRirdhxCQ6P8j162aGktbM+4d77n9NhSgk98tBG3SNw7ADDmjjfxwRr3zMYabLfsr8HfFFPwIz6y2i3odChSVRsNkkVDGXbl2ujRR6IWqdf3srdhrU2lgpecXgOiE+usaD851u4+HJ29eJU3NQu8JlnMBoBfv7ICW/YbxtM2j1lqkj+TkhaL0hFCHBFCXCaEuFa1YEtEM4jo4epq/9V7xsT8+7XCM9/48cm442vDQUSY0L8T7rtoTNyCmuXrH9y1BGcN6woAGNKtBI9cMi5tIn218zBmSyJdXlmyHZc89qlv+fw1e3DWn95TPvQ1df6ziHrFE9NkKrCdKURFWArlov/7BLd5uJ4sDtTUY9gv53qW1TUZ3+G1pdvxxEcb3HKafdQ1NkVdXX5s3lcTPfbSqw0geYVv/y2SHZjsA7RX/bW7D0d93X7t/UsxG7vlxaXR0F0vOV/6fBu+/dBHrut2ht46F4s27jPbcJf/6JnPcer/vAMA0TU0J9c+tVh6r/9v/nr819NGKhOvmUqmXDqZUPhbAVTa3vcyr2mT6/nwE+HWGUNx4fhKnHZMFwDAoK4luPiEPr71LdXfpUMexvYxohHyQgEUhFs2y4ZOzp8DklC7qlmv4gvJRiy/U8AsGqOKNHlTylJK2yXheoDbb2zHetivfWoxFjt81PY+7pm32ncmcrDWuE8n3f22p8ttza7D2LyvBne9thKjbn/TVQ4AD7+3FnWNTVEl5OTnzy3Bwg37fEMcP163D1P/9B6EENh6wPu3Hf7LudH1BC/r+7dzVkZ97l5rBwDw838vkYYu/v3jTbjk0U99+wCAT9bv8/28s06Dj9FhRWrJFHONYr+GZbCoZjvpJBNP+QIAA4moLwxFfwGA72SgHwZA1w75uPO8Y7XrlxVGABhHMFqzgkgoEHfYyncn9METH8mtyZbA74G1ePDdtb5l81buQnlRxLf8aEMTtlcflW7EenbhZkwd3s23/O8fb0T30gJf5QQA33jgQ/zkzEG+5bfPXo7LJ3qHPQLAuj1HcOUTC6VrHsfe9gb+fMEo3/LHP9ygjPr47ZyVEAK+M7NXlmzHK0vkboqVOw7hsQ824PZX/Gc7767ajRue/QJ7Dvu7/BZs2Id/fLLJt/xIfRPmSvYY1DY042h9U5xrxgvZuotlMKj0rmwdoaa+UZpW3JoBeM0EEnXl6ZKSwieipwFMBlBBRFsA/FII8QgRXQdgLoAggEeFEMsTbDcnjjhsDUZVluGBi8Zg8uAueNGMDIiEAijKMxR+gIDbzx2O4ryQKwLE4rwxPfH84oQmbRnBfsavk3vmrfYMwbMz4c63XPH2dm54bklcpImT/TUNuOmFpdI+Fm7cL1VOc5fvxNzlO6VtvLliJ6YM7Sqt85wkbYUuqpmKDm+ukH+Xg7UNUmUPxGLe/Xhj+Q5XhI+dusYmHHOr9wY7i2Vbq3HrS/5rN3WNza504E4O1NRjrWRPyxGJy9HqA2hDFr4Q4kKf63MAzEmh3dkAZo8bN+7KZNtg/Dl7hLEL17Lw80JBFJgLuyHzWshjI43Fd8b3xvQR3V1hcxb9Kop8z9/NNlbtkJ8fsHRr6utIupFNMlTK2B7SmCyqqCgd/MJkLbakQc4lit9E5yznN1bs9HShWew7UiedQQLA6X+QBycs3rQfZ/zxXd/y+qZmPDp/PV6TRJGlm6xMrcCLti2DtYAbCQZQaO7MDZuKPhjw/9PoWBSJ5uH3YlTvsvQJmWFUB8as2O4/i9Bl497UFf6GvfIBNB2KdKdH2opE8fPfW2xOw+CXjjZWSGaHgNoQACBV9gCUkVu1DU2+7q8M6fvsVPi8aNsyWAo/LxyIpmW48uR+ANxx/XY6FkZcO3vtXH1y/7Tn5m8t1qdhprIxDQrqUK3cH62KStIhHYpUFRufDrdROmZMzh3MTnROjlMhczkC8t80U3H4fABKDpNns/AjoQA23DU9upDljOu3U2o7eMXJ8l+dhaK8EN694VScdPdbSZ/AlS2kw5Waqel5utnVRs4jkEVv6aKaiei4hTJJWwrLZNoIdgvfwtpNOaxHB9/PBQPka+Hbr0eCyf958aHpTC6T5IZoJVmp8NmH3zLELHy3P/6kgZ3x7g2TlZ91Yl/sDaeg8Hum4TAXhmmrtKXUCinDPvyWwVL0eT6ZNZ2586eN6IY1vzFS+vpZ+Pa0sX6DwkMXj1XKpjofN5XZA8NkO+zSYdKOpeh1lec3x1ZGwzZlUToWXhb+tZP746xh3fDYpcehoth/Y5RXjnY7fjmBEqGyXD2LkCxlMEzG4LBMJu1E4/B9LHw7G+6ajlOHdIl9VhKlI6tjDQKnDumCX5wz1Pez3Uvlhz8ENLK9De3uvw4BAMd0k5cDQK+OqUUb+Z1LkAglaRjcdJDtvWBalpyy8Nml0zLY4/ATxc9dY8fLwg/bE7lJ+v32cZW46Hj/9Ng6Iv/5glG49MQq3/K+ndXHParCS7t2kM9EdNYiVPf/1MFdpOUAcP9FY6Tlv/7acFx9Sj/f8lGVZVjz22k4RjJIvvrDSbj93GG+5d1L8/GvqycoZVXtGtZBNQNU2QM6v4tsBpou/FJm5JSFz7QMpQVhFEWC0j/+n00ZhO9OcCdj0xkkPBW+PYpHMmiUF0Xwm6+P8C0PEOHzW8/E1GH+uW7KCiM4ZbB/iu2iSMjzCEA7KuX07XGV0vJQIIARHsfu2Xn6quOl36OiOE8p57QR3T1/J4uS/BC+Nqqnb3ls8PfXlPnhoPR3nz6iO46rch8PaGf9ndPQSaJIJw7ohMcuO07axqCuxVI5AWCk7YhLLwZ2VZ/8lurpcDr4rVXl1MYrpmUoygvh45tOx1kSZXPdaQNx+7nDXde9zjB14jULCHmkanby4g8mxp3P60WACGWFEZRLlEdhRK6gZHsNLC6eUIUfn+Gf/GxsVTk+uel033IiYPZ/TcLJg/wHni4l+ThnpPeh84Ceyw2QD8J5oYA0asr6rWS/q7Vfw7c8FHCd9eqEiKTrPxef0Ac9Sv0NkJMHdcYbPz5F2s9PzxyE74yXH5502hD1rKlfhXwGGFYMOoDaTeaX4C+nXDpMy1GSH9ZS3sng9UDoxOmPqixTtm0p67BE9oJwUKrkdBQ+AIRDEqtXoUittQaZRZoX1lPGKsIKZSxrx1KgsjuSF5IrfN0wXJkc4aBcTus+ysaVb4+vRL7HgekWPz5jEL47oUoq498uH48Oir0gN087RloOqNea/BR+Trl0eNG2beHnbvBSDiFbjh6ZglJh6WqZpRcIkDRFhO4ipcxyLogEpZae9XVDktxEeSH5TEQnIgqQK9xIMChV1tY3kN1PY9Dwl0VnIV9VLxwMSGc01neU6dFIMIB8SR+yvwnAWGg/ZVBn6d/Hk1eMx9Th/rOyvhVFmH3dJEh+dgDAZB+XY05Z+Lxo2zZxRpN4+vBtD1sqsfTnSvzRdmR96ET6AHJFqppFkKlKZQ9+flhuOetG+sgsY9UswroVsjsSVrh0dGciagvff1CxvoPstwsFA8gLSwYmzb87WR/hYAAFkj5G9CzFiF6l0jZW3H4W+vq4jdJwwqMnnEuHSQsrbj/L9cctC8sE9BWEk9nXTYqmflDp7LS4dCRt5CsUfnQmIlGlkaDKpaNr4Sfvf49Z+P7tR0IBqbLUdenI1wFIOSAAKmVMKVn40UFaUi0cJOlMxOojKJEzGCDfQUPn7N1kYIXPJE0wQBjZy5iFFUbcf0peysGuFHRdAE46FsXWHWSK1OjPv7ybItZfp438cFA+cGiYzkQk7UPbhy8bNMIBaR8xH77cBaZatNVB6hZSuJ4sN4tsYAoH5Ba+7sAkc2+FAvK1hnBAvQgeDgQ8nxsgS0+8YnIbK82CH54+fLtLJ0mFbx9IvJ7Jfp2L8POzBgPwf7jv/c5o7XhwmZwFksVBwG7hy5EpIZVFaiFVxkG5dU6uFx51SM/6ViH7PuEQSX3nOq6nQICkbjBVdI3Vh2wgDwVJPiAE1bMEmZxZeeIVk9uoQvC8FIBd6ej60GXterVwYv9O0QU1PyV41rBuMauWvNMg/+XC0a7+nMhcB0DsO6ruldTdkoa1hrxwUNqOjiIFVO4YTctZ0ks4KA/ttIpUfzs66wD+8ll9+NdRtWGVq9yGfuU5tWjLUTrtA5WFn6/pm3Zij+5J1oevitC5adoQzDDPu5Vb3yqFb/yfioWvOyzKLHj1QqXaVaLuQ09S2a1Px4IqIF/oVv1mFnKXjp4rUZXm22/HcBrOs/EkKxU+R+m0D/qUF7rijO2KrbQwjBe+f6K0jWevmYCBXeJ3PCaS88Vv+m5/mL1q2I949HJBlOSrd+na+1EPTJLQTk0LPxXrO2bhJz8T0bbwU7CcdQcm+TqByqVjuWPkMxEZ1qDy+GXjpftK8sNBz78jkUsWPtM+mDq8GxbdckbcNecUdnTvjtI2jqsqj1raFnEunSTcSk682ojbEZxC+KgqgZtOH7qeL9lAqFr4jboxFF81HRuv5C4dvS+r+t2lFr7iS1otzxzVw7eO9V0/vfl0jPY4w9naEFhZXohrTukv7c+LnNp4xbQPiOQLW7o4Bwn7+1RcJTLs0RXJtnHL9GNwk7kbMxXfuO6ERmaRqgat6FqDQtJ0LNpKLXztmYgcqQ/f7GPu9Sfj2F7+XoSeZQW+aTOsWV+Xknx0LPTYLWv7kslkIc1UWCYrfCbjPHzx2KgvMxn1Lw97TOGzkibiT+5KbtA6oV+nqCL3Gvgm9OuEt3822ezD/1Ec26dc2k9ZoXlvFREhMnTHZfnAlLp1rhqYrE/eco48rUEwQK7ZpYVlfQ/uVoIzjpFHavl9p7BtlpDswqyMDOl7VvhM5pkyrBv6a6Qi9kO2eWXmSPe0W2Wluup7VLc/pF4uAJ0eVAqwb+ei6E5LP4W/7rfTpKd/nTWsK1770Ula/cmIWs4pLNrqdl8lSTlt3YdHLx0XHci8+jixfwX+dvl4aT+divM8la19FuGpjG2X/JS1fV3H6+/TfiUZhc8uHaZdkIxOkj0ww3qUYsltU1KQyKdP+6JuknrUrhQ8F4ZtDXvNIojUlvmkgZ3R3cou6VG1W4f86IAgo1+FXipgacik+b8suygAjKsqx7PXeOfNt37r04Z0xYn9O3n0IVe0Othnb6pBUkvhZ8DCZ4XPtGmsXa35kh2QfigfqDSsEzjR3ewkI06ZKGYRXopU53up6nQqjkgPNQGMMw+uP2MgAMNyVvHX747ztPQtUX50xkDcNG2ItI0h3UqU/aiUsWqB2Q/7bMprwqJjnatcOnbRk1H4HIfPtGnuOv9Y/PFbIzGshzzUtqPHNN6ycqcM7erpl03mgbLj5QKKV8bJtatarFPJrZO22q6wvGo7labXQSuTB3eJhhFec0o//OGbI+PKT+hXjvduODX6/syhXT3PA7bHlCdrOavqxCnSNGzc85LTPvj69RGImyXI+2OFr4Dj8NsfHfLDOG9ML2mdW88Zihe+P9F13XroivNC6OSxUcX50KbD4Lc/6F4RHzrRR/YHvW8n9xpGojMXrxBP+3cfV1XuWitxDhoPXjwW3xwb/zvETUSIXMc6DuhSjN6Oa857/saPT0Yf23f0UqQXjq+UljtR+caTPcfBPntTDco6swivv4U411NSLp2EP6JFVip8Jjf5fyf0QZVHuljrofR7wFUP7co7pnou7lpUVbgXEe0P6YAuxbjzvBFx/Tg3xnw46zSc7jhFyd7GtZP74+dTB/uWA8B93xkTF/boLH/6qhNcZ/TaFWdxXgjzfjo5vg+PW+NUts5+nPfZS/EGHRbuoK4lvuWAMZj/+msjfMu98Pq97co12UVq++dUcqhi9nVIZiaSUxY+k5v4KW7rwVeV+5Eflmdg/Pv3jscPTo3fHOP04V84vrf04OweZQXo57Cu7fKGggHXYeRORTD92O4YYrPinV+rtCCMoT3irXwvhWW/5lWuUuhORep1f1VK0/mZ4rxQvFyOPo7tVeraceqlKGfYjoJM1ZXnJScQP+NJtgv24TOMAj/FbblcZYpdban5l3cpycfZjtOLvKxH1SKrUz6XYlWUG/361/eSy6uNoMIKdn4k0QHA+Ezstc6Cs7OK615p3M85Pzwpbk+C14BwzrH+p1BFZZHI6SzPxMZBHTg9MtNu+fe1E/Dmil2+5VZeG+lhEkRogv9DopoFOBWO11ReNbt3t+FUavH1PV0lcZazu0NnVImnwlIMGq6dywrl7NmGYmHTNaj41LH0mtfP4/yuztvhfP/dCX1w+7nD3Q3ZKIoE485ByNR5znGDSjIWPit8pr0ytk+5dDeppVBkD04gAKDJeF2c5/6zVvlRZekbtNtQWN86i8vxrhJ5uV+deP+6xixB5dNXWN9abiOfgae5Sfj2oZpp6MhpJxIMYPntU6V9AMDd3xjpuiZDpc6TsvA5Hz6Tq1hKTXogRSCAUEDgZ2cNdi1sqj7rVa6jxFTlbgtfZ6HSJkMSytr5OZX/3dkn4KVIPeRUbExTzRpibQhfOZ3uFNUMSXV7vcqdazX/882ROFNxMI4zyknVVzKLtuzSYXKWqEtHZuGTkefFLzNhonHhKgvfW4HJ+9SKTFEoax0545WNuw+XW0jlelK5dJKU0xnpo2rDNQBoDKikKHcPoG457Izt0xG//6Z6BnDuqJ7R18lY+I2805bJVaKLtjIffoAUB0Yr+lD4360+UmlDx9BTRtg4B5VkonScA5EqSkexaOu9eUnep7Nf1SK5l5xuuV1NKMsTdQvp6O6515+MyvJYqK/Xb3DlSX2lbQiRGT8+K3wm67EeQpnSDgZIEcUj/1MPOkxhHUXpKldYtVqZO5WLoXKl5+xHJ/rF7TrRUPgJKmtPt1CCA1Oi+wcAZ4ilug+1W0hnlhb/3jnwX3piFW6ePlTZTias/BZT+ETUj4geIaLnWqpPpn1gPdgypR0MyA+/TtTC13FTOHErj8SsR0OOxGTwduko3C2K7+p+7yGnyh2jNWh4y+zXbzKhnbLPG32kwcJ3XFO5nnTWGlbeMTXptNwytBQ+ET1KRLuIaJnj+lQiWkVEa4holqwNIcQ6IcQVqQjL5CZRha8Iy5Ra+LbPdvFIN5yORVvVoOIaEBR1krF6gcTj8JXRRQo5dO6VahagNWgo9wu425B93pDBUUep8BO38HUGPzuRYAD5ikPnk0V30fZxAPcCeMK6QERBAPcBOBPAFgALiOhlAEEAdzo+f7kQwj/QmmEkxMIy/esEAoSgZAZsPex/vmAUpgx1JxDTUvhKH3Fi1qGXuPFJudSKVLkYqqOMXUrP+V5+L3RcJarBS+e7utcv1HKqyp2zxlR/Y686TndhOvpIFi2FL4R4j4iqHJfHA1gjhFgHAET0DIBzhRB3AjgnWYGI6CoAVwFA7969k22GaUfounRkocvWoNGhIIyCiDsZmmpB0HlNZZ37yajCXkVn0PFM72tvw6NLpW880XUCT5ePv0yebWjE4aveqyxinVmEqg0dXeysk4k+kiUVH35PAJtt77eY1zwhok5E9CCA0UR0o189IcTDQohxQohxnTt3TkE8pr0QU/iSOkTSnbCWtejnFnJaYbI2AB/rXOmbVbt0VNa5jntA2YZilqDyQTv7TY+F7ypWLi7rLAzLZDL6cNaRt6G3aKu6v4l9Pp20WBy+EGIvgGt06hLRDAAzBgwYkFmhmDaB9cCowjKFxH6xFnT9FnbTcdiIelFX2YUjSsejDS3Xk70NtTJW72D1kEPRh95uXXm5So5EF22T2c3rbsN9zXmegsvCT3RxWVqaGqlY+FsBVNre9zKvpQznw2fsxMIy5QpfplCjg4ZPJaelp3JBeLp0lIu6iQ0qyS7aJu4WUrhOFLMEr3ulE+6o/q7yNnXklLXn9ZlE8yXp1HG5dBSfz1aXzgIAA4moLxFFAFwA4OV0CMUnXjF2rO3vsrDLAJFiJ6580NDJe55oArZEy406sdfeilbdZqIuHWcV0ugj0Rh6z5lGgoNGqq4SnQynyawDuPuR95Hofo50ohuW+TSAjwAMJqItRHSFEKIRwHUA5gL4EsC/hBDL0yEUW/iMHUtZqNIjq1w+gP/DZjVtFXvlvle5+VUuH71cOnJFqvJjO/vxdEFQvKJ1tZngLELLHaNYj9BrwylDoorUfS3RfEfJ+PDd5fLPZ9Kloxulc6HP9TkA5qRVIoZxoFpwBSwfvsTCjy78etchMlxC351QhdtmDvPtQ4Z6MU5ebtRRKEGXwpe34R1t5N+eV7+eG6+UylreplFHPjCpfN9pWbR1DRqqNuTlgIZLJg2GQbJkZWoFdukwdnQtfNXGLPv/XoQCAS23EOC3M1T+oOq4IBJNSpaqS8c7Tt/xXuEK8bKcVZkune1qZfVUtJmMS83lflEuqLaAhZ9BEz8rFT67dBg7vToW4MLxlTixfyffOqqdtiGFhQ8YMfqlBWH/Pmyf9TpQ3Cq/9MQqzP/vU13lWu4ARXpkVQQI4FgH8GxDYVknujCcZPRLolE6zirWrMz+XoaOIlXO0jQ0ZjIDjz2ldyZ22FpwemQm6wkFA7jzvGOldSKhAKjJv1zl0gGA5689EZ2KI8o2/vG94zGil9sYsR7UTkUR9OroPhjd6l+W61xl9eqFO+rPEvRcOqnNInzbSHDg8ZMjdoiKuw07yUTYONFRxkqPjse122YOw8Cuxbj5hWVabqNkyUoLn106TKLMOnsIbpp2jG957BAV/zZ6dypEkcdpWdE2zIf9+H6dUJLvngnouJ7sRX0rij3K06CMVYOGypWiseCa6MYrzygcVaSPzuCmkCO+rrRYs43kBo3HLzsuGgjg9/dhuYt03EbJkpUKn106TKIM71mKkZVlvuXDepRiVGUZupUWJN1HLNLHu9y6LlMKASJ065CPl6+biOkeB27HWd8eYUHp2MFqb8JvxqNy2SgHJsUhK87P6WTLVC3spho2a8ghL/fq4pIT+8S34aFVJw/ugn4VRdI+Yn8/KimTJysVPsOkm0FdS/DiDyZ6nnerC5HxMPpN6wOKAcEoM/YLHNurzLcPi3TsYE0mX767jrxctfjsL6e8D/csIXG3kOrzTms6GbfQsb3K4tZsfAcWhUFgXc6kD58VPsNoEiT55i6dw9aDAZJakXYF5rXRzHkpqWyZGm6QRCJo9Nwx7j7iIpKU7hpPMZXZRe0KXS+k0l3pxR9MjAYM+Opy+2zFp+3oXg9F3zkXpcM+fCYbMZS1xF2jcRSjMUuQtWGUXXxCH1zicRi71g5WpbvFrqx95FAqY3+ZvD6jmomowlz9BlHVTEQlQ0VJRFlnVGUZJg6oMOT0Udeq+2G/7mvhKz6fDrJS4bMPn8lGAopY/9gD7d9GMKCXAmL6sd091yRUueztbXjVB9SWtSWnrI9EZhFe711yKsp93WgadWJ13de6lOTj05tOj/rX/X4aq2n/9Ru5i8t+3d+Hn6MWPsNkI0qXjkboZ4DkLh3ro74WbYIhk6o4fT8lab+clP9dZ9E2gXUCv1saVMTy3zhtSPS133ft0iE/mq8paes8AQvffw1I/vl0wAqfYTTpVBxBeZEkTl/xQFt1VCkijHo+n9fZBUver2NtJOYqSWZzV3p2Fdv785bTnvTO66tcNrEvXv3hJN9yZ/vKCBofjamzVkAk/21jYZmZgxU+w2hyzSn98fz3T/Qt10njHFD48ElhSVrXy4si6NOpEGEPJ3wirhKdmYTKR++ZhTLBKB1VG37346dTBmF07zJpHZV1bu/f37+usM41Zk2qsF2VyycdZKXC50VbJhvJDwc9s2haRN0xCgtf57B133BJ8/rMkT3w7g2nJuXSibecveUg1aChWPhNeIOYhxwdbGku/OQ8b0wvHFdV7pLJSxb5/ghvuS0S8eH7W/Dxbbn7UMuZKlmp8HnRlmmLWEpN9ryqDmpR+fD1lFeszDu0U8elA2kdlfWt49IJKPoYXVmG2ddNQkVxnnSQVCtjd3/uOgpXmubMy5AnSR9+rlr4DNMW0cnXQ6Rw+SgGjdjCnlyO0oIwbjx7CKYM6+Yu13CVqCxWZS4dp4Xv0c9pQ7pG3TGem6KIMKJXKYIBvQHO3zrXt/CVytjn86Txu6h8+DqDeaqwwmeYNKHzwCpj+VUuHY1BJWgOKlef0h9dO+T7tiFrJ5GNV54bnjQWl78xthdumDLYLPcUI/pZaairUpH6yxCVV9VGdCD28/Hr9GHVlbeRSThbJsOkCb3kaSQ9OSuW5C05t4BVR25pGv9fNrEKJ/TzTjmtzqhpq+vRV3EkhOkjumP1rkP4audhyVnCei4q+fc1/le7UnybsLmFkpsl6PQRG3j8ZMhRC58XbZm2iMqXbJXpuSfk5bKsn4GA4rhHs2xi/wqc5eHysfq3lL5XJJAqH08gQLjvojEY26c82p53P5pusDRENqXiFlIpa92BWNZH9O8ng1o5KxU+L9oybZEOZsrkDh6pky0M6zt1xSNTgD3KCtCjzD8rqNYGsQDhpIEV+N35I3BM9xJX+eTBnXH5xL4oigQVJ5FZcit84ym4dFSDpO5ALJPDcsP4pjbW6EMVdhkdNHItPTLDtEV6dyrEnB+ehElm3hUvggFSRItY//v78EvyQnEhi06uP30gnrtmgm95WWHE/N+/jXAggKJICN8+rren1TqgSwlunTFUebSkyoLXc+moLOf4tpLpIzpL8HWlWfX8Pm/VU/ehXgfwbSJl2IfPMGlkaA/30Yd2SGHhq6xvIsKrPzwJXTpI9gMECAGJlTigSzHeu+FU9O7kfSoXANw6Y6j0uMe4vjTCHVVWrXTQCBACQnJKmEKhq3z89jqphmXKdLXSXWfdxwz68FnhM0wLEgyo/NXqxT+ZotZF1cZEySzFTiQY8PTxW6gUur5LR+7jt7flJ4OeWyg5l43Wbl5HXXcfajlThRU+w7QgqoiTiKk8I7JV2SziD98aid7l/oOHykWluwtWL0upSlmr20jWv643cMnrqAaEdMAKn2FakKtP7o+CiL8yP2t4N4SChC4e8fPZyEkDO0vLlRvJNK1vItnB79DrQzKGqgYF1YY31RqAvW1l+Ki/mCnDCp9hWhCvc2ztlBaEcd6YXi0kTeaxImxUKYFVbq7UIpvkilanDVXyNKNMYZ0r3EI6M5FUaRvzRoZh2iSlBWHp4q+OMk41tYLOTlv1piidNuQBlWq3kDVVkTSSIllp4RPRDAAzBgwY0NqiMAyTApeeWIVzJLMaPWWs59JR72D1l1PX3aKK5U8lEoggL08HWWnh88YrhmkfFESCqJQu6ppRPFL/uuJYSCss0zfzp0YEjSpkUsNHT1CktFC4hUgxA0gHWanwGYbJDfT867rKOh0RNMm3ofLhW+sV/ovL8fUyASt8hmFajaK8EAIEpZ9fJy1CKjHy6vQMauvbcD35y6k8VUtDzlTJSh8+wzC5QeeSPLzx41PQt6LIt47aN64bh+8vh+oAcdWgYpXpbbzy/3ymYYXPMEyrMqBLsbQ8QCTPDqobh6/hFvLdFKXZhmrHsKyNlgjLZIXPMExWEwoS/GN0EkiPLF1wja/rbsPqSyaHXuinn1dIJ61GqrDCZxgmq/nJmYPQLNH4qeay12tDbx1Anh459T5ShRU+wzBZjXWIih9B3bQHadl45S+Haq1B1UbOxuEzDMPoouOy0T1pTLnxSpECQseHr9xpm8E4/Baz8InoawCmA+gA4BEhxBst1TfDMO2XyYO64GdTBqFvJ1mkjzpkUmsnrkSOVKN0dGYRqaJl4RPRo0S0i4iWOa5PJaJVRLSGiGbJ2hBCvCiEuBLANQC+nbzIDMMwMUoLw7jutIHKTJUql45Ozn25y0YxaATks4hYH/5tpIquhf84gHsBPGFdIKIggPsAnAlgC4AFRPQygCCAOx2fv1wIsct8fYv5OYZhmBbBiKDxL9c/ZEU+aMg2ZqkUevTc3NZetBVCvEdEVY7L4wGsEUKsAwAiegbAuUKIOwGc42yDjGHtLgCvCSEWpyQ1wzBMAgzuVoK+Ff7x/oGA+tQt4395HanbSKHQVZu/0kEqPvyeADbb3m8BcLyk/n8BOANAKRENEEI86FWJiK4CcBUA9O7dOwXxGIZhDF6+bpK0PNWc+4Dho08lEqhdpUcWQtwD4B6Neg8DeBgAxo0bJ9tvwTAMkxYCSpeP8b9y4TeFA91Vm7/SQSphmVsBVNre9zKvpQwRzSCih6urq9PRHMMwjBS1D18jll9xUItq4TeWoC1zpKLwFwAYSER9iSgC4AIAL6dDKM6HzzBMSxIMkDKOH1D7+WXKWp0tM76vTKDl0iGipwFMBlBBRFsA/FII8QgRXQdgLozInEeFEMvTIRSfeMUwTEty0fG9MbZPR99yreRp8Lfejc/G/+8uz54onQt9rs8BMCetEhntzgYwe9y4cVemu22GYRgnA7qUYECXEt9yXQtf68Qrn3kAuV6kH06twDAMo0AvxbJelE4qaZxTJSsVPi/aMgyTTehtvFLk0lGcrdsCUZnZqfB50ZZhmGzCOghdfti6XrZM/0NWjP9zzsJnGIbJJgZ1KcFN04bgpIGdfevo+vDVB7UkL6eKrFT47NJhGCabCAQIV53cH0V5/nEu6t268f+7Px99lZyQGmSlwmeXDsMwbQ39nPvy8lZPj8wwDMPIMXz48nLj/9ZLj8wKn2EYJg0YUTrJp2fImvTILQ3vtGUYpq1xXFU5SvL9Vap6p631f44pfN5pyzBMW+MX5wyVllP0f3nytEzCLh2GYZgWIHbEoU85x+EzDMO0D6I+fB+fTjTXTq4t2nIcPsMw7Q1y/O8qjx5xmDkZslLhcxw+wzDtjUjIULdhn/wMsTj9HFu0ZRiGaW9MH9EdZQVhdC7J8yyPzgByzcJnGIZpbxTlhTBlWDff8pxNj8wwDJNr5Gx6ZIZhmFwjZ9Mjc5QOwzC5Rs6GZXKUDsMwuUbsCMQcs/AZhmFyjZZIj8xhmQzDMFlAIEC4edoxmDzY/1StVGGFzzAMkyVceXK/jLbPLh2GYZgcgRU+wzBMjsAKn2EYJkfISoXPcfgMwzDpJysVPsfhMwzDpJ+sVPgMwzBM+mGFzzAMkyOwwmcYhskRSAjR2jL4QkS7AWxM8uMVAPakUZx0k+3yASxjOsh2+YDslzHb5QOyS8Y+QgjP7bpZrfBTgYgWCiHGtbYcfmS7fADLmA6yXT4g+2XMdvmAtiEjwC4dhmGYnIEVPsMwTI7QnhX+w60tgIJslw9gGdNBtssHZL+M2S4f0DZkbL8+fIZhGCae9mzhMwzDMDZY4TMMw+QI7U7hE9FUIlpFRGuIaFYrylFJRG8T0QoiWk5EPzKvlxPRm0S02vy/o3mdiOgeU+4lRDSmheQMEtFnRPSK+b4vEX1iyvFPIoqY1/PM92vM8qoWkq+MiJ4jopVE9CURTcime0hEPzZ/32VE9DQR5bf2PSSiR4loFxEts11L+J4R0SVm/dVEdEkLyPh783deQkQvEFGZrexGU8ZVRHSW7XrGnncvGW1lPyUiQUQV5vtWuY8JI4RoN/8ABAGsBdAPQATAFwCGtpIs3QGMMV+XAPgKwFAAdwOYZV6fBeB35utpAF4DQABOAPBJC8n5EwD/APCK+f5fAC4wXz8I4Frz9fcBPGi+vgDAP1tIvr8B+J75OgKgLFvuIYCeANYDKLDdu0tb+x4COBnAGADLbNcSumcAygGsM//vaL7umGEZpwAIma9/Z5NxqPks5wHoaz7jwUw/714ymtcrAcyFsSm0ojXvY8LfqbU6zsiXASYAmGt7fyOAG1tbLlOWlwCcCWAVgO7mte4AVpmvHwJwoa1+tF4GZeoFYB6A0wC8Yv6x7rE9dNH7af6BTzBfh8x6lGH5Sk2FSo7rWXEPYSj8zebDHDLv4VnZcA8BVDmUaUL3DMCFAB6yXY+rlwkZHWVfB/CU+TruObbuY0s8714yAngOwEgAGxBT+K12HxP5195cOtYDaLHFvNaqmFP30QA+AdBVCLHdLNoBoKv5ujVk/xOAnwNoNt93AnBACNHoIUNUPrO82qyfSfoC2A3gMdPt9H9EVIQsuYdCiK0A/gfAJgDbYdyTRciue2iR6D1r7WfpchgWMySytLiMRHQugK1CiC8cRVkjo4z2pvCzDiIqBvBvANcLIQ7ay4Qx5LdKXCwRnQNglxBiUWv0r0kIxpT6ASHEaABHYLgjorTyPewI4FwYA1MPAEUApraGLInQmvdMByK6GUAjgKdaWxY7RFQI4CYAt7a2LMnS3hT+Vhj+NYte5rVWgYjCMJT9U0KI583LO4mou1neHcAu83pLyz4RwEwi2gDgGRhunT8DKCOikIcMUfnM8lIAezMoH2BYQ1uEEJ+Y75+DMQBkyz08A8B6IcRuIUQDgOdh3NdsuocWid6zVnmWiOhSAOcAuMgcmLJJxv4wBvcvzOemF4DFRNQti2SU0t4U/gIAA80oiQiMhbGXW0MQIiIAjwD4UgjxR1vRywCslfpLYPj2revfNVf7TwBQbZuCpx0hxI1CiF5CiCoY9+ktIcRFAN4G8A0f+Sy5v2HWz6iVKITYAWAzEQ02L50OYAWy5B7CcOWcQESF5u9tyZc199BGovdsLoApRNTRnMlMMa9lDCKaCsPFOFMIUeOQ/QIzyqkvgIEAPkULP+9CiKVCiC5CiCrzudkCIzBjB7LoPkpprcWDTP2DsVr+FYzV+5tbUY5JMKbNSwB8bv6bBsNnOw/AagD/AVBu1icA95lyLwUwrgVlnYxYlE4/GA/TGgDPAsgzr+eb79eY5f1aSLZRABaa9/FFGJEOWXMPAfwKwEoAywA8CSOSpFXvIYCnYawpNMBQSlckc89g+NHXmP8uawEZ18Dwd1vPy4O2+jebMq4CcLbtesaedy8ZHeUbEFu0bZX7mOg/Tq3AMAyTI7Q3lw7DMAzjAyt8hmGYHIEVPsMwTI7ACp9hGCZHYIXPMAyTI7DCZxiGyRFY4TMMw+QI/x9aQZCb/7fWBQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "energie = []\n", "for i, t in enumerate(t_eval):\n", " r_x, r_y, v_x, v_y = split(res.y[:,i])\n", " energie.append(energie_cinetique(r_x, r_y, v_x, v_y))\n", "energie = np.array(energie)\n", "plt.semilogy(t_eval, energie)" ] }, { "cell_type": "code", "execution_count": 13, "id": "63e56e79", "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWGUlEQVR4nO3dfbAldX3n8fcnA7izhHJARmAGVIwsLKwr4K1RV5JSQR6mXEGLJMMfCUarJrpSFfbBLCxVaOluqSHGVBYjNVEqJusqeRCckkEYlC1iakHu4PAMYWRxmQvCKAKysirsd/84fc3hcs69PffcPveB96vq1u3T/Tvd3+l75nxO9+/XfVJVSJLUxi8tdgGSpOXD0JAktWZoSJJaMzQkSa0ZGpKk1gwNSVJrCxIaSS5P8liSO/vmHZRke5L7m98HDnnuuU2b+5OcuxD1SJK6sVBHGn8OnD5j3gXAN6rqKOAbzePnSXIQ8GHgDcAG4MPDwkWStPgWJDSq6kbg8RmzzwS+0Ex/AThrwFNPA7ZX1eNV9SNgOy8MH0nSErFPh+s+pKoeaaa/DxwyoM164KG+x7ubeS+QZDOwGWD//fd//THHHLOApUrSyrdjx44fVNXaUdbRZWj8QlVVkpHuV1JVW4AtABMTEzU5ObkgtUnSi0WS7426ji5HTz2a5DCA5vdjA9pMAUf0PT68mSdJWoK6DI2twPRoqHOBrw5ocy1wapIDmw7wU5t5kqQlaKGG3H4J+J/A0Ul2J3kf8Ang7UnuB05pHpNkIsnnAKrqceBjwC3Nz0ebeZKkJSjL8dbo9mlI0t5LsqOqJkZZh1eES5JaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LUmqEhSWrN0JAktWZoSJJaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LUmqEhSWrN0JAktWZoSJJaMzQkSa11GhpJjk6ys+/nqSTnz2jzliRP9rW5uMuaJEnzt0+XK6+q+4DjAZKsAqaAKwc0/buqekeXtUiSRjfO01MnA9+tqu+NcZuSpAU0ztDYBHxpyLI3JbktyTVJjhtjTZKkvTCW0EiyH/BO4K8HLL4VeGVVvQ74r8BVQ9axOclkksk9e/Z0VqskabhxHWmcAdxaVY/OXFBVT1XV0830NmDfJAcPaLelqiaqamLt2rXdVyxJeoFxhcY5DDk1leTQJGmmNzQ1/XBMdUmS9kKno6cAkuwPvB343b557weoqsuAs4EPJHkWeAbYVFXVdV2SpL3XeWhU1f8BXjZj3mV905cCl3ZdhyRpdF4RLklqzdCQJLVmaEiSWjM0JEmtGRqSpNYMDUlSa4aGJKk1Q0OS1JqhIUlqzdCQJLVmaEiSWjM0JEmtGRqSpNYMDUlSa4aGJKk1Q0OS1JqhIUlqzdCQJLVmaEiSWus8NJI8mOSOJDuTTA5YniR/kmRXktuTnNh1TZKk+dlnTNt5a1X9YMiyM4Cjmp83AJ9tfkuSlpilcHrqTOAvqucmYE2Swxa7KEnSC40jNAq4LsmOJJsHLF8PPNT3eHcz73mSbE4ymWRyz549HZUqSZrNOELjpKo6kd5pqA8m+bX5rKSqtlTVRFVNrF27dmErlCS10nloVNVU8/sx4Epgw4wmU8ARfY8Pb+ZJkpaYTkMjyf5JDpieBk4F7pzRbCvw280oqjcCT1bVI13WJUman65HTx0CXJlkelv/vaq+nuT9AFV1GbAN2AjsAn4C/E7HNUmS5qnT0KiqB4DXDZh/Wd90AR/ssg5J0sJYCkNuJUnLhKEhSWrN0JAktWZoSJJaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LUmqEhSWrN0JAktWZoSJJaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LUmqEhSWqts9BIckSSG5LcneSuJL83oM1bkjyZZGfzc3FX9UiSRtfld4Q/C/z7qro1yQHAjiTbq+ruGe3+rqre0WEdkqQF0tmRRlU9UlW3NtM/Bu4B1ne1PUlS98bSp5HkVcAJwM0DFr8pyW1Jrkly3Czr2JxkMsnknj17uipVkjSLLk9PAZDkl4G/Bc6vqqdmLL4VeGVVPZ1kI3AVcNSg9VTVFmALwMTERHVXsTTcVd+Z4pJr7+PhJ55h3ZrVfOi0oznrBA+g9eLR6ZFGkn3pBcYXq+orM5dX1VNV9XQzvQ3YN8nBXdYkzddV35niwq/cwdQTz1DA1BPPcOFX7uCq70wtdmnS2HQ5eirA54F7quqPhrQ5tGlHkg1NPT/sqiZpFJdcex/P/Py558175ufPccm19y1SRdL4dXl66s3AbwF3JNnZzPtPwCsAquoy4GzgA0meBZ4BNlWVp560JD38xDN7NV9aiToLjar6FpA52lwKXNpVDdJCWrdmNVMDAmLdmtWLUE2PfSwaN68Il1r60GlHs3rfVc+bt3rfVXzotKMXpR77WLQYDA2ppbNOWM/H3/1a1q9ZTYD1a1bz8Xe/dtE+2dvHosXQ+ZBbaSU564T1S+b0j30sWgweaUjL1LC+lMXsY9HKZ2hIy9RS62PRi4Onp6Rlavo0maOnNE6GhrSMLaU+Fr04eHpKktSaRxrSLLx4Tno+Q0MaYvriuelrIaYvngNGCg6DSMuZoSH16X9D/6WE52bcCm364rn5vsl3FUTSuNinITVm3pZjZmBMG+XiOa/i1nJnaEiNQW/og4xy8ZxXcWu5MzSkRps37lEvnvMqbi13hobUGPbGvSpZsBsUehW3ljs7wqXGW49Zy3+76X+/YP45bziC/3zWaxdkG17FreXO0JAaX7vtkYHzb7h3z4JuZ7lcxe3QYA1iaEj03iCfeObnA5e9GDupHRqsYTrv00hyepL7kuxKcsGA5S9JckWz/OYkr+q6Jmmm2Ya8DuvruOo7U7z5E9/kyAuu5s2f+OaK+sY8hwZrmE5DI8kq4DPAGcCxwDlJjp3R7H3Aj6rqNcCngU92WZM0yGxHE4M6qVf6V606NFjDdH2ksQHYVVUPVNXPgC8DZ85ocybwhWb6b4CTk6TjuqTnGXY0ceA/3Xfg6ZiV/kncocEapuvQWA881Pd4dzNvYJuqehZ4EnjZzBUl2ZxkMsnknj0L2zEpDRsK++F/fdzA9nN9Eh/nqasutuXQYA2zbDrCq2oLsAVgYmJi8P0dpHna26Gw69asZmpAcKxbs3pgJ/K/vWIn51+xk/VrVvPWY9Zyw717FmRUUlcd1g4N1jBdh8YUcETf48ObeYPa7E6yD/BS4Icd16UXibmGjc53WOmHTjv6eW/W8I+fxAedupr+lDP1xDPPuxZk1Df52U6TjfoGv1yGBmu8uj49dQtwVJIjk+wHbAK2zmizFTi3mT4b+GbVkDvFSXthrs7qUTqzzzphPR9/92tZv2b1C64W39vO4lH6Quyw1rh1eqRRVc8mOQ+4FlgFXF5VdyX5KDBZVVuBzwN/mWQX8Di9YJFGNten8FE/pQ/7JD7s1NVs5vsmP9tpMqkLnV+nUVXbquqfVdWvVNV/aeZd3AQGVfV/q+rXq+o1VbWhqh7ouia9OMz1KbyrT+mDOpHnUjCvTmw7rDVu3rBQK9Zcw0a7Glbaf+oKoO348flc6zHbaTKpC1mO3QcTExM1OTm52GVoiZs5sgh6n8Kn31TnWj5offPpNJ/5vOnRU8NOYa1fs5q/v+Bt8/gXS7NLsqOqJkZZx7IZcivtrbmGje7NsNJRhrYO6/s48oKrGfSRzU5sLWWGhla0uYaNzrV8+ihh0FHBqENb7cTWcmSfhjRE/5DcYUY5KrATW8uRRxrSEG2+M3yUowKvutZyZGhIQ8x1FDHsqGBvOsy96lrLjaEhDTHbRXrrh4SBX16klc4+DWmIYX0Of/ybx/P3F7ztRXnLdMkjDWmI+fQ5eC8orXSGhjSLve1zcBitVjpPT0kLyGG0Wuk80pAWkMNotdJ57ympI/O9V5XUFe89JS1RDr3VSmWfhtQBh95qpTI0pA449FYrlaEhdaCrL3iSFlsnoZHkkiT3Jrk9yZVJ1gxp92CSO5LsTGLPtlYMh95qperqSGM78C+q6l8C/wBcOEvbt1bV8aP26EtLiV/DqpWqk9FTVXVd38ObgLO72I60lHkHW61E4+jTeC9wzZBlBVyXZEeSzbOtJMnmJJNJJvfs2bPgRUqS5jbvI40k1wOHDlh0UVV9tWlzEfAs8MUhqzmpqqaSvBzYnuTeqrpxUMOq2gJsgd7FffOtW5I0f/MOjao6ZbblSd4DvAM4uYZcdl5VU83vx5JcCWwABoaGJGnxdTV66nTg94F3VtVPhrTZP8kB09PAqcCdXdQjSVoYXfVpXAocQO+U084klwEkWZdkW9PmEOBbSW4Dvg1cXVVf76geSdIC6Gr01GuGzH8Y2NhMPwC8rovtS5K64RXhkqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmudhUaSjySZar4jfGeSjUPanZ7kviS7klzQVT2SpNF18h3hfT5dVX84bGGSVcBngLcDu4Fbkmytqrs7rkuSNA+LfXpqA7Crqh6oqp8BXwbOXOSaJElDdB0a5yW5PcnlSQ4csHw98FDf493NvBdIsjnJZJLJPXv2dFGrJGkOI4VGkuuT3Dng50zgs8CvAMcDjwCfGmVbVbWlqiaqamLt2rWjrEqSNE8j9WlU1Slt2iX5M+BrAxZNAUf0PT68mSdJWoK6HD11WN/DdwF3Dmh2C3BUkiOT7AdsArZ2VZMkaTRdjp76gyTHAwU8CPwuQJJ1wOeqamNVPZvkPOBaYBVweVXd1WFNkqQRdBYaVfVbQ+Y/DGzse7wN2NZVHZKkhbPYQ24lScuIoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJaq2Tr3tNcgVwdPNwDfBEVR0/oN2DwI+B54Bnq2qii3okSQujk9Coqt+cnk7yKeDJWZq/tap+0EUdkqSF1UloTEsS4DeAt3W5HUnSeHTdp/GrwKNVdf+Q5QVcl2RHks0d1yJJGtG8jzSSXA8cOmDRRVX11Wb6HOBLs6zmpKqaSvJyYHuSe6vqxiHb2wxsBnjFK14x37IlSSNIVXWz4mQfYAp4fVXtbtH+I8DTVfWHc7WdmJioycnJ0YuUpBeRJDtGHXDU5empU4B7hwVGkv2THDA9DZwK3NlhPZKkEXUZGpuYcWoqybok25qHhwDfSnIb8G3g6qr6eof1SJJG1Nnoqap6z4B5DwMbm+kHgNd1tX1J0sLzinBJUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktSaoSFJas3QkCS1ZmhIklozNCRJrRkakqTWDA1JUmuGhiSpNUNDktTaSKGR5NeT3JXk/yWZmLHswiS7ktyX5LQhzz8yyc1NuyuS7DdKPZKkbo16pHEn8G7gxv6ZSY4FNgHHAacDf5pk1YDnfxL4dFW9BvgR8L4R65EkdWik0Kiqe6rqvgGLzgS+XFU/rar/BewCNvQ3SBLgbcDfNLO+AJw1Sj2SpG7t09F61wM39T3e3czr9zLgiap6dpY2v5BkM7C5efjTJHcuUK1dOhj4wWIXMYflUCNY50KzzoW1XOo8etQVzBkaSa4HDh2w6KKq+uqoBbRVVVuALU1Nk1U1McdTFt1yqHM51AjWudCsc2EtpzpHXcecoVFVp8xjvVPAEX2PD2/m9fshsCbJPs3RxqA2kqQlpKsht1uBTUlekuRI4Cjg2/0NqqqAG4Czm1nnAmM7cpEk7b1Rh9y+K8lu4E3A1UmuBaiqu4C/Au4Gvg58sKqea56zLcm6ZhX/Efh3SXbR6+P4fMtNbxml7jFaDnUuhxrBOheadS6sF02d6X3glyRpbl4RLklqzdCQJLW2ZENjud2ipNnGzubnwSQ7h7R7MMkdTbuRh7/trSQfSTLVV+vGIe1Ob/bvriQXLEKdlyS5N8ntSa5MsmZIu0XZn3Ptn2YQyBXN8puTvGpctfXVcESSG5Lc3fxf+r0Bbd6S5Mm+18PF466zqWPWv2N6/qTZn7cnOXHM9R3dt492Jnkqyfkz2izavkxyeZLH+q9fS3JQku1J7m9+Hzjkuec2be5Pcu6cG6uqJfkD/HN6F6L8D2Cib/6xwG3AS4Ajge8CqwY8/6+ATc30ZcAHxlj7p4CLhyx7EDh4EffrR4D/MEebVc1+fTWwX7O/jx1znacC+zTTnwQ+uVT2Z5v9A/wb4LJmehNwxSL8rQ8DTmymDwD+YUCdbwG+Nu7a9vbvCGwErgECvBG4eRFrXQV8H3jlUtmXwK8BJwJ39s37A+CCZvqCQf+HgIOAB5rfBzbTB862rSV7pFHL9BYlzbZ/A/jSOLbXkQ3Arqp6oKp+BnyZ3n4fm6q6rv7xbgE30buOZ6los3/OpPe6g97r8OTmtTE2VfVIVd3aTP8YuIdZ7rqwxJ0J/EX13ETvGq/DFqmWk4HvVtX3Fmn7L1BVNwKPz5jd/xoc9h54GrC9qh6vqh8B2+ndL3CoJRsas1gPPNT3eORblCywXwUerar7hywv4LokO5pboyyG85pD/MuHHLK22cfj9F56nzIHWYz92Wb//KJN8zp8kt7rclE0p8dOAG4esPhNSW5Lck2S48Zb2S/M9XdcSq/JTQz/ULgU9uW0Q6rqkWb6+8AhA9rs9X7t6t5TrWSJ3KKkrZb1nsPsRxknVdVUkpcD25Pc23xKGEudwGeBj9H7T/oxeqfS3ruQ22+rzf5MchHwLPDFIavpfH8ud0l+Gfhb4PyqemrG4lvpnWZ5uunfuorexbjjtiz+jk3f6DuBCwcsXir78gWqqpIsyPUVixoatcxuUTJXvUn2oXer+NfPso6p5vdjSa6kd6pjQf9ztN2vSf4M+NqARW328cha7M/3AO8ATq7mBOyAdXS+Pwdos3+m2+xuXhcvpfe6HKsk+9ILjC9W1VdmLu8PkaraluRPkxxcVWO9+V6Lv+NYXpMtnAHcWlWPzlywVPZln0eTHFZVjzSn8h4b0GaKXl/MtMPp9SMPtRxPTy3lW5ScAtxbVbsHLUyyf5IDpqfpdfaO9W69M84Dv2vI9m8BjkpvBNp+9A7Ht46jvmlJTgd+H3hnVf1kSJvF2p9t9s9Weq876L0Ovzks+LrS9KF8Hrinqv5oSJtDp/takmyg954w1nBr+XfcCvx2M4rqjcCTfadexmnomYSlsC9n6H8NDnsPvBY4NcmBzanqU5t5wy1GT3/L0QDvond+7afAo8C1fcsuojd65T7gjL7524B1zfSr6YXJLuCvgZeMoeY/B94/Y946YFtfTbc1P3fROw0z7v36l8AdwO3Ni+qwmXU2jzfSG23z3UWqcxe9c607m5/LZta5mPtz0P4BPkov5AD+SfO629W8Dl+9CPvwJHqnIW/v248bgfdPv06B85p9dxu9AQf/ahHqHPh3nFFngM80+/sO+kZUjrHO/emFwEv75i2JfUkvyB4Bft68b76PXh/aN4D7geuBg5q2E8Dn+p773uZ1ugv4nbm25W1EJEmtLcfTU5KkRWJoSJJaMzQkSa0ZGpKk1gwNSVJrhoYkqTVDQ5LU2v8HstP07LyucZAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "l, = ax.plot(split(res.y[:,0])[0], split(res.y[:,0])[1], 'o')\n", "ax.set_xlim(-10, 10)\n", "ax.set_ylim(-10, 10)\n", "\n", "def animate(i):\n", " r_x, r_y, v_x, v_y = split(res.y[:,i])\n", " l.set_data(r_x, r_y)\n", "\n", "ani = FuncAnimation(fig, \n", " animate, \n", " frames=len(res.t),\n", " interval=100)\n", "\n", "\n", "HTML(ani.to_jshtml());" ] }, { "cell_type": "code", "execution_count": null, "id": "98262899", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 39, "id": "2cfdbc54", "metadata": {}, "outputs": [], "source": [ "# Question 8\n", "\n", "from numba import jit\n", "\n", "alpha = 0.2\n", "\n", "@jit\n", "def join_numba(r_x, r_y, v_x, v_y):\n", " return np.concatenate((r_x, r_y, v_x, v_y))\n", "\n", "@jit\n", "def split_numba(y):\n", " N = len(y)//4\n", " return y[:N], y[N:2*N], y[2*N:3*N], y[3*N:4*N]\n", "\n", "@jit\n", "def force_piege_numba(r_x, r_y):\n", " return -k_x*r_x, -k_y*r_y\n", "\n", "@jit\n", "def force_coulomb_numba(r_x, r_y):\n", " N = len(r_x)\n", " f_x = np.zeros(N)\n", " f_y = np.zeros(N)\n", " for i in range(N):\n", " tot_x, tot_y = 0, 0\n", " for j in range(N):\n", " if i!=j:\n", " d2 = (r_x[i]-r_x[j])**2 + (r_y[i]-r_y[j])**2\n", " coef = kappa*q**2/(np.sqrt(d2)**3)\n", " tot_x += coef*(r_x[i]-r_x[j])\n", " tot_y += coef*(r_y[i]-r_y[j])\n", " f_x[i] = tot_x\n", " f_y[i] = tot_y\n", " return f_x, f_y\n", "\n", "\n", "@jit\n", "def force_friction_numba(r_x, r_y, v_x, v_y):\n", " return -alpha*v_x, 0*v_y\n", "\n", "@jit\n", "def force_totale_numba(r_x, r_y, v_x, v_y):\n", " piege = force_piege_numba(r_x, r_y)\n", " coulomb = force_coulomb_numba(r_x, r_y)\n", " friction = force_friction_numba(r_x, r_y, v_x, v_y)\n", " return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]\n", "\n", "\n", "@jit\n", "def f_numba(t, y):\n", " r_x, r_y, v_x, v_y = split_numba(y)\n", " f_x, f_y = force_totale_numba(r_x, r_y, v_x, v_y)\n", " return join_numba(v_x, v_y, f_x/m, f_y/m)" ] }, { "cell_type": "code", "execution_count": 40, "id": "e511c564", "metadata": {}, "outputs": [], "source": [ "N = 20\n", "\n", "r_x_0 = np.random.normal(size=N)\n", "r_y_0 = np.random.normal(size=N)\n", "v_x_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "v_y_0 = np.random.normal(scale=1/np.sqrt(m), size=N)\n", "\n", "t_max = 1500\n", "t_eval=np.linspace(0, t_max, 400)\n", "\n", "# Ici, on apelle f_numba un première fois pour la compilation\n", "f_numba(0.0, join(r_x_0, r_y_0, v_x_0, v_y_0));" ] }, { "cell_type": "code", "execution_count": 41, "id": "c2d145ec", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 168 ms, sys: 0 ns, total: 168 ms\n", "Wall time: 167 ms\n" ] } ], "source": [ "%%time\n", "res = solve_ivp(f_numba, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0), \n", " t_eval=t_eval)" ] }, { "cell_type": "code", "execution_count": null, "id": "dfb78d7a", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "55a0d2f2", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 44, "id": "39df2180", "metadata": {}, "outputs": [], "source": [ "# Question 9\n", "particle_type = np.array([0]*20 + [1]*5)\n", "masses = np.array([9, 3]) # The masse of particles given by masses[particle_type]\n", "colors = np.array(['blue', 'red'])\n", "\n", "@jit\n", "def force_friction_numba(r_x, r_y, v_x, v_y, particle_type):\n", " return -alpha*v_x*(particle_type==0), 0*v_y\n", "\n", "@jit\n", "def force_totale_numba(r_x, r_y, v_x, v_y, particle_type):\n", " piege = force_piege_numba(r_x, r_y)\n", " coulomb = force_coulomb_numba(r_x, r_y)\n", " friction = force_friction_numba(r_x, r_y, v_x, v_y, particle_type)\n", " return piege[0]+coulomb[0]+friction[0], piege[1]+coulomb[1]+friction[0]\n", "\n", "@jit\n", "def f_numba(t, y, particle_type):\n", " r_x, r_y, v_x, v_y = split_numba(y)\n", " f_x, f_y = force_totale_numba(r_x, r_y, v_x, v_y, particle_type)\n", " return join_numba(v_x, v_y, f_x/masses[particle_type], f_y/masses[particle_type])" ] }, { "cell_type": "code", "execution_count": 45, "id": "64f87801", "metadata": {}, "outputs": [], "source": [ "N = len(masses[particle_type])\n", "r_x_0 = np.random.normal(size=N)\n", "r_y_0 = np.random.normal(size=N)\n", "v_x_0 = np.random.normal(scale=1/np.sqrt(masses[particle_type]), size=N)\n", "v_y_0 = np.random.normal(scale=1/np.sqrt(masses[particle_type]), size=N)\n", "\n", "t_max = 1500\n", "t_eval=np.linspace(0, t_max, 400)\n", "res = solve_ivp(f_numba, [0, t_max], join(r_x_0, r_y_0, v_x_0, v_y_0), \n", " t_eval=t_eval, args=(particle_type,))\n" ] }, { "cell_type": "code", "execution_count": 46, "id": "020cf474", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "
\n", " \n", "
\n", " \n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdHElEQVR4nO3de5BcdZ338fcnmVwnWQgYEiAEQV2QiwQyRqPBRcEAWSTgIga3BJSqgA+4suujgpTIyl68sVQJaiqICi5y0wWihIQg7KLlAzKBhEsgEBCXxJgECCEXcv8+f/zOkE5Pd8/JdJ+eGfJ5VXXNufz6nG/OdPo753c7igjMzMzy6NfTAZiZWd/hpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuTUkaUj6saSVkp4s2baXpHmSnst+jqjy3nOyMs9JOqcR8ZiZWTEadafxU+Cksm2XAL+JiHcBv8nWdyJpL+DrwPuACcDXqyUXMzPreQ1JGhHxIPBq2eapwA3Z8g3AaRXeeiIwLyJejYjVwDw6Jx8zM+slWgo89qiIWJ4t/wUYVaHM/sBLJetLs22dSJoOTAdobW0df+ihhzYwVDOzt7758+e/HBEj6zlGkUnjTRERkuqaryQiZgIzAdra2qK9vb0hsZmZ7S4k/aneYxTZe2qFpH0Bsp8rK5RZBhxQsj4m22ZmZr1QkUljFtDRG+oc4K4KZeYCkyWNyBrAJ2fbzMysF2pUl9ubgf8HHCJpqaTzgG8CH5X0HHBCto6kNkk/AoiIV4ErgUey1zeybWZm1gupL06N7jYNM7NdJ2l+RLTVcwyPCDczs9ycNMzMLDcnDTMzy81Jw8zMcnPSMDOz3Jw0zMwsNycNMzPLzUnDzMxyc9IwM7PcnDTMzCw3Jw0zM8vNScPMzHJz0jAzs9ycNMzMLDcnDTMzy81Jw8zMcnPSMDOz3Jw0zMwst0KThqRDJC0oeb0u6eKyMsdJWlNS5vIiYzIzs+5rKfLgEbEYGAcgqT+wDLijQtHfRsQpRcZiZmb1a2b11PHA8xHxpyae08zMGqiZSWMacHOVfRMlLZR0j6TDmxiTmZntgqYkDUkDgVOB2yvsfhQ4MCKOAq4B7qxyjOmS2iW1r1q1qrBYzcysumbdaZwMPBoRK8p3RMTrEbEuW54NDJD0tgrlZkZEW0S0jRw5sviIzcysk2YljbOoUjUlabQkZcsTspheaVJcZma2CwrtPQUgqRX4KHB+ybYLACJiBnAG8DlJW4E3gGkREUXHZWZmu67wpBER64G9y7bNKFm+Fri26DjMzKx+HhFuZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuTlpmJlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmllvhSUPSi5KekLRAUnuF/ZL0PUlLJD0u6ZiiYzIzs+5padJ5PhwRL1fZdzLwruz1PuCH2U8zM+tlekP11FTgxkgeAvaUtG9PB2VmZp01I2kEcK+k+ZKmV9i/P/BSyfrSbNtOJE2X1C6pfdWqVQWFamZmtTQjaUyKiGNI1VAXSvpQdw4SETMjoi0i2kaOHNnYCM3MLJfCk0ZELMt+rgTuACaUFVkGHFCyPibbZmZmvUyhSUNSq6ThHcvAZODJsmKzgLOzXlTvB9ZExPIi4zIzs+4puvfUKOAOSR3n+nlEzJF0AUBEzABmA1OAJcAG4DMFx2RmZt1UaNKIiBeAoypsn1GyHMCFRcZhZmaN0Ru63JqZWR/hpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuTlpmJlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuRWWNCQdIOkBSYskPSXpCxXKHCdpjaQF2evyouIxM7P6FfmM8K3AFyPiUUnDgfmS5kXEorJyv42IUwqMw8zMGqSwO42IWB4Rj2bLa4Gngf2LOp+ZmRWvKW0akt4OHA08XGH3REkLJd0j6fAax5guqV1S+6pVq4oK1czMaig8aUgaBvwSuDgiXi/b/ShwYEQcBVwD3FntOBExMyLaIqJt5MiRhcVr1kxbt8L27T0dhVl+hSYNSQNICeOmiPiv8v0R8XpErMuWZwMDJL2tyJjMeoPFi+G442Dw4PT61Kdg9eqejsqsa4U1hEsScD3wdET8R5Uyo4EVERGSJpCS2CtFxWTWG7z6KkycCK+9BhGwbRv88pewaBE89hhIPR2hWXVF9p76IPBp4AlJC7JtXwXGAkTEDOAM4HOStgJvANMiIgqMyazH/fSnsHFjShgdNm+G55+H3/0Ojj22x0Iz61JhSSMifgfU/JspIq4Fri0qBrPe6Ikn4I03Om+PgGefddKw3s0jws2arK0Nhg6tvO+II5obi9muctIwa7JPfxqGD4f+/XdsGzQIxo2DCRPyH+fBB+G889LxZs/eubrLrChOGmZN9ld/BY88AqedBkOGwB57wPnnw9y5+RvBv/pVOPlk+MlP4D//E848MyUPJw4rmvpiu3NbW1u0t7f3dBhmPeL551M11saNO29vbYU5c2DSpJ6Jy3o/SfMjoq2eY/hOw6yPmTOn8vb16+Guu5obi+1+nDTMmmjjRli+PI3N6K7W1p3bQzoMGJCqvsyK5KRhtoueew6+9jW46KL0V3+eaUC2boV//EfYay84+GAYNQquv7575z/ttMptFy0t8Pd/371jmuXlpGFvHREwYwb89V/DyJHwyU+mBoAGuukmOOoo+OY34fvfhzPOgKlTu75z+OIXYebMND5j40Z45RX4h3+AWbN2PYY994Q774Rhw1IvrOHDU4P6zJkpIZkVyQ3h9tbxxS+mpLFhQ1rv1y/V1zzxBIwZs6Pc5s3wox/BjTemOp3p09Of6P1q/w21bl26Q+g4fIfW1tSL6ROfqPy+jRvTHUalAX3HHAPz5+/Cv7HEhg1w332wZQuccELqhWVWSyMawoucRsSseV55BX7wg527FG3fnr5Zr7oKrr56x7bJk1Of145v/8ceg3vugZ//vOYp/ud/UhVQufXr4eabqyeNV2rMpvbSSzVPWdPQoXDqqd1/v1l3uHrK3hoWLUoj5Mpt3pwmdOowZ0760770dqGj29HChTVPMXBg9X1DhlTfN2pUmsm2kvHja57SrNdx0rC3hrFjYdOmztv79UttHB3uvz/VM5Xbti3dStTwN39TuddSa2samV1NS0tqAymfOmToUPjXf615SrNex0nD3hoOPDA9oKL8bmPwYPjSl3asjx5d+c/+AQNgn31qnuJ//xcOPXTHektLOt3nPw8f+cjOZdetgy98AY4+OjWUv/e9aeT2uHGw995w4onw29+mNo2ibdqUmnWWLy/+XLYbiIg+9xo/fnyYdbJ2bcSnPhUxaFB6jR0bMXv2zmWWL48YOjQi9bXa8dpzz4j166seeuXKiBEjIvr12/GWAQMijjtuR5k774yYODFizJidy0GEFHHNNQX9u2uYOTNi+PD0GjQo4sQTI1avbn4c1jsA7VHn9697T9lbz4YN6U/9kSMrT+Z0//2pO+6mTek7fcSI1KZx9NFVD3nllfBv/9Z56o4hQ+DRR+FXv4Irrujcs6pUSwu8/HLzejn95jepobw0poEDUzXbvfc2JwbrXdx7yqySoUOrzz0OqS7pL39JvaYGDID3vKfLmQL/8IfOCQPS29vbu04YsKPZpFk9nr7znc4xbd6cqsWWLt25F7JZXk4atnvq3z892CKH115Lo8Ar2bo1HapSV9xKWls7b9u+HR54IJ3jiCPggx9szCNfly6tvH3gQFixwknDusdJw6wLp5wCf/xj5+0DBqS8M2lS+gu+K4MHp6qhUi+/DB/6UBqvsW1b6ux1xBFp0N6wYfXFffzx6UmAW7bsvH3bNnj3u+s7tu2+Cu89JekkSYslLZF0SYX9gyTdmu1/WNLbi47JLK+nnkq1WJWSwjveAXffDQcckL74a43j6N8/jR8svyO54AJYsiQ1wbzxRhoysmBBel5Gvb7yldR+MmDAjm0d3Xxr1d6Z1VJo0pDUH/g+cDJwGHCWpMPKip0HrI6IdwJXA98qMiazXfHHP+78pVtq8eKUOP7lX9KI8BNOSF1whw1LX9ZnnQUf/zh8/euwdm3nu4ytW9PcU+V3Aps2wc9+Vn/s++2XEtD556ehKscdB7ffnroCm3VX0dVTE4AlEfECgKRbgKnAopIyU4ErsuVfANdKUvTFbl32ljNuXOUxg5A6Xq1cCf/+7/CnP6W7jlWrUpXTO99ZPdmUvr/aDLnliaS79t8frrmmMccyg+Krp/YHSmfXWZptq1gmIrYCa4C9yw8kabqkdkntq1atKihce6vaujVVy4wenapmJk9OVU9dGTMm3THUqs7ZsCFNWHj33amX77vf3XXCgFSmUqN3Swt87GNdv7+WFSvg979PP80aqc+MCI+ImRHRFhFtI0eO7OlwrI+ZPj2Ns1ixIrUdzJsHEyemO4SuXHddSjgHHVR5GhFIjct/93fpHOXmzIFjj02D1s86KzVOlx57xIgdSam1NQ1Mv+qqXf83QkqOn/lMOteUKennuec27s7FrNCR28BEYG7J+qXApWVl5gITs+UW4GWyKdurvTwi3HbF8uVpNHT5IPABAyI+//ms0Pr1EV/7WsSBB6aR5JddFrFuXadjnXZa5+OUvgYPjli2bEf566/feQB6v35pdPYzz+wos3p1xPe+F3H++RHXXVfxtLlddlnEkCE7xzRkSNpuRgNGhBedNFqAF4CDgIHAQuDwsjIXAjOy5WnAbV0d10nDdsUDD0TssUflL/n3vS8itm2LGDcuon//HTsGDowYPz7tKzF+fO2kMWxYxA03pLJbtqTZScrL9OsXceaZxfxbK52vY5YUs0YkjUKrpyK1UVyU3U08nSWEpyR9Q1LHuNjrgb0lLQH+CejULdesHgcfXLkxu6UFjjyS1CCxYMHOj9/bvDlNtz5nzpubtm1L3W9rkXaMr1i6tHK10PbtO8/W3kivv75r2812VeGD+yJiNjC7bNvlJcsbgSqPrzGr39ixcPLJaZxE6VQgAwemh/0x9duV3/jGG2mOkClTgDTwrqWl9kA+KZ0L0my2pXnoMJ7iXH7KHqxhUetpsP2kLp8WuKva2tKUJ+X83A5rlD7TEG5Wj5//PDUQDxmSvqePPDJN2nfoocCLL1Z/4+jRby5KaZ7DSs96Gjgwjc24++4dD2QaPhzOPDOtn8d1PMJ7uZirmc51fPelT6Y506v1ue2ma69NjekdDfb9+6f1a69t6GlsN+ZZbm23sn17qjLa6Yt/1Kg04KKSNWvSc8ZLVju660rpeAcfDP/8z/C3f9s5oWzcCBefs5qrb9uPIZTNeDhsWHpO+emnN+Yfl1m8GL797VSVNm5cGhl+yCENPYX1UZ7l1mwX9etX4U7hoos6z3supT60JQkD0t3EQw+lR4w/+ywcfnjNGdUZPBhmfPIBtt8zENaWJY116+C22xqeNA45BK6/vqGHNHuTk4bZpZemDPCLX6R6ps2b4cMfhhtueLPIY4+lSQT32APOOAMmTEivXAYPpl+lWWslTwJlfY6rp8w6LF2aeky94x3pReqwet55cOutqVprwID0XT9rVudHvFa1aVNqG3nttZ23Dx2aRhl+4AMN/WeYVdOI6ik3hJt1GDMmNVhkCQNScrjttjRVyJYt6ef69WkiwjzToQOpPuxXv0pVXcOHp7aMQYPSHU6dCeP3v0/t6UcemWrZXnqp6/eY1cPVU2Y1/OQnKUmUi0hPwDv++JwHmjQJli9P/X7XroWPfjTNJliH229PU4S88UaKZ/FiuOmm9PjZgw6q69BmVTlpmNVQq0fsLtfsDh2aGtcbYNs2uPDCnR/numVLGsR3+eWNmVrdrBJXT5nVcPbZlR/RCmkSwp6ydGnlO6Dt2+H++5sfj+0+nDTMavj4x9PjXltbUwP44MHphuGWWyoP8iva9u1pUOJ111VvU9lnn+bGZLsXV0+Z1dCvX3oq30MPwdy5aRrzadPSeMBm27Ah9QRetCgtlz+HA1Jy+/KXmx+b7T6cNMy6IKVnb0yc2LNxfPOb8PjjO49BhBTf8OGpTeNLX0pJzawoThpmfcSNN3ZOGJDGjtx6a3oK4PDhzY/Ldi9u0zDr46Q0lYkThjWDk4ZZH3H22akhvpSUnkneE20stnty0jDrIy65JI38HjZsx8Oe9torTftu1ixu0zDrI4YOTb245s1LD1oaOzZNnlhtHIlZEZw0zPqQfv3gxBPTy6wnFJI0JH0H+BiwGXge+ExEvFah3IvAWmAbsLXe2RfNzKxYRbVpzAOOiIj3AM8Cl9Yo++GIGOeEYWbW+xWSNCLi3ojYmq0+BIwp4jxmZtZczeg99Vngnir7ArhX0nxJ02sdRNJ0Se2S2letWtXwIM3MrGvdbtOQdB8wusKuyyLirqzMZcBW4KYqh5kUEcsk7QPMk/RMRDxYqWBEzARmQnpyX3fjNjOz7ut20oiIE2rtl3QucApwfFR5pmxELMt+rpR0BzABqJg0zMys5xVSPSXpJODLwKkRsaFKmVZJwzuWgcnAk0XEY2ZmjVFUm8a1wHBSldMCSTMAJO0naXZWZhTwO0kLgT8Ad0fEnILiMTOzBihknEZEvLPK9j8DU7LlF4Cjiji/mZkVw3NPmZlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuTlpmJlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluThpmZpZbYUlD0hWSlmXPCF8gaUqVcidJWixpiaRLiorHzMzqV8gzwktcHRHfrbZTUn/g+8BHgaXAI5JmRcSiguMyM7Nu6OnqqQnAkoh4ISI2A7cAU3s4JjMzq6LopHGRpMcl/VjSiAr79wdeKllfmm3rRNJ0Se2S2letWlVErGZm1oW6koak+yQ9WeE1Ffgh8A5gHLAcuKqec0XEzIhoi4i2kSNH1nMoMzPrprraNCLihDzlJF0H/LrCrmXAASXrY7JtZmbWCxXZe2rfktXTgScrFHsEeJekgyQNBKYBs4qKyczM6lNk76lvSxoHBPAicD6ApP2AH0XElIjYKukiYC7QH/hxRDxVYExmZlaHwpJGRHy6yvY/A1NK1mcDs4uKw8zMGqenu9yamVkf4qRhZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuTlpmJlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluThpmZpabk4aZmeXmpGFmZrkV8rhXSbcCh2SrewKvRcS4CuVeBNYC24CtEdFWRDxmZtYYhSSNiPhkx7Kkq4A1NYp/OCJeLiIOMzNrrEKSRgdJAs4EPlLkeczMrDmKbtM4FlgREc9V2R/AvZLmS5pecCxmZlanbt9pSLoPGF1h12URcVe2fBZwc43DTIqIZZL2AeZJeiYiHqxyvunAdICxY8d2N2wzM6uDIqKYA0stwDJgfEQszVH+CmBdRHy3q7JtbW3R3t5ef5BmZrsRSfPr7XBUZPXUCcAz1RKGpFZJwzuWgcnAkwXGY2ZmdSoyaUyjrGpK0n6SZmero4DfSVoI/AG4OyLmFBiPmZnVqbDeUxFxboVtfwamZMsvAEcVdX4zM2s8jwg3M7PcnDTMzCw3Jw0zM8vNScPMzHJz0jAzs9ycNMzMLDcnDTMzy81Jw8zMcnPSMDOz3Jw0zMwsNycNMzPLzUnDzMxyc9IwM7PcnDTMzCw3Jw0zM8vNScPMzHJz0jAzs9ycNMzMLDcnDTMzy62upCHpE5KekrRdUlvZvkslLZG0WNKJVd5/kKSHs3K3ShpYTzxmZlaseu80ngQ+DjxYulHSYcA04HDgJOAHkvpXeP+3gKsj4p3AauC8OuMxM7MC1ZU0IuLpiFhcYddU4JaI2BQRfwSWABNKC0gS8BHgF9mmG4DT6onHzMyK1VLQcfcHHipZX5ptK7U38FpEbK1R5k2SpgPTs9VNkp5sUKxFehvwck8H0YW+ECM4zkZznI3VV+I8pN4DdJk0JN0HjK6w67KIuKveAPKKiJnAzCym9oho6+ItPa4vxNkXYgTH2WiOs7H6Upz1HqPLpBERJ3TjuMuAA0rWx2TbSr0C7CmpJbvbqFTGzMx6kaK63M4CpkkaJOkg4F3AH0oLREQADwBnZJvOAZp252JmZruu3i63p0taCkwE7pY0FyAingJuAxYBc4ALI2Jb9p7ZkvbLDvEV4J8kLSG1cVyf89Qz64m7ifpCnH0hRnCcjeY4G2u3iVPpD34zM7OueUS4mZnl5qRhZma59dqk0demKMnOsSB7vShpQZVyL0p6IitXd/e3XSXpCknLSmKdUqXcSdn1XSLpkh6I8zuSnpH0uKQ7JO1ZpVyPXM+urk/WCeTWbP/Dkt7erNhKYjhA0gOSFmX/l75QocxxktaUfB4ub3acWRw1f49Kvpddz8clHdPk+A4puUYLJL0u6eKyMj12LSX9WNLK0vFrkvaSNE/Sc9nPEVXee05W5jlJ53R5sojolS/g3aSBKP8NtJVsPwxYCAwCDgKeB/pXeP9twLRseQbwuSbGfhVweZV9LwJv68HregXwf7so0z+7rgcDA7PrfViT45wMtGTL3wK+1VuuZ57rA/wfYEa2PA24tQd+1/sCx2TLw4FnK8R5HPDrZse2q79HYApwDyDg/cDDPRhrf+AvwIG95VoCHwKOAZ4s2fZt4JJs+ZJK/4eAvYAXsp8jsuURtc7Va+80oo9OUZKd+0zg5macryATgCUR8UJEbAZuIV33pomIe2PHbAEPkcbx9BZ5rs9U0ucO0ufw+Oyz0TQRsTwiHs2W1wJPU2PWhV5uKnBjJA+Rxnjt20OxHA88HxF/6qHzdxIRDwKvlm0u/QxW+w48EZgXEa9GxGpgHmm+wKp6bdKoYX/gpZL1uqcoabBjgRUR8VyV/QHcK2l+NjVKT7gou8X/cZVb1jzXuJk+S/ors5KeuJ55rs+bZbLP4RrS57JHZNVjRwMPV9g9UdJCSfdIOry5kb2pq99jb/pMTqP6H4W94Vp2GBURy7PlvwCjKpTZ5eta1NxTuaiXTFGSV854z6L2XcakiFgmaR9gnqRnsr8SmhIn8EPgStJ/0itJVWmfbeT588pzPSVdBmwFbqpymMKvZ18naRjwS+DiiHi9bPejpGqWdVn71p2kwbjN1id+j1nb6KnApRV295Zr2UlEhKSGjK/o0aQRfWyKkq7ildRCmip+fI1jLMt+rpR0B6mqo6H/OfJeV0nXAb+usCvPNa5bjut5LnAKcHxkFbAVjlH49awgz/XpKLM0+1zsQfpcNpWkAaSEcVNE/Ff5/tIkEhGzJf1A0tsioqmT7+X4PTblM5nDycCjEbGifEdvuZYlVkjaNyKWZ1V5KyuUWUZqi+kwhtSOXFVfrJ7qzVOUnAA8ExFLK+2U1CppeMcyqbG3qbP1ltUDn17l/I8A71LqgTaQdDs+qxnxdZB0EvBl4NSI2FClTE9dzzzXZxbpcwfpc3h/tcRXlKwN5Xrg6Yj4jyplRne0tUiaQPpOaGpyy/l7nAWcnfWiej+wpqTqpZmq1iT0hmtZpvQzWO07cC4wWdKIrKp6cratup5o6c/ZG+B0Uv3aJmAFMLdk32Wk3iuLgZNLts8G9suWDyYlkyXA7cCgJsT8U+CCsm37AbNLYlqYvZ4iVcM0+7r+DHgCeDz7UO1bHme2PoXU2+b5HopzCamudUH2mlEeZ09ez0rXB/gGKckBDM4+d0uyz+HBPXANJ5GqIR8vuY5TgAs6PqfARdm1W0jqcPCBHoiz4u+xLE4B38+u9xOU9KhsYpytpCSwR8m2XnEtSYlsObAl+948j9SG9hvgOeA+YK+sbBvwo5L3fjb7nC4BPtPVuTyNiJmZ5dYXq6fMzKyHOGmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmltv/ByjCjOkN208gAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots()\n", "X = split(res.y[:,0])[0]\n", "Y = split(res.y[:,0])[1]\n", "l = ax.scatter(X, Y, c=colors[particle_type])\n", "ax.set_xlim(-10, 10)\n", "ax.set_ylim(-10, 10)\n", "\n", "def animate(i):\n", " r_x, r_y, v_x, v_y = split(res.y[:,i])\n", " off = np.array([r_x, r_y]).T\n", " l.set_offsets(off)\n", "\n", "ani = FuncAnimation(fig, \n", " animate, \n", " frames=len(res.t),\n", " interval=100)\n", "\n", "\n", "HTML(ani.to_jshtml())" ] }, { "cell_type": "code", "execution_count": null, "id": "155a4bb3", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "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.8" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 5 }