{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Greičio integravimas Euler ir RK4 metodais\n", "Pradiniu laiko momentu $t_0=0$ kūnas buvo koordinačių pradžios taške $x_0=0$. Jei jo greičio nuo laiko priklausomybė yra $v(t)= t^4$, raskite ir pavaizduokite atsumo priklausomybė nuo laiko iki laiko momento $t_1=1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pirmiausia nubrėžkime greičio grafiką" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfjElEQVR4nO3deXzU1b3/8dcneyCEEMISEwKI7AgKYXHHWhWwlmvdtXWp1dpqt9tfW6+9tZu9rb3db1VEitZWodVaxRartVURUTZlX0OAJIAkJBASQpaZOb8/EmwaAjPATL6zvJ+Pxzwy3/mezHwOk7xzOPP9fo855xARkdiX5HUBIiISHgp0EZE4oUAXEYkTCnQRkTihQBcRiRMKdBGROBE00M1srplVmtm6Y+w3M/uVmZWY2RozGx/+MkVEJJhQRuhPAtOOs386MLTtdhfw6KmXJSIiJyolWAPn3CIzG3ScJjOBp1zrGUrvmlmOmeU75/Yc73nz8vLcoEHHe1oREelo5cqV+5xzfTrbFzTQQ1AAlLfbrmh77LiBPmjQIFasWBGGlxcRSRxmtvNY+8Lxoah18lin1xMws7vMbIWZraiqqgrDS4uIyBHhCPQKYEC77UJgd2cNnXOznXPFzrniPn06/R+DiIicpHAE+gLglrajXaYAtcHmz0VEJPyCzqGb2TxgKpBnZhXAt4FUAOfcLGAhMAMoARqA2yNVrIiIHFsoR7ncGGS/A+4JW0UiInJSdKaoiEicUKCLiMQJBbqISBf65WtbWVpaHZHnVqCLiHSRsuoGfv7aFpZtr4nI8yvQRUS6yB9XlJNkcE1xYUSeX4EuItIFfP4Az64sZ+rwvuT3zIzIayjQRUS6wJtbqth7sInrJw4I3vgkKdBFRLrA/OXl5GWl85ERfSP2Ggp0EZEIqzzYyD83VXLNhEJSkyMXuwp0EZEIe+69CvwBF9HpFlCgi4hEVCDg+OPyciYPzmVwXveIvpYCXUQkgt4prWZHdQM3TIrs6BwU6CIiEfXM0jJyuqUyfUx+xF9LgS4iEiFVdU28sv4Drh5fSEZqcsRfT4EuIhIhz64sxxdw3DipqEteT4EuIhIBgYBj3rIyppyeyxl9s7rkNRXoIiIR8FbJPsprDnPT5IFd9poKdBGRCHhm6U5yu6dx+eh+XfaaCnQRkTDbe7CR1zZWcu2EQtJTIv9h6BEKdBGRMHtmaRn+gOOmyV3zYegRCnQRkTBq8QeYt6yMi4b1YWDvyJ4Z2pECXUQkjP6+YS+VdU3cck7XfRh6hAJdRCSMnnpnBwU5mUwdHrnL5B6LAl1EJEy27q3j3dIaPjllIMlJ1uWvr0AXEQmT3727k7TkJK6L0JqhwSjQRUTCoL7Jx/Pv7eKKsfn0zkr3pAYFuohIGDz/XgX1TT4+5cGHoUco0EVETlEg4Pjtkh2MK+zJ2QNyPKtDgS4icooWl+xjW9UhbjtvEGZd/2HoEQp0EZFT9OSSHeRlpTPjzMgvYnE8CnQRkVOwY98hXt9cyU2Ti7r0ui2dUaCLiJyCp97ZSbIZn+zi67Z0JqRAN7NpZrbZzErM7L5O9vc0s5fMbLWZrTez28NfqohIdKlv8vHsinKuGJtP3+wMr8sJHuhmlgw8DEwHRgE3mtmoDs3uATY458YBU4GfmllamGsVEYkqz60op67Jx63nDvK6FCC0EfokoMQ5V+qcawbmAzM7tHFAD2v9eDcLqAF8Ya1URCSK+AOOJ5bsYHxRDuOLenldDhBaoBcA5e22K9oea+/XwEhgN7AW+JJzLhCWCkVEotA/Nu5lZ3UDd5x/utelfCiUQO/soErXYftyYBVwGnAW8Gszyz7qiczuMrMVZraiqqrqBEsVEYkecxZvpyAns0uXmAsmlECvAAa02y6kdSTe3u3A865VCbAdGNHxiZxzs51zxc654j59+pxszSIinlq3q5Zl22u4/bxBpCRHz8GCoVSyHBhqZoPbPui8AVjQoU0ZcAmAmfUDhgOl4SxURCRa/GbxdrqnJXPdxAHBG3ehlGANnHM+M7sXeAVIBuY659ab2d1t+2cB3weeNLO1tE7RfMM5ty+CdYuIeGLvwUZeWr2bT50zkOyMVK/L+TdBAx3AObcQWNjhsVnt7u8GLgtvaSIi0eeJt3cQcI7bzx3sdSlHiZ7JHxGRKFff5OPppTuZPiafot7dvC7nKAp0EZEQzV9WRl2jj7sujJ5DFdtToIuIhKDFH2Du4u1MHpzLOA+veX48CnQRkRD8Zc1udtc28tmLonN0Dgp0EZGgnHM89mYpQ/tmMXVYX6/LOSYFuohIEG9t3cemD+q488LTSUrybkWiYBToIiJBPPrGNvplpzPzrNO8LuW4FOgiIsfxftl+3imt5jPnn+75ikTBKNBFRI7j0Te20TMzlRujYEWiYBToIiLHsHVvHa9u2Mut5w4iKz2kE+s9pUAXETmGR9/cRmZqMrdHyYpEwSjQRUQ6UbG/gQWrdnPjpCJ6dY+NFTUV6CIinZi9qBQz+MwF0XcRrmNRoIuIdFB5sJH5y8u5enwhp+Vkel1OyBToIiIdPLaoFH/A8fmpZ3hdyglRoIuItFNd38TTS3cyc9xpUXmJ3ONRoIuItDNn8XaafAE+f3Fsjc5BgS4i8qEDDc08tWQHV5yZzxl9s7wu54Qp0EVE2sx9eweHmv3c+5HYG52DAl1EBIDahhaeWLydaaP7M6J/ttflnBQFuogIMGdxKXVNPr700aFel3LSFOgikvAONDTzxNs7mHFmf0bmx+boHBToIiI8/lYph5p9fOmSYV6XckoU6CKS0GoONfPk2zuYcWY+w/v38LqcU6JAF5GE9vhbpTS0+PnyJbE7d36EAl1EElZlXSNPvr2DK8eextB+sT06BwW6iCSwR17fRrM/wFcuje258yMU6CKSkHYdOMwzS8u4dkIhg/O6e11OWCjQRSQh/eq1rQB8MQ7mzo9QoItIwimtque59yq4eUpRTF3vPBgFuogknJ/9fQvpKUkxd73zYBToIpJQ1lbU8pc1e/j0eYPp0yPd63LCSoEuIgnlob9tole3VD570elelxJ2IQW6mU0zs81mVmJm9x2jzVQzW2Vm683szfCWKSJy6t7aWsXikn184SND6ZGR6nU5YZcSrIGZJQMPA5cCFcByM1vgnNvQrk0O8AgwzTlXZmZ9I1SviMhJCQQcP3p5E4W9Mrl5SpHX5UREKCP0SUCJc67UOdcMzAdmdmhzE/C8c64MwDlXGd4yRUROzUtrdrN+90G+etkw0lOSvS4nIkIJ9AKgvN12Rdtj7Q0DepnZG2a20sxu6eyJzOwuM1thZiuqqqpOrmIRkRPU5PPzk1c3MzI/m5njOsZX/Agl0K2Tx1yH7RRgAnAFcDnwLTM76lxa59xs51yxc664T58+J1ysiMjJeGrJTsprDnP/jBEkJXUWafEh6Bw6rSPyAe22C4HdnbTZ55w7BBwys0XAOGBLWKoUETlJ+w8183//3MrU4X24YGh8DyRDGaEvB4aa2WAzSwNuABZ0aPMicIGZpZhZN2AysDG8pYqInLhf/XMr9U0+7p8x0utSIi7oCN055zOze4FXgGRgrnNuvZnd3bZ/lnNuo5n9DVgDBIA5zrl1kSxcRCSY7fsO8bt3dnL9xCKGxcHlcYMJZcoF59xCYGGHx2Z12P5f4H/DV5qIyKn50csbSU9J4iuXxs8FuI5HZ4qKSFxasm0fr6zfy+emDqFvjwyvy+kSCnQRiTs+f4DvvbSBwl6ZfOaC+DvF/1gU6CISd+YvL2fTB3V8c8ZIMlLj8ySizijQRSSu1Da08NNXNzPl9FymjenvdTldSoEuInHlF//YQu3hFh742GjM4vckos4o0EUkbmzZW8dT7+zkhklFjDot2+tyupwCXUTignOOB15cR4+MFL522XCvy/GEAl1E4sJLa/bwbmkNX7t8OL26p3ldjicU6CIS8+qbfPzgrxsYU5DNDRPj81rnoQjpTFERkWj2f//Yyt6DTTz6yQkkx/HVFIPRCF1EYtqWvXX8ZvF2risuZHxRL6/L8ZQCXURiViDg+Oaf15KVkcJ90+P/aorBKNBFJGY9t7KC5Tv2c//0keQm6Aeh7SnQRSQmVdc38T8vb2TSoFyumVDodTlRQYEuIjHpfxZuor7Rx4NXjYnrZeVOhAJdRGLOkpJ9/Om9Cu688PSEWLgiVAp0EYkpjS1+/uvPaxnYuxtfuiQxFq4IlY5DF5GY8ovXtrKzuoFn7pycUJfGDYVG6CISM9btquXxt0q5vngA5w7J87qcqKNAF5GY4PMHuO/5NeR2T+P+GTrmvDOachGRmPDYolLW7TrIozePp2e3VK/LiUoaoYtI1Nv8QR2/fG0rV5yZz/Qz870uJ2op0EUkqvn8Ab723GqyMlL43szRXpcT1TTlIiJR7bFFpaypqOXhm8bTOyvd63KimkboIhK1tuxtnWqZcWZ/rhirqZZgFOgiEpWafQG+PH8VPTJS+N7MMV6XExM05SIiUemX/9jChj0Hmf2pCeRpqiUkGqGLSNRZuXM/j76xjWsnFHLZ6P5elxMzFOgiElUamn189Y+ryO+ZyQNXjvK6nJiiKRcRiSrf/8sGdtY0MO/OKfTI0AlEJ0IjdBGJGi+v3cO8ZeXcfdEQppze2+tyYo4CXUSiwp7aw9z3/FrGFvbkKx8d5nU5MSmkQDezaWa22cxKzOy+47SbaGZ+M7smfCWKSLzzBxz/+YfVtPgD/PKGs0lL0VjzZAT9VzOzZOBhYDowCrjRzI76pKKt3UPAK+EuUkTi26w3t/FOaTXf+fhoBud197qcmBXKn8FJQIlzrtQ51wzMB2Z20u4LwJ+AyjDWJyJxbvmOGn729y18bGw+12qx51MSSqAXAOXttivaHvuQmRUAVwGzwleaiMS7/Yea+eK89xnQK5MffuJMzLTY86kIJdA7+xd2HbZ/AXzDOec/7hOZ3WVmK8xsRVVVVYglikg8cs7x1WdXU13fzK9vGq9DFMMglOPQK4AB7bYLgd0d2hQD89v+uuYBM8zM55x7oX0j59xsYDZAcXFxxz8KIpJAHn+rlH9uquS7Hx/NmIKeXpcTF0IJ9OXAUDMbDOwCbgBuat/AOTf4yH0zexL4S8cwFxE5YmlpNQ/9bTPTx/TnlnMGel1O3Aga6M45n5ndS+vRK8nAXOfcejO7u22/5s1FJGSVBxu5d977DMztxo+vGat58zAK6dR/59xCYGGHxzoNcufcbadelojEoxZ/gHufeZ/6Rh+/v2Oy5s3DTNdyEZEu89DLm1i2o4ZfXH8Ww/v38LqcuKPTsUSkS7y4ahdzFm/nlnMG8h9nFwT/BjlhCnQRibh1u2r5+nNrmDQ4l299TJfEjRQFuohEVHV9E5/93Upyu6fxyM3jSU1W7ESK5tBFJGKafQE+//R7VNU38dzd52gpuQjTn0oRiQjnHN9esI6l22v48dVjGVuY43VJcU+BLiIRMfftHcxbVs49Fw/Rh6BdRIEuImH3+uZKfvDXDVw+uh9fvXS41+UkDAW6iITVht0Huffp9xjRP5ufX38WSUk6E7SrKNBFJGz21B7m008uJzszlbm3TaRbmo676Er61xaRsKhrbOH2J5ZT3+Tj2bvPoX/PDK9LSjgKdBE5ZUcOT9xaWc8Tt01kZH621yUlJE25iMgpCQQcX39uNW9t3ccPrzqTC4f18bqkhKVAF5FT8sOXN/LCqt187fLhXDdxQPBvkIhRoIvISXt8USmPv7WdW88ZyOenDvG6nISnQBeRk/KH5WX8YOFGrhibzwNXjtZCFVFAgS4iJ+wva3Zz3/NruWhYH35+3Vkk61jzqKBAF5ET8vqmSr48fxUTB+Yy65MTSEtRjEQLvRMiErK3S/Zx9+9XMiK/B3NuKyYzLdnrkqQdBbqIhOSdbdXc8dvlDM7rzlOfnky21gONOgp0EQlq+Y4a7vjtcgb06sbvPzOZ3O5pXpcknVCgi8hxLdtew21zl9G/ZwZP3zlZi1REMQW6iBzTkm37uHXuMvr1zGDenVPo20PXZ4lmupaLiHTqra1V3PnUCopyu/H0Z6bQp4dG5tFOI3QROcprG/Zyx29XMKh3d+bdqTCPFQp0Efk3L7y/i8/+fiUj+/dg3p1T6K0585ihKRcR+dDv3t3JAy+uY/LgXObcOpGsdEVELNG7JSI45/jVP0r4+WtbuGREXx6+eTwZqTppKNYo0EUSnD/g+NaL63hmaRmfGF/AQ1ePJTVZs7GxSIEuksAON/v50vz3eXXDXj43dQhfv3y4rpoYwxToIgmqqq6Jzzy1gjUVB/j2laO4/bzBXpckp0iBLpKAtu6t4/Ynl7OvvolHb57AtDH9vS5JwiCkiTIzm2Zmm82sxMzu62T/zWa2pu22xMzGhb9UEQmHt7ZW8YlHl9DYEuAPd52jMI8jQUfoZpYMPAxcClQAy81sgXNuQ7tm24GLnHP7zWw6MBuYHImCReTkOOd44u0dPPjXDQzt24M5txYzILeb12VJGIUy5TIJKHHOlQKY2XxgJvBhoDvnlrRr/y5QGM4iReTUNPn8PPDCev6wopxLR/Xj59efpWPM41Ao72gBUN5uu4Ljj77vAF4+laJEJHw+qG3kc0+v5P2yA9x78Rn856XDSNKScXEplEDv7J13nTY0u5jWQD//GPvvAu4CKCoqCrFEETlZS0urueeZ92ho9vPwTeO5Ymy+1yVJBIXyoWgFMKDddiGwu2MjMxsLzAFmOueqO3si59xs51yxc664T58+J1OviITAOcfji0q5ac5SsjNSefGe8xTmCSCUEfpyYKiZDQZ2ATcAN7VvYGZFwPPAp5xzW8JepYiE7EBDM//v2dW8trGSaaP78+Nrx2q5uAQRNNCdcz4zuxd4BUgG5jrn1pvZ3W37ZwEPAL2BR9rOMvM554ojV7aIdGblzv18cd77VNY18u0rR3HbuYN05mcCMec6nQ6PuOLiYrdixQpPXlsk3vgDjkdeL+EX/9hKfs8Mfn3TeM4akON1WRIBZrbyWANmHbckEuN2HTjMV+avYtmOGj4+7jQevGqMplgSlAJdJEY55/jTe7v47oL1BJzjZ9eN46qzCzTFksAU6CIxaF99E/c/v5ZXN+xl0qBcfnLtOIp666zPRKdAF4khzjkWrN7Ndxas51Czn2/OGMmnzx9Msk4UEhToIjHjg9pG/vuFtby2sZKzi3L48dVjGdqvh9dlSRRRoItEOX/A8bt3dvCTV7fgCwT47ytGcvt5GpXL0RToIlFs3a5a7v/zWtZU1HLB0Dwe/I8xDOzd3euyJEop0EWi0P5Dzfz075t5ZmkZud3T+dWNZ3Pl2HwdwSLHpUAXiSI+f4D5y8v5yaubqWv0ccs5g/jKpcPomanjyiU4BbpIlFi0pYoH/7qBLXvrmTw4l+/OHM2I/tlelyUxRIEu4rGNew7y0N828cbmKopyuzHrk+O5fHR/Ta/ICVOgi3ikvKaBn/19Cy+s2kV2RirfnDGSW84dSHpKstelSYxSoIt0sT21h3n49RL+sLycJDPuvmgId184hJ7dNE8up0aBLtJF9tQe5rE3S3lmaRkOx/UTB3DvxUPp3zPD69IkTijQRSKsrLqBR9/cxnMrywk4uHZCIfdcfAYDcnXtFQkvBbpIhKzbVctji0pZuHYPyWZcP3EAn71wiIJcIkaBLhJGgYDj9c2V/GbxdpZsqyYrPYU7zh/MHecPpl+2plYkshToImFwsLGFP62s4LdLdrCjuoH+2Rn81/QR3Di5SItNSJdRoIucgnW7anl66U5eeH83h1v8jC/K4auXDWfamP6kJid5XZ4kGAW6yAk62NjCglW7+eOKctZU1JKRmsTMcQXcPKWIsYU5XpcnCUyBLhICf8CxuGQfz79XwSvrP6CxJcCI/j34zpWjuGp8oa61IlFBgS5yDM45VlfU8tLq3by0ejeVdU1kZ6Rw9fhCrp84gDMLeur0fIkqCnSRdpxzrKmo5eV1H7Bw7R7KahpIS07iouF9uHp8AReP6KtT8yVqKdAl4bX4AyzfUcNrGyp5Zf0H7DpwmJQk45whvfnCR87gstH9NaUiMUGBLglpX30Ti7ZU8cbmKt7YXMnBRh9pKUlccEYeX/7oUC4d1Y+cbmlelylyQhTokhCafH5W7tzP4q37WFyyj7W7anEO8rLSuWx0fy4d1Y8LhubRLU2/EhK79NMrcanJ52dNRS3vbqvm3e3VrNy5n8aWAClJxtlFOXzlo8O4eHhfRp+WTZIWW5Y4oUCXuLCvvon3yw7wXtl+VuyoYXVFLc2+AAAj87O5cVIR5w3JY8qQ3mSl68de4pN+siXmHGxsYcPug6ypOMCailrWVNRSVtMAQEqSMaagJ7eeM5DiQblMGpRLr+6aC5fEoECXqBUIOMr3N7Dpgzo27alj0wcHWb/74IfhDVCQk8m4AT355JQixhf1YkxBTzJSdVihJCYFunjucLOfnTWHKK06RGlVPduqDrG1so6SynoaW1qnTcygKLcbYwqyuX7iAEadls2ZBT3Jy0r3uHqR6KFAl4jzBxyVdY3sPnCYiv2HKatuoHx/A2U1DeysbmBPbeO/tT+tZwZn9OvBzZN7M7RvFiPysxnWL0tHoIgEEdJviJlNA34JJANznHM/6rDf2vbPABqA25xz74W5VokyzjkONvrYV99EVV0TlXVNVB5spLKuiT21jXxQe7jtayO+gPu37+3TI50BvTI5Z0hvBvXuzqC87pye153T+3RXcIucpKC/OWaWDDwMXApUAMvNbIFzbkO7ZtOBoW23ycCjbV8lBjjnaGj2U9foo66xhYONLdQebrs1tHDgcAsHGlqoOdTM/oZmag613qrrm2n2B456vvSUJPplZ5DfM4Pigb3Iz8mkICeTgl6ZFOZkUtirG5lpmucWCbdQhkKTgBLnXCmAmc0HZgLtA30m8JRzzgHvmlmOmeU75/aEveI45JzDF3D4A21f/Q5fIIA/4GgJOFp8AXyBAC1+R4s/QIs/QJOvdbvZF6DZF6DJ56fJF6CpxU+jL0Bji5/Gltavh5v9NLT4Odzso6HZz6FmP4eafB/e6pt8dBhAHyU7I4Ve3dPI6ZZG/+wMRuVnk5uVRp+sdPLabv2y0+nbI4PszBRdtErEA6EEegFQ3m67gqNH3521KQDCHuhvbqniwb9sCNouSD61tnH/anVUe/fvd4+0bb1/5HGHc+22nSPgWh8PuH9tB5wjEGhtG3AOv3MEAuB3rSEeKanJRmZqMplpyWSmJtMtLYVuaclkZ6RwWs8MuqenkJWeQo+M1q9ZGSlkZ6SSnZlKj4wUcjJTyemWRnZGCilarEEk6oUS6J0NtTqmUChtMLO7gLsAioqKQnjpo2WlpzC0X1ZIba3Tso5q1Nnd1u12o0yj9UiLf923f32Ptb5WkrW2MYykpH9tJ5mRZIYZJNuRfUZyUuu2mZGa3Pp4shkpyUmt22akJSeRktz6WFqykZqc9OH+9JQk0pKTSUtJIiM1ibSUJNJTkslIbf2arDMgRRJKKIFeAQxot10I7D6JNjjnZgOzAYqLi09qaDphYC8mDJxwMt8qIhLXQvl/9HJgqJkNNrM04AZgQYc2C4BbrNUUoFbz5yIiXSvoCN055zOze4FXaD1sca5zbr2Z3d22fxawkNZDFktoPWzx9siVLCIinQnpgF/n3EJaQ7v9Y7Pa3XfAPeEtTUREToQOXRARiRMKdBGROKFAFxGJEwp0EZE4oUAXEYkT1v709y59YbMqYOdJfnsesC+M5cSKROx3IvYZErPfidhnOPF+D3TO9elsh2eBfirMbIVzrtjrOrpaIvY7EfsMidnvROwzhLffmnIREYkTCnQRkTgRq4E+2+sCPJKI/U7EPkNi9jsR+wxh7HdMzqGLiMjRYnWELiIiHUR1oJvZNDPbbGYlZnZfJ/vNzH7Vtn+NmY33os5wC6HfN7f1d42ZLTGzcV7UGU7B+tyu3UQz85vZNV1ZX6SE0m8zm2pmq8xsvZm92dU1hlsIP989zewlM1vd1ueYv3qrmc01s0ozW3eM/eHJMudcVN5ovVTvNuB0IA1YDYzq0GYG8DKtCwdNAZZ6XXcX9ftcoFfb/emx3u9Q+tyu3T9pvfLnNV7X3UXvdQ6t6/cWtW339bruLujz/cBDbff7ADVAmte1n2K/LwTGA+uOsT8sWRbNI/QPF6d2zjUDRxanbu/Dxamdc+8COWaW39WFhlnQfjvnljjn9rdtvkvrClGxLJT3GuALwJ+Ayq4sLoJC6fdNwPPOuTIA51ys9z2UPjugh7Wu85hFa6D7urbM8HLOLaK1H8cSliyL5kA/1sLTJ9om1pxon+6g9S97LAvaZzMrAK4CZhE/QnmvhwG9zOwNM1tpZrd0WXWREUqffw2MpHUZy7XAl5xzga4pzzNhybKQFrjwSNgWp44xIffJzC6mNdDPj2hFkRdKn38BfMM552+/eHeMC6XfKcAE4BIgE3jHzN51zm2JdHEREkqfLwdWAR8BhgB/N7O3nHMHI1ybl8KSZdEc6GFbnDrGhNQnMxsLzAGmO+equ6i2SAmlz8XA/LYwzwNmmJnPOfdCl1QYGaH+jO9zzh0CDpnZImAcEKuBHkqfbwd+5Fonl0vMbDswAljWNSV6IixZFs1TLom6OHXQfptZEfA88KkYHqm1F7TPzrnBzrlBzrlBwHPA52M8zCG0n/EXgQvMLMXMugGTgY1dXGc4hdLnMlr/R4KZ9QOGA6VdWmXXC0uWRe0I3SXo4tQh9vsBoDfwSNuI1edi+KJGIfY57oTSb+fcRjP7G7AGCABznHOdHvoWC0J8r78PPGlma2mdiviGcy6mr8JoZvOAqUCemVUA3wZSIbxZpjNFRUTiRDRPuYiIyAlQoIuIxAkFuohInFCgi4jECQW6iEicUKCLiMQJBbqISJxQoIuIxIn/Dxn2DQNNgga1AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "[0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13\n", " 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27\n", " 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4 0.41\n", " 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5 0.51 0.52 0.53 0.54 0.55\n", " 0.56 0.57 0.58 0.59 0.6 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69\n", " 0.7 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8 0.81 0.82 0.83\n", " 0.84 0.85 0.86 0.87 0.88 0.89 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97\n", " 0.98 0.99]\n" ] } ], "source": [ "t=np.linspace(0,1, 100, endpoint=False)\n", "def v(t):\n", " return t**4\n", "\n", "\n", "\n", "plt.plot(t, v(t))\n", "plt.show()\n", "print(t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suskaičiuokite atsutmą laiko momentu t=1 naudojantis Eulerio metodu (Ats.:$v=t^n$, $x=1/(n+1)$)\n", "$$x=x_0+v(t_0) dt$$\n", "Kiek taškų reikia vieno procento tikslumui? " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9751666650000002\n" ] } ], "source": [ "t=np.linspace(0,1, 100, endpoint=False)\n", "dt= t[1]-t[0]\n", "#print(dt)\n", "\n", "x0=0\n", "\n", "for t0 in t:\n", " x0 = x0 + v(t0)*dt\n", "\n", "ats=1/5\n", "print(x0/ats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suskaičiuokite atsutmą laiko momentu t=1 naudojantis Runge Kutta metodu\n", " $$v_1 = v(t_0,x_0) $$\n", " $$v_2 = v(t_0+ \\frac{1}{2}dt,x_0+ \\frac{dt}{2}v_1)$$\n", " $$v_3 = v(t_0+ \\frac{1}{2}dt,x_0+ \\frac{dt}{2}v_2)$$\n", " $$v_4 = v(t_0+ dt,x_0+ dt v_3)$$\n", "\n", "$$x=x_0+(v_1+2v_2+2v_3+v_4) \\frac{dt}{6}$$\n", "Kiek taškų reikia vieno procento tikslumui? " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000000004166667\n" ] } ], "source": [ "t=np.linspace(0,1, 100, endpoint=False)\n", "dt= t[1]-t[0]\n", "#print(dt)\n", "\n", "x0=0\n", "\n", "for t0 in t:\n", " v1 = v(t0)\n", " v2 = v(t0+dt/2)\n", " v3= v(t0+dt/2)\n", " v4= v(t0+dt)\n", " x0 = x0 + (v1 +2*v2+2*v3+v4)*dt/6\n", "\n", "ats=1/5\n", "print(x0/ats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Keplerio orbitos\n", "Išspręskite žemės judėjimo aplink Saulę judėjimo lygtį. Saulės masė daug didesnė nei Žemės, todėl galite laikyti, kad ji nejuda koordinačių pradžios taške. Tinkamai pasirinkite vienetus.\n", "\n", "\\begin{align}\n", " \\frac{dv^x}{dt} &= - \\frac{x}{\\sqrt{x^2+y^2}} \\frac{GM}{r^2}\\\\\n", " \\frac{dv^y}{dt} &= - \\frac{y}{\\sqrt{x^2+y^2}} \\frac{GM}{r^2}\\\\\n", " \\frac{dx}{dt} &= v^x\\\\\n", " \\frac{dy}{dt} &= v^y\n", ".\\end{align}\n", "\n", "- Išbandykite įvairias pradines sąlygas apskritiminei, eliptinei ir hiperbolinėms trajektorijoms.\n", "- Nubrėžkite orbitų grafikus.\n", "- Patikrinkite Eulerio ir RK4 metodų konvergavimą apskritiminei orbitai\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAD4CAYAAAAKL5jcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAn+klEQVR4nO3deXxU5b3H8c8zM5ns+0o2SCAEwhYgsgqCC4JV0eveaqm2Bb3S9tb2Vm9bt+6t3mu1ihata11qXdhEURQUlX0nrCFAVrIQsu8zz/0jY0oxgcCc5GRmfu/Xa14zkzmc8+OQ+fKc55zzPEprjRBCAFjMLkAI0X9IIAghOkkgCCE6SSAIITpJIAghOtnMLuBMYmJi9KBBg8wuQwivsnXr1kqtdWxXn/XrQBg0aBBbtmwxuwwhvIpS6lh3n8khgxCikwSCEKKTBIIQopMEghCikwSCEKKTBIIQopMhgaCUel4pVa6U2tPN50op9YRSKk8ptUspNc6I7QohjGXUdQgvAk8CL3fz+Rwgw/WYCDztehb9XGu7k4r6Fsprm6ltbqehpZ36lo7nxlYHAKfeQm+xKILtNoLsVoL9O54jguzEhfoTE+KP3SaN0v7MkEDQWn+mlBp0hkXmAi/rjt+cDUqpCKXUAK11qRHbF+dPa015XQv5FQ3kV9ZzpKKBI5UNFFc3UV7XQlVDq6HbiwzyIy40gJSoINJjg0mLCWZQdDCDY4OJDfVHKWXo9sS56asrFZOAwlPeF7l+9rVAUErNB+YDpKam9klxvqSkuoldRTXsLq5md3Etu4uqOdnY1vm5v81CWkwwyZGBjB8YSVxoAPFh/sSF+RMW4Eewv40Q1yPQbsXi+gJ/9T12ODWNrY7OFkR9SzvVja2U17VQXttCeV0zZbUtFFQ18NmhClrbnZ3bjgnxZ2RSGKOSwhmZFE52SgTxYQF9un98XV8FQlex3+VQTVrrxcBigJycHBnOyU0l1U2sP3yCDfknWJ9/gqKTTQBYLYqh8aFclhVP1oAwBseFkBYTTGJ4IBbL+f8v7WeFAD8rUcH2sy7rcGpKqps4eqKBvPJ69hTXsqe4hs8OVuB0/cunxQQzKT2KSenRTEyLJiFcAqI39VUgFAEpp7xPBkr6aNs+pd3hZFtBNav3lbF6bxn5lQ0ARAT5MSktmu9emEZ2SgTDB4QR4Gc1tVarRZESFURKVBDTMv51r01Tq4O9pbVsLzjJhvwqVuwq5fVNHQ3MofEhXDo8nkuz4slOjnArvMTXKaPGVHT1IazQWo/s4rNvAAuBK+joTHxCaz3hbOvMycnRcnPT2bU5nHx+qJLlu0pYs7+ck41t+FkVkwfHcNHQWCanRzMsIdRjvzwOp2ZfaS3rD5/gk/3lbDpahcOpiQmxc8mweK7OTmRSejRWD/379TWl1FatdU6XnxkRCEqp14EZQAxQBjwI+AForZ9RHT1FTwKzgUbgdq31Wb/pEgjd01qzs6iGJduLWb6zhBMNrYQF2LhkeDyXDo9n+tAYQgP8zC6zV9Q0trH2YDmr95WzZn859S3txIf5c/WYRK4Zm0TWgDDpnDyDXg+E3iKB8HW1zW28vbWIv284xuGKBuw2C5cNj2dudiIzMuN87rRec5uD1fvKWLK9hLUHyml3arIGhHHb5IHMzU4kyN6v7/A3hQSCF9hXWsvL64+xZHsxTW0OslMiuGVCCnNGDSDMS1sC56qqoZUVu0p4bWMB+4/XERpg4/rxydw2aSDpsSFml9dvSCB4KK01G/KrWLQ2j3WHKvG3WZibnchtkwYxKjnc7PL6La01W46d5JX1x3h/TyntTs3lWQn858zBjE6OMLs8050pEKQ91Q85nZqP95ezaG0e2wuqiQnx597Zw7hlQgoRQWc/nefrlFJcMCiKCwZFUV43nFfWH+OlL4/yQe5xpmXEcNeMwUxOj5Z+hi5IC6GfWXeogj9+sJ89xbUkRway4KLB3DA+2fRThJ6urrmNVzcW8Ny6I1TWtzA5PZp75wwjOyXC7NL6nBwyeIBdRdX88YP9fJF3gqSIQH582VCuyU7EZvWtTsLe1tzm4PVNBTz5SR4nGlq5YlQCP52V6VN9DBII/djxmmZ+t3Ify3aWEBVsZ+HMIXxrUir+NmkR9Kb6lnae/SyfZ9fl09Lu5NaJqdwzK5PwQO/voJVA6Ida25288MURHv/4EA6nZv70dOZPT/faawf6q4q6Fh7/+CCvbSwgKtjO/8wZzn+MS/Lq/gUJhH5m/eET3L90D3nl9Vw6PI4HrhxBanSQ2WX5tD3FNfxyyR52FFZzwaBIfnPNKDITQs0uq1dIIPQT9S3t/H7lPl7dWEBqVBAPXpXFJcPjzS5LuDidmn9uLeQP7++nvqWd/7p0KAump3tdP46cduwHvsyr5Gdv76K4uonvXZjGTy/PlDMH/YzForjpglQuHR7PA0tzeWTVAT7MPc6jN4whI947Wwun867o64ea2xw8tCyXbz63ET+rhX8umMwvr8ySMOjHokP8eepb43jym2MpqGrkG3/5nOfW5dOfW9NGkRZCL8qvqGfha9vZW1rLd6YM4t7Zwwi0SxB4iitHJzIxLZqfv7ub37y3jw35VTx6w2ivvjhMWgi9ZOmOYq76y+eU1jTx/HdyeOjqERIGHig21J/Ft43nwauy+PRgOd944nO2FZw0u6xeI4FgsNZ2Jz9/dzc/emMHIxLDWfmjaVw8TDoOPZlSitunpvHWnVOwWODGZ9bz4hdHvPIQQgLBQJX1Ldz63EZe21jAXTMG89r3JzIgPNDssoRBxqREsOIH05iRGcdDy/fy83f3/NuYkN5A+hAMkltSw/yXt1JZ38ITt4zl6jGJZpckekF4oB+LbxvPox8eYNHaw+RX1PP0reN7NIakJ5AWggFW7y3j+qfX43Bq3rpzioSBl7NYFD+bPYw/35TN9sJq5j71OfkV9WaXZQgJBDe9ubmQBX/fypC4EJYtnCrjFPiQa8Ym8Y/5k2hscXD9M+vZXVRjdkluk0A4T1prnlqTx8/e3sWUwdG8MX8ScTKHgM8ZmxrJP++cTKCflZsXr+fLvEqzS3KLBMJ50Frz6xX7eGTVAa7JTuRv8y4g2F+6Y3xVemwIb981haTIQL7zwmY+2HPc7JLOmwTCOdJa89CyXJ7/4gh3TE3j/27M9rmBTcXXJYQH8OaCyYxMCmPha9s8NhTkN/kcaK15cFkuL60/xoLp6dx/5XCPnetAGC8iyM7L353I6ORwFr62jQ9zPS8UJBB6SGvNA0tzeXn9MRZclM59c4Z59T3z4vyE+Nt48Y4JjEgK5+7XtrF6b5nZJZ0TCYQe+tOqA7yywRUGsyUMRPfCAvx4+Y4JDB8Qxn++uo31h0+YXVKPSSD0wHPr8nl67WFunZQqYSB6JDywIxQGRgcx/5Ut7D9ea3ZJPSKBcBbvbi/iN+/t44pRCTx89UgJA9FjEUF2XrxjAsF2G/Oe30RxdZPZJZ2VBMIZrDtUwX//s+M6g8duypbJRMU5S4oI5MU7LqCx1cG85zdR09hmdklnJIHQjfyKeu5+dRtD4kL4623jZRRkcd6GJYTx7LdzOHaigR++sR2Hs//eJSmB0IWaxja+99IWbFYLz347R0ZCFm6blB7Nw1eP5NODFfxp1X6zy+mWXF53mnaHk4Wvb6PwZCOvfm8SKVEyGrIwxjcnppJbUsNfP80na0AYc7OTzC7pa6SFcJr//egg6w5V8ptrRjIhLcrscoSXefCqEVwwKJJ7397FvtL+d+ZBAuEUaw6U8/Taw9wyIZWbLkg1uxzhhew2C4u+NZ7QAD9+8Pp2mlodZpf0byQQXEprmvjJmzsZlhDKg1dlmV2O8GKxof48dmM2hyvq+dWKXLPL+TcSCHT0G/zw9e00tzl46lvjZIh00esuzIjhzosG8/qmQt7bVWp2OZ0kEIC/fpbP5qMn+d21oxjsQ7MAC3Pdc9lQslMiuO+dXZTW9I+Llnw+EPaV1vLn1Qf5xugBXDO2//X6Cu/lZ7Xw+M3ZtDmc/PLdPf1iFGefDoTWdic/eXMn4YF+/HruSLPLET5oYHQwP52Vycf7y1m6o8Tscnw7EJ5ck8fe0lp+d+0orxk1V3ie26emMTY1goeW51JR12JqLT4bCHnldSxak8c12YnMGpFgdjnCh1ktikeuH01ji4Nfrdhrai2GBIJSarZS6oBSKk8pdV8Xn89QStUopXa4Hg8Ysd3zpbXm/iW5BNmt/PJKOcUozDckLpQ7Zwxm+c4SNuabN36C24GglLICTwFzgCzgFqVUV9+ydVrrbNfjV+5u1x3LdpawPv8E/z17GDEh/maWIkSnuy4aTGJ4AA8t32vaDVBGtBAmAHla63ytdSvwBjDXgPX2irrmNn7z3j5GJ4fzzQlyNaLoPwLtVn7xjSz2ldby+qYCU2owIhCSgMJT3he5fna6yUqpnUqp95VSI7pbmVJqvlJqi1JqS0VFhQHl/btFaw9TUdfCr+eOlPENRL9zxagEJqVH8eiHB6hp6vuxE4wIhK6+Vae3d7YBA7XWY4C/AEu6W5nWerHWOkdrnRMbG2tAef9SWtPE858f4ZrsRMakRBi6biGMoJTi/iuzqG5s49nP8vt8+0YEQhGQcsr7ZODfTqhqrWu11vWu1ysBP6VUjAHbPiePfXQQreEnszL7etNC9NiIxHCuHD2A5784QmV9356GNCIQNgMZSqk0pZQduBlYduoCSqkE5RqMUCk1wbXdPu1KPXC8jre2FnHb5IEyxoHo93582VCa2xwsWnO4T7frdiBorduBhcAqYB/wptY6Vyl1p1LqTtdi1wN7lFI7gSeAm3UfX6f52EcHCbbbWDhzSF9uVojzMjg2hOvGJfP3jcf69D4HQ65D0Fqv1FoP1VoP1lr/1vWzZ7TWz7heP6m1HqG1HqO1nqS1/tKI7fbUwbI6Psg9zu1TBxEpVyQKD/HDSzJwODXPrTvSZ9v0iSsVF63JI8hu5fapaWaXIkSPpUQFcdXoAbyxqaDPzjh4fSAcrWxg2c4Sbp00UFoHwuPMnz6YhlYHr2481ifb8/pAWLwuH5vVwvculNaB8DxZiWFMy4jhhS+O0tLe+8OteXUg1DS28e62Yq7NTiIuLMDscoQ4L/Onp1NR18KKnb0/spJXB8KbWwppanMwb8ogs0sR4rxdOCSGtJjgPrmc2WsDweHUvLzhKBMGRZGVGGZ2OUKcN6UUt0xIYcuxkxwsq+vVbXltIKw9UE5hVZO0DoRXuG5cMnarpddbCV4bCP/YXEhMiD+zRsSbXYoQbosO8efykQm8s62Y5rbe61z0ykCoamhlzYFyrslOxM/qlX9F4YNuGJ9MTVMbnx40/i7gr3jlt2X5zhLaHJrrxiebXYoQhpkyOJqoYDsrenEeB68MhHe2FZE1IIzhA6QzUXgPm9XC7JEJrN5bRmNre69sw+sC4UhlAzuLaviPcTLHgvA+V41OpKnNwSf7y3tl/V4XCKtyjwMwe6SMpCy8z4S0KGJD/Xl/9/FeWb9XBsLIpDCSI2XMA+F9rBbFzMxYPjtUQbvDafj6vSoQymqb2V5QzeVZ0joQ3mtGZhx1ze1sK6g2fN1eFQgf7S0D4HI5XBBe7MKMGKwWxdoDxvcjeFUgfHawgqSIQDLiZAZn4b3CAvwYPzCSNQeMvx7BawLB4dSszz/BtIwYXMM3CuG1pmfEsK+0lurGVkPX6zWBsLu4hrrmdqYO6fPBnIXoczmDogDYeuykoev1mkD4Iq8S6LiaSwhvNyY5Aj+rYosEQtc25J9gWEIo0TJXo/ABgXYrI5PC2XK0ytD1ekUgOJ2aHQXVjB8YaXYpQvSZnIGR7CysobXduOsRvCIQ8isbqGtpl+nZhE8ZmRROq8PJ4Yp6w9bpFYGwo7AagLESCMKHZLlu3tt/vNawdXpJIJwk1N/G4Fi5/kD4jrSYYOxWC/tKjRtWzSsCYV9pHcMTw7DI9O7Ch9isFjLiQ9hXKi2ETlprDpbVMTReWgfC92QmhBo68KrHB0J5XQt1ze0MjQ81uxQh+tzAqGDKalsMG2fR4wPhUFlHD+sQuX9B+KCUqEAAiquNmSHa4wMhr7yjuSSBIHxRSlTHuB+FVY2GrM/jA6HoZBMBfhZi5QpF4YNSXAMBFZ6UFgIAJTVNJEYEyh2OwifFhfpjtSiO10ggAFBc3UxieKDZZQhhCotFERnkR1VDmzHrM2QtJiqtbiIxQmZ2Fr4rMsjOyQZjxkXw6EBwODUV9S0kyFTvwodFBtupMmigFI8OhNqmNrSGiCC72aUIYZqoIDtV0kKAmqaO46aIID+TKxHCPMH+Nppa5cIkql2BEB4ogSB8l91moaW9HwWCUmq2UuqAUipPKXVfF58rpdQTrs93KaXGGbHdGgkEIfC3WWgxaJAUtwNBKWUFngLmAFnALUqprNMWmwNkuB7zgafd3S5AY0vHhJdBdpsRqxPCI/nbLIaNmmREC2ECkKe1ztdatwJvAHNPW2Yu8LLusAGIUEoNcHfDbU4NgN0mFyUJ32W3WWh1ONFau70uI/5rTQIKT3lfBEzswTJJwNcmuldKzaejFUFqauoZN/zV3HY2i0d3hQjhlhmZcUQF29Ea3L1g14hA6KqE06OqJ8t0/FDrxcBigJycnDNGXruj42M/mwSC8F3jB0YaNsCwEd+kIiDllPfJQMl5LHPO2pxftRDkkEEIIxgRCJuBDKVUmlLKDtwMLDttmWXAt11nGyYBNVrrrx0unCurq33kcLp/7CSEMOCQQWvdrpRaCKwCrMDzWutcpdSdrs+fAVYCVwB5QCNwu7vbhY7OFMDQcemF8GWGnK/TWq+k40t/6s+eOeW1Bu42Ylun6gwEhwSCEEbw6N44u1VaCEIYyaMDIcDPCkCTQQNMCuHrPDoQvrqpqbrRmMEhhPB1Hh0Ika7bnk8adC+4EL7OowPhXy0ECQQhjODRgRDib8NmUZyUQwYhDOHRgaCUIjbUn/LaFrNLEcIreHQgACRFBFJcbcwkFUL4Os8PhMhAw6axEsLXeX4gRARSWt0s9zMIYQCPD4TkyCDanZrjtc1mlyKEx/P4QEiPDQYgr7ze5EqE8HweHwhD40MBOHi8zuRKhPB8Hh8IUcF2YkL8OVgmgSCEuzw+EAAyE0IkEIQwgHcEQnwYB8rqaJNxEYRwi1cEQnZqBM1tTg5IP4IQbvGKQBiXGgHA9oKT5hYihIfzikBIiggkNtSf7QXVZpcihEfzikBQSjE2JYJt0kIQwi1eEQgAE9KiOHqikRK5r0GI8+Y1gXBhRgwAnx+qNLkSITyX1wRCZnwoMSH+rMuTQBDifHlNICiluHBINF/mVeKUOx+FOC9eEwgA0zJiOdHQyu7iGrNLEcIjeVUgXDwsDqtFsSr3uNmlCOGRvCoQIoPtTE6P5oM9x+mYPU4IcS68KhAALh+ZQH5lg4yPIHzO3z4/wiE3b/LzvkDIikcpWLlbDhuE7zhS2cCvV+xl9b5yt9bjdYEQFxbAhEFRLNlRLIcNwme8tbUQi4L/GJfk1nq8LhAArhufzJHKBrmUWfgEh1PzzrZipg+NJT4swK11eWUgXDFqAIF+Vt7aWmx2KUL0ui/yKimtaeaG8Slur8srAyHE38ackQms2FVCs0wVL7zc3zccIzLIj0uGx7m9Lq8MBIDrc5Kpa25n+c4Ss0sRotcUnWxk9b4ybpmQSoCf1e31eW0gTE6PZmh8CC9+eVQ6F4XXemXDMZRS3DppoCHr89pAUEoxb8ogcktq2XpMOheF92luc/CPzYXMyoonMSLQkHV6bSAAXDs2ibAAGy98edTsUoQw3D+3FFLd2Ma8KYMMW6dXB0KQ3cZNF6TwwZ7jFFbJDNHCe7Q5nDzzaT7jUiOYmBZl2Hq9OhAA7rgwDatSPPPpYbNLEcIwS7YXU1zdxMKLh6CUMmy9bgWCUipKKfWRUuqQ6zmym+WOKqV2K6V2KKW2uLPNczUgPJDrc5L555YiymRCWOEFHE7N02sPkzUgjJmZ7p9qPJW7LYT7gI+11hnAx6733Zmptc7WWue4uc1zdtdFg3FozeLP8vt600IYbtnOYvIrG7h7prGtA3A/EOYCL7levwRc4+b6ekVKVBBzsxN5deMxyqWVIDxYS7uD//3wICMSw5gzMsHw9bsbCPFa61IA13N37RcNfKiU2qqUmn+mFSql5iultiiltlRUVLhZ3r/86JIMHE7Nnz8+ZNg6hehrr24ooOhkE/fOHobFYmzrAHoQCEqp1UqpPV085p7DdqZqrccBc4C7lVLTu1tQa71Ya52jtc6JjY09h02c2cDoYL41cSD/2FzI4QoZK0F4nrrmNp5ck8fUIdFMc40ybrSzBoLW+lKt9cguHkuBMqXUAADXc5c3Y2utS1zP5cC7wATj/go9t/DiIQTYLPzpg/1mbF4Ityxae5iqhlbunT3M8L6Dr7h7yLAMmOd6PQ9YevoCSqlgpVToV6+BWcAeN7d7XmJC/Flw0WBW5Zax6UiVGSUIcV4OV9Tz3Lp8rhuXzOjkiF7bjruB8AfgMqXUIeAy13uUUolKqZWuZeKBz5VSO4FNwHta6w/c3O55+960NJIiAnlg6R7aZfp44QG01jy0LJcAPyv3zRnWq9tyKxC01ie01pdorTNcz1Wun5dora9wvc7XWo9xPUZorX9rROHnK8hu4/4rs9h/vI4X5ZJm4QFW5R5n3aFK7rlsKLGh/r26La+/UrErl4+IZ0ZmLH9efUguVhL9Wm1zGw8v38uwhFBuM+iOxjPxyUBQSvHw1SNodTh5eHmu2eUI0a3fvbePstpm/njdaGzW3v+6+mQgQMdpyB9dksHK3cdZsUsGURH9z7pDFbyxuZDvT09nTEpEn2zTZwMBYMH0dMYkh3P/kj1U1LWYXY4Qnepb2rnv7d2kxwbz40uH9tl2fToQbFYLj94whoYWB79csltGVhL9xoNLcymtaeKR60cbMjRaT/l0IABkxIdyz6yhrMot462tRWaXIwRLthfz9rYifnBxBuMHGjfWQU/4fCAAfH9aOpPSo3hgaS555e5NhSWEO46daOCXS/ZwwaBIfnDxkD7fvgQCYLUoHr95LIF2Kwtf2y5DtwtTtLQ7+OHr27Eo+PPNY/vkrMLpJBBc4sMC+N8bx7D/eB2/WrHX7HKED3po2V52FtXwp+vHkGTQoKnnSgLhFDMz41hwUTqvbSzgzc2FZpcjfMhrGwt4fVMBd88czOxeGOegpyQQTvPfszKZlhHDL5bsZusxuQFK9L6tx07y4LI9zMiM5Z7LMk2tRQLhNDarhb/cMpbEiEAWvLKN0poms0sSXqywqpE7/76VxIhAHr9pLNZeGPTkXEggdCEiyM5z386hqbWd77+8hYaWdrNLEl6oprGN21/cTEubg7/NyyE8yM/skiQQupMRH8pfvjmWvSW13PXqNtrkVmlhoJZ2B/Nf2ULBiUYWfzuHIXGhZpcESCCc0cXD4vndtaP47GAF9769S65kFIZwODU//ecuNh6p4pEbRjMpPdrskjrZzC6gv7t5QipltS08tvog8WEB3Du7dweoEN5Na80v3t3N8p0l3DdnGHOzk8wu6d9IIPTADy8ZQlldM0+vPUyw3crCizPMLkl4IK01Dy/fyxubC1k4cwh3XjTY7JK+RgKhB5RS/HruSJpbHTz64UGsFgt3zeh//5ii/9Ja88iqA7z45VHumJrGT2b13R2M50ICoYesFsUjN4yh3an54wf7sVkU35+ebnZZwgNorfn9+/tZ/Fk+t0xI5f4rh/faqMnukkA4B1aL4v9uHIPDqfntyn20Opz854zB/fYfV5jP6dQ8sGwPf99QwG2TBvLw1SP69e+LBMI5slkt/PnmbGxWxSOrDlDT1Mb/zOm9cfKF52p3OPnZ27t4Z1sxCy5K575enE/BKBII58HPauGxG7MJD/Rj8Wf5VDe28rtrR5lyd5ronxpa2vnB69v5ZH8591w2lB8YPG17b5FAOE8WS8dArRGBfjzxSR5VDW08fnM2wf6yS31deW0zd7y0mb0ltfx67ghumzzI7JJ6TP5Lc4NSintmZfLw1SP4ZH8ZNzyzXu598HGHyuq4dtGX5Fc08Ny8HI8KA5BAMMS8KYP423cuoKCqkblPfsGuomqzSxIm+DD3ONcu+pJWh5N/zJ/MxcPizS7pnEkgGGRmZhxv3zUFu83CjX9dzzvbZHxGX+F0av7vwwPMf2Ur6bHBLLl7KqOSw80u67xIIBgoMyGUJXdPJTslgnve3MnP390tw7F5uZrGNr770mae+CSPG8Yn8+aCyaaNdmQE6QEzWEyIP3//7kQe/fAgz3x6mN1FNSz61jhSooLMLk0YbPPRKv7rjR2U1Tbz62tGcuvEVI84k3Am0kLoBTarhfvmDOPZb+dw9EQDVzy+jne3F8ndkl6i3eHksY8OctNf12OzKt66awq3TRro8WEAEgi96rKseFb+cBqZCaH8+B87WfjadqobW80uS7jh2IkGbnl2A49/fIhrspN474fTyO6jadb6ghwy9LKUqCD+sWAyz3x6mMc+OsiWY1X84brRzMyMM7s0cQ4cTs0LXxzh0Q8P4Gex8NhNY7h2bLLZZRlOWgh9wGpR3D1zCEvunkpYgB+3v7CZH7y+XeaT9BB55XVc/8yX/Oa9fUwdHMNH91zklWEAoPrzcW1OTo7esmWL2WUYqqXdwdNrD7NozWEC/Cz8zxXDuSknBYvJg2uKr2toaefJNXk8ty6fYH8bD189gqvHJHp8X4FSaqvWOqfLzyQQzJFXXs/P393NpiNVZKdE8MBVWYxLjTS7LEHH7crLd5Xyu/f2cby2mevGJXPfnGHEhvqbXZohJBD6KadT8/a2Iv606gAVdS1cPSaRe+cM8+jz2J5ue8FJfv/+fjYdqWJkUhgPXz2S8QO9K6glEPq5+pZ2nll7mGfX5QNw+9Q05k9PJyrYbnJlviOvvI5HVh1gVW4ZMSF2/uvSodwyIdX0eRJ6gwSChyg62cgjqw6wbGcJQX5Wbp+axvempRERJMHQW45UNrBoTR5vbysiyG5j/vR0vnthmlfftSqB4GEOltXx+OpDvLe7lFB/G/OmDGLelEFecwzbH+wtqWXR2jxW7i7Fz2rhWxMHcvfMwUSHeP8+7rVAUErdADwEDAcmaK27/PYqpWYDjwNW4Dmt9R96sn5fDYSv7Cut5fHVh1i19zh+VgvXZifxvWlpZMT3j0k9PI3Tqfk8r5IXvjjCmgMVhPjbuG3yQO6YmuZTYdubgTAccAJ/BX7aVSAopazAQeAyoAjYDNyitT7rnOu+Hghfya+o52+fH+GtrUW0tDu5aGgs35yYysXD4vCTUZrOqra5jbe3FvHK+mPkVzYQE2LnO1MGcdvkQYQHmj99Wl87UyC4daCktd7n2sCZFpsA5Gmt813LvgHMBc4aCKJDemwIv712FD+Zlckr64/x2qZjLHhlK3Gh/tyYk8JNF6TIzVOncTo1G49U8fa2IlbuLqWx1cG41Agevzmb2SMT8LdZzS6xX+qLnpMkoPCU90XAxO4WVkrNB+YDpKam9m5lHiYq2M6PLs3g7pmDWXOggjc2FbBobR5PrskjZ2AkV2cnMmfkAJ9q/p4ur7yOpTtKeGdbMcXVTYT427hqdCK3ThrosWMU9KWzBoJSajWQ0MVHv9BaL+3BNrpqPnR7nKK1Xgwsho5Dhh6s3+fYrBYuy4rnsqx4SqqbeGdbESt2lfLA0lweWpbLlMExXD4inpnD4kiO9O6Wg9aa3cU1fLDnOKtyj3O4ogGLggszYvnZ7ExmZSUQaJfWQE+dNRC01pe6uY0iIOWU98lAiZvrFC6JEYEsvDiDhRdncLCsjuU7S1i+s4T7l+bC0lwy40OZOSyOGZmxZKdEEODn+V+OyvoWvsir5PNDlXyeV0lpTTNWi2JSehTfmTKIWSMSiA8LMLtMj2TIaUel1Fq671S00dGpeAlQTEen4je11rlnW690Kp4frTWHKxpYs7+cT/aXs/loFe1Ojd1mITslgklpUUxIi2Z0SjhhAf27U01rTdHJJrYVnGR7QTWbjlSxt7QWgPBAP6YOiebiYfFcOjxOrtfood48y3At8BcgFqgGdmitL1dKJdJxevEK13JXAH+m47Tj81rr3/Zk/RIIxqhtbmPD4RNsOlLFpqNV7Cmuwen6Z0+LCWZEYhijksIZPiCM9NhgEsMDTbnZqrXdSX5lPQeO13GwrI4Dx+vYUVhDZX3HXaGBflayUyK4MCOGC4fEMDIp3CuvJOxtcmGS+Dd1zW1sPXaSPcU17C6uYU9xLcXV/xo+3t9mIS0mmLSYYJIiAokPCyAuzJ+EsABiQ/0JDfAjNMCGv83S4zv/mlodVDW2crKhlRMNrVQ1tFBS3UxhVSNFJ5soPNlI8ckm2l1JZbMo0mKCGZUcztjUSMalRpAZHyqT4RhAAkGcVVVDK/uP13K0spEjlfUcqWwgv7KB0upmmroZKNbPqggN8MPfZsGiFBYLWJXCohQt7U6a2xwdj3YnDmfXv2cxIXaSI4NIjgxkYHQQQ+NDyUwIJT0mBLtNvvy9odeuQxDeIyrYzpTBMUw5bZZ7rTV1Le2U1zZTVttCRV0Ldc1t1Da3U9/STl1zGy1tTpwanFrjcGqcWuNvsxLgZ+l8Dva3ER1sJzLYTnSwnahgOwPCA+UMQD8jgSDOSClFWIAfYQF+DImTS6a9nbTJhBCdJBCEEJ0kEIQQnSQQhBCdJBCEEJ0kEIQQnSQQhBCdJBCEEJ369aXLSqkK4Ng5/rEYoLIXyukrnly/J9cOnl3/udQ+UGsd29UH/ToQzodSakt312l7Ak+u35NrB8+u36ja5ZBBCNFJAkEI0ckbA2Gx2QW4yZPr9+TawbPrN6R2r+tDEEKcP29sIQghzpMEghCik8cHglLqBqVUrlLKqZTq9rSLUmq2UuqAUipPKXVfX9Z4JkqpKKXUR0qpQ67nyG6WO6qU2q2U2qGUMnVcubPtS9XhCdfnu5RS48yosys9qH2GUqrGtZ93KKUeMKPOriilnldKlSul9nTzufv7XWvt0Q86JprNBNYCOd0sYwUOA+mAHdgJZJldu6u2PwH3uV7fB/yxm+WOAjH9oN6z7kvgCuB9OibpmQRsNLvuc6h9BrDC7Fq7qX86MA7Y083nbu93j28haK33aa0PnGWxzvkltdatwFfzS/YHc4GXXK9fAq4xr5Qe6cm+nAu8rDtsACKUUgP6utAu9Offg7PSWn8GVJ1hEbf3u8cHQg91Nb9kkkm1nC5ea10K4HqO62Y5DXyolNrqmv/SLD3Zl/11f/e0rslKqZ1KqfeVUiP6pjRDuL3fPWKQ1b6eX9JoZ6r/HFYzVWtdopSKAz5SSu13/Y/R13qyL03d32fQk7q20XGtf71rgqElQEZvF2YQt/e7RwSC9vD5Jc9Uv1KqTCk1QGtd6mrelXezjhLXc7lS6l06mr9mBEJP9mV/nc/zrHVprWtPeb1SKbVIKRWjtfaEm57c3u++csiwGchQSqUppezAzcAyk2v6yjJgnuv1POBrLR6lVLBSKvSr18AsoMue5j7Qk325DPi2q9d7ElDz1WGRyc5au1IqQbmmo1JKTaDjO3Kizys9P+7vd7N7Tg3oeb2WjmRsAcqAVa6fJwIrT+uBPUhHL/MvzK77lLqigY+BQ67nqNPrp6NXfKfrkWt2/V3tS+BO4E7XawU85fp8N92c/emntS907eOdwAZgitk1n1L760Ap0Ob6nf+u0ftdLl0WQnTylUMGIUQPSCAIITpJIAghOkkgCCE6SSAIITpJIAghOkkgCCE6/T9BxIx9iTVWcQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Euler method\n", "\n", "n=1000\n", "t=np.linspace(0,2*np.pi, n, endpoint=False)\n", "dt= t[1]-t[0]\n", "\n", "xcord = np.zeros(n)\n", "ycord = np.zeros(n)\n", "tcord = np.zeros(n)\n", "#print(dt)\n", "\n", "x0=1\n", "y0=0\n", "vx0=0\n", "vy0=1 # v = dx/dt = r0 /(r/vI) = vI\n", "\n", "def pag(x,y, t):\n", " r=np.sqrt(x**2+y**2)\n", " ax = -x/r**3\n", " ay = -y/r**3\n", " return ax,ay\n", "\n", "i=0\n", "for t0 in t:\n", " \n", " ax0,ay0 = pag(x0,y0, t0)\n", " x0 = x0 + vx0*dt\n", " y0 = y0 + vy0*dt\n", " vx0 = vx0 + ax0*dt\n", " vy0 = vy0 + ay0*dt\n", " xcord[i]=x0\n", " ycord[i]=y0\n", " tcord[i]=t0\n", " i=i+1\n", " #v1 = v(t0)\n", " #v2 = v(t0+dt/2)\n", " #v3= v(t0+dt/2)\n", " #v4= v(t0+dt)\n", " #x0 = x0 + (v1 +2*v2+2*v3+v4)*dt/6\n", " \n", "plt.plot(xcord,ycord)\n", "ax=plt.gca()\n", "ax.set_aspect(\"equal\")\n", "plt.show()\n", "\n", " " ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAD4CAYAAAD2OrMWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsGElEQVR4nO3deXxU9b3/8dcn+76HJARCQti3sAkICuKO1iKtWsWl1qq1xettbXvVX+/1tle7qW3V6tWKetVWrdaKu4K4YBEBCUtYAgmEACEh+75n5vv7YyY0xYQsM5n183w85sFk5iyfmZD3nHPmnO9HjDEopdRABbi7AKWUd9HQUEoNioaGUmpQNDSUUoOioaGUGpQgdxcwFElJSSYzM9PdZSjls3Jzc6uMMcm9PeeVoZGZmcm2bdvcXYZSPktEjvT1nO6eKKUGRUNDKTUoGhpKqUHR0FBKDYqGhlJqUJwSGiLyrIhUiMiePp4XEXlURA6KSJ6IzO7x3MUicsD+3N3OqEcpNXyctaXxHHDxaZ5fBoy3324FngAQkUDgcfvzU4BrRGSKk2pSSg0Dp5ynYYz5TEQyTzPJcuAFY7sOf7OIxIlIGpAJHDTGFAGIyF/t0+5zRl1qeFmthtqWDupaO6lv7aS+pZO61g6a2rrosBgsViudFkOXxWA1htDgAEKDAgkNCiA0KICIkCDiI4NJjAwlITKE+IhgggJ1j9nTuerkrnTgWI+fS+yP9fb4/N4WICK3YttKISMjY3iqVF/RZbFSXN1MQXkTByuaKKlt4XhdK8drWymtb6Ojy+q0dYlAYmQIoxMiGB0fQUaC7ZaVHMnE1GhiwoKdti41dK4KDenlMXOax7/6oDFPAU8BzJ07V0cOGgadFiv5ZQ3kHqll17E69p9opKiymQ7LP4NhRHQoI+PCmZoey4VTU0mNCSMhMoTYiGDiwoOJDQ8mOiyYkMAAggKFwAAhODAAATosVto7rbR3WWjvstLU3kVtcwfVzR3U2P+taGjjWG0LO47V8u7uMizWf/6qR8WHMzkthslpMcxIj2VuZjxxESFueKf8m6tCowQY3ePnUUApENLH48oF2rss5B6pZWNhFduKa8k7Xkdbpy0gUmJCmZIWw5KJyUxMiWZCSjTjRkQRFhw45PWFBQTa5x/YFkOnxUpZXRsHKxvJL2tkX1kD+WUNrM8vp3vAuYkp0czNjGdeVgJnZicyIjpsyPWpgXFVaLwF3G4/ZjEfqDfGlIlIJTBeRLKA48DVwEoX1eSXjla38PH+cj4rrGJzUTUtHRaCAoSp6bFcMy+DOWPimZ0Rz8i4cHeXSnBgABmJEWQkRnDupJSTj7d2WMgrqePL4hq2Ftfy5s5SXtxyFIDp6bEsnZjMOZNGkDMqjsCA3jZmlSPEGWOEisjLwDlAElAO/Df2jxNjzJMiIsBj2L5haQG+Y4zZZp/3EuBhIBB41hjzy/7WN3fuXKMXrA1ccVUz7+4u473dZewtbQAgMzGCs8cns3hCMmdmJxIV6pXXLgJgsRryyxrYUFDJJ/sr2H60FquBhMgQLpqawmU5I5mflagBMggikmuMmdvrc944sLCGRv+qmtpZs/04a3YcZ1+ZLShmZcRx6fQ0LpiSwpjESDdXOHzqWjr4rLCK9fvKWZ9fTkuHhRHRoVw6I43lM9PJGRWL7XNM9UVDw09YrIbPCit5Zesx1ueX02U1zBwdx9dmpLFsehrpHrDL4WqtHRY+2l/O27tK+WR/JR0WK5NSo7l2fgbLZ6XrNzJ90NDwcfUtnby09Sh//qKY0vo2EiND+MbsdL51xmjGjYh2d3keo6Gtk3d2lfHiliPsLW0gPDiQr+eM5IaFY5g6Mtbd5XkUDQ0fdbS6hWc/P8yr247R0mFh0bhErps/hvMmpxASpCdJ9cUYQ15JPS9tOcpbu0pp7bRw9vgkbluSzcLsRN11QUPD5xSWN/LwR4W8v7uMwADhspyRfPesLP20HIL61k5e3HKEZzcWU9XUzrT0GG5bks0l09II8OMDpxoaPqKosolHPyrkzV2lRAQHcv2Zmdy4MJPUWD03wVFtnRbW7DjOU58VcbiqmUmp0fz4womcP3mEX255aGh4ubL6Vn63roDXt5cQGhTIDQvH8L3F2SRE6tmQzmaxGt7JK+UPHxZQXN3CrIw4fnrRRBZmJ7m7NJfS0PBSrR0WnvqsiCc3HMJiDNcvGMNtS7JJjg51d2k+r9Ni5bXcEh5ZX8iJhjaWTkzm3sumkpXku19V96Sh4WWMMbydV8Zv3suntL6NS6ancs+yyYxOiHB3aX6nrdPC85uK+ePHB+nosnLz2Vncfu44IkK892S4gdDQ8CJHqpu55/XdbDpUzZS0GO69bAoLxia6uyy/V9HQxm/e38/rO46TFhvGf31tCsumpfrs8Q4NDS/QZbHy7OeH+f2HBQQFBHDXskmsnJehpz57mNwjNdz75l72ljZw8dRU7rt8mk/uLmpoeLj9Jxr46d/y2H28nvMnp3D/5dP0GxEP1mWx8vRGW8BHhATy88umsnzmSJ/a6jhdaOgZQG5kjOHZjYf5+mOfU1bfymMrZ7H6hjkaGB4uKDCA25Zk894dZzM2KZIfvrKTW17Ipaa5w92luYRuabhJZWM7P/nbLjYUVHL+5BH89pszSIzyvc1cX2exGv7v88M88MEB4iODeeTqWT5xDEq3NDzMhoJKlj3yGZuLqrlv+VRW3zBXA8NLBQYIN589ljWrFhIZEsTK1Zv5w4cF/zLimK/R0HAhYwyPfVzIjf+3lcTIUN7+t7O4/sxMn9oX9ldTR8by9r+dxeWz0nnko0JWrt5MVVO7u8saFhoaLtLU3sVtf8nloXUFXDZjJG+sWsSEFL0C1ZdEhgbx+6tm8rsrc9h5rI7lj33O3tJ6d5fldM5qlnTahkci8lMR2Wm/7RERi4gk2J8rFpHd9ue8+0BFHw5XNXP545+zPr+C/7x0Mo9cPZPwkKGPtak82zfnjOK12xZiNYYrnviCd/PK3F2SUzl8INTe8KgAuADbAMJfAtcYY3rtXSIilwE/Msaca/+5GJhrjKka6Dq96UBo7pEabn7eVuvjK2ezcJx/XcPgzyoa2/j+X7aTe6SWO84bz4/OH+81u6LDfSB0HvaGR8aYDqC74VFfrgFedsJ6Pd4He06wcvUWYsODeWPVIg0MPzMiOoyXbpnPlXNG8ehHhdz19zy6LM7rE+MuzgiNvhohfYWIRGAbXPjvPR42wDoRybU3ROqViNwqIttEZFtlZaUTyh5ez31+mO+/mMuUkTH8/fsLfXpMTtW30KBAHrhiBnecN55Xt5Vw219yae2wuLsshzgjNAbc8Ai4DPjcGFPT47FFxpjZ2Pq5rhKRxb3NaIx5yhgz1xgzNzk52bGKh5Exhj98WMDP397H+ZNTeOnmBfp1qp8TEe68YAL3XT6Nj/ZXcN0zW6hr8d4TwZwRGn01QurN1Zyya2KMKbX/WwGswba745WMMTyw9gCPfFTIlXNG8eR1c/SApzrp+gVj+N+Vs9ldUs81q7dQ66VnkDojNL7E3vBIREKwBcNbp04kIrHAEuDNHo9Fikh0933gQmCPE2pyOWMM97+bzxOfHmLl/Ax++80ZerGZ+opl09N4+ttzOVTZxMqnvTM4HA4NY0wXcDuwFsgHXjXG7BWR20Tkth6TrgDWGWOaezyWAmwUkV3AVuBdY8wHjtbkasYYfvH2Pp7ZeJgbF2byy8un+fX4kur0Fk9I5ukbvDc49NoTJ/jdugP88eODfPesLP7z0sle87Wacq/PCiq5+YVtjEuO4uVbFxAb7jk9WPTak2H0zMbD/PHjg3xr7mgNDDUoiycks/qGuRRWNHLLC9to6/SOb1U0NBzwWm4J972zj2XTUvnVN6ZrYKhBWzIhmd9dNZOth2v44V93esWFbhoaQ/TpgQru+nseZ41L4uGrZ+pBTzVkX88Zyb1fm8IHe09w75t78PRDBr49OuowKSxv5N9e2sGElGj+dP0cQoP0a1XlmJvOyqKisZ0nNxxiVHwE3z8n290l9UlDY5Bqmju46fkvCQ0O5JlvzyUyVN9C5Rx3XTyRktoWHli7n4mpUZw7KcXdJfVKd08Gob3Lwm1/zqW8oZ3VN8xhpB92YVfDR0R48IocpqTF8O8v7+RgRaO7S+qVhsYg3P9OPluLa3joyhxmZcS7uxzlg8JDAnnqhrmEBgdwywu51Ld0urukr9DQGKC3d5Xy581HuOXsLL6eM9Ld5Sgflh4XzpPXzaGktoUf/22nxx0Y1dAYgKLKJu7+ex6zM+L4j4snubsc5QfmZiZwz7LJrM+v4LlNxe4u519oaPSjrdPCD17cTkhQAI+tnE1woL5lyjW+syiT8yeP4Nfv7WfPcc8ZNlD/Avrxm/f3s/9EI7//1kw98KlcqvvAaEJkCLe/tJ2m9i53lwRoaJzWpoNVPLepmBsXZrJ04gh3l6P8UHxkCI9eM4ujNS3c93avI2i6nIZGHxrbOvnpa3mMTYrkLj2OodxoXlYCty7O5pVtx/iswP2j1mlo9OH+d/Ipq2/loatydCAd5XY/PH882cmR3PP6brfvpmho9OIfhZW8su0Yty3JZraej6E8QFhwIA9emUNZfSu/fi/frbVoaJyivcvCvW/uJSspkjvOG+/ucpQ6aXZGPN89K4sXtxxl6+Ga/mcYJq5qlnSOiNT3aJh070DndbWnNhRxuKqZX3x9KmHBuluiPMudF0wkPS6ce9/c47Z2CA6Hhr1Z0uPYRhOfAlwjIlN6mfQfxpiZ9tv/DHJelzha3cJjnxzk0ulpLJ7guSOeK/8VHhLIzy6dzP4Tjby09ahbanBHsyRnzet0//POPoIChP/6mttyS6l+LZuWyqJxiTy09gDVbmgy7cpmSWeKyC4ReV9Epg5y3mFvlrSlqJr1+eWsOnccqbFhTl++Us4iIvz8sqm0dFh4aF2By9fvqmZJ24Exxpgc4I/AG4OY1/bgMDZLMsbw6/f3kxoTxk2Lspy6bKWGw/iUaK5bMIZXtx2jqLLJpet2SbMkY0yDMabJfv89IFhEkgYyryt8sOcEO4/VcecFE/Tgp/Iaq5aOIzQogN9/6NqtDZc0SxKRVLGPuisi8+zrrR7IvMOt02LlwbUHmJASxTfnjHLlqpVySHJ0KDctyuKdvDKXXtDmqmZJVwB77E2RHgWuNja9zutoTYPx1s5Siqqa+cmFE3VwYOV1blk8ltjwYB5ad8Bl63TKAJf2XY73TnnsyR73HwMeG+i8rmK1Gv7304NMSo3mgimeOR6jUqcTGx7M95aM5YEPDpBXUseMUXHDvk6/PiN07d4THKpsZtXScdqzRHmt6xeMITo0iD9tKHLJ+vw2NIwxPPbJQbKSIrlkepq7y1FqyKLDgrl2wRje31NGcVVz/zM4yG9DY+PBKvaWNnDbkrF6LEN5vZsWZRIUEMDqfwz/1obfhsbzm46QGBnC5bN6PZdMKa8yIiaMb8xO52+5JcN+lqhfhkZJbQsf7y/n6nmjtTua8hk3nZVFR5eV13JLhnU9fhkaL26xXeizcv4YN1eilPNMSInmjMx4Xt56FOswNpL2u9Bo67TwypfHOH9yCuk6ULDyMSvnZ1Bc3cIXRdXDtg6/C42P8iuoae7gugW6laF8z7JpacRFBPPiliPDtg6/C401O0pIiQll0bgkd5eilNOFBQfyzdmjWLe3nNrmjmFZh1+FRk1zB58eqGT5zHT9mlX5rBWz0umyGj7Ye2JYlu9XofFuXildVsPlM/VrVuW7po6MYWxSJG/tHJ4Lxv0qNNbsOM6k1GimjIxxdylKDRsR4bKckWw+XE15Q5vTl+83oVHR0Mb2o3VcqqeMKz9wWc5IjIF38sqcvmy/CY2P9lcAcMFUvZpV+b5xI6KYlBrNh/ucf1zDb0Jj/b5yRsWHMzEl2t2lKOUSSyeNYFtxLQ1tnU5drqv6nlwrInn22yYRyenxXLGI7Lb3Q9nmjHpO1dLRxcaDVZw/OUUvgVd+49xJI+iyGjYWVjl1ua7qe3IYWGKMmQHcBzx1yvNL7f1Q5jpaT282HaymvcvK+ZN110T5j1mj44gJC+IT+665s7ik74kxZpMxptb+42ZsAwi7zKZD1YQGBXBGlvZlVf4jKDCAxROS+eRAJcY471oUV/Y96fZd4P0ePxtgnYjkisitfc3kSN+TzUXVzBkTr1e0Kr9z9vgkqpraOVTpvMF5XNX3xDahyFJsoXFXj4cXGWNmY9u9WSUii3ubd6h9T+paOsg/0cCCsYkDnkcpX3FGZgIAXxY7r2G0S/qeAIjIDOBpYLkx5uQleMaYUvu/FcAabLs7TrP1cA3GoKGh/FJWUiRJUSF86cQu867qe5IBvA5cb4wp6PF4pIhEd98HLgT2OKGmfxZXXENIUAA5o2OduVilvIKIMHdMAl8e8aDQGGDfk3uBROB/T/lqNQXYaO+HshV41xjzgaM19ZRXUs/ktBg9nqH81hlZCRyraXXaKeWu6ntyM3BzL/MVATmnPu4sVqthX2kDy2eNHK5VKOXxZoyybWXvLa0nJcbx5uY+fUbokZoWGtu7mDZSd02U/5qUajsLel9pg1OW59Oh0d3fclq6hobyX9FhwYxOCCe/rNEpy/Pp0MgvayAoQJig15soPzc5NYb8Mt3S6NfhqmYyEiIICfLpl6lUvyanxXC4upm2TovDy/Lpv6bDVc1kJkW6uwyl3C4rKRJjbD1/HOWzoWGM4Uh1C5mJGhpKjU6IAOBojYZGn8ob2mnttJCVrKGhVEZ3aFRraPTpmH0zrPvNUsqfJUWFEB4cyNGaVoeX5bOhUdFga4KbEhPq5kqUcj8RYXRC+MkPU0f4bGhUNtpOmU2O0tBQCiApKtQpHeV9NzSa2gkKEOIjQtxdilIeISEyhBondF3z2dCoaGgnKSqUAO2kphQAiZEhVGto9K22pYP4SN3KUKpbQmQojW1ddHRZHVqOz4ZGS4eFyBC9HF6pbgmRwQDUtTq2teGzodHcYSEi1ClX/ivlE8KCbR+i7Z26pdGr1o4uIoJ1S0OpbqHdodHl2PUnrmqWJCLyqP35PBGZPdB5h6qlw0KE7p4odVKo/cLNNndvaQywWdIyYLz9divwxCDmHZKOLqte3apUD92h0e4BB0L7bZZk//kFY7MZiBORtAHOOyQi4MT+MEp5vWO1tlPI80rqHFqOq5ol9TXNgBstDbZZUoAIpvf2K0r5pZRo29nRo+Mdux7LVc2S+ppmwI2WBtssKUAEq2aGUid1Hwh19PwlZ3wnOZBmSX1NEzKAeYfMqvsnSp3UaT+WERzo2FnSLmmWZP/5Bvu3KAuAemNM2QDnHZKQoAA6LRoaSnXrtHSHhmN/9g5vaRhjukSku1lSIPBsd7Mk+/NPYuuJcglwEGgBvnO6eR2tCSAqNIimtk5nLEopn9BhDw1Hv1V0VbMkA6wa6LzOEBUaRFN7l7MXq5TXqm+1fYjGhgc7tByfPZEhOiyIxjYNDaW61bVoaJxWlIaGUv+irqWTyJBAh49p+GxoxEfYBhwx+g2KUoBt9yTOCYNS+WxopMaE0dppoUG3NpQCoKa5nbgIx3ZNwIdDIyXW1h27vKHNzZUo5RlK69pIiw13eDk+Gxpp9tAoq9fQUMoYw/G6VkbFa2j0KTXGFhon6h3v86CUt2to7aKpvYv0OA2NPqXEhBEYIE5pQ6eUtzteZ/vwTNctjb6FBAUwJjGCQxXN7i5FKbc7WmP7O9Ddk35kJ0dxqLLJ3WUo5XYF5U2IwLgRUQ4vy+dDo7i6mS6LYyMVKeXtDpQ3Mjo+gogQx68c8enQGDciik6L4Yge11B+ruBEIxNSop2yLJ8OjUmptjdpX2mDmytRyn3auywcrmpmYqrjuybg46ExMTWa0KAAdh2rc3cpSrlNwYkmuqyGSakxTlmeT4dGcGAAU0fGsMvBgVSV8mbbj9YCMHtMvFOW59OhAZAzOo7dx+v1YKjyW9uP1pISE8pI+1nSjnIoNEQkQUQ+FJFC+79fiTIRGS0in4hIvojsFZF/7/Hcz0XkuIjstN8ucaSe3swcHUdbp5UD5Y3OXrRSXiH3SC2zM+IRcWxs0G6ObmncDXxkjBkPfGT/+VRdwI+NMZOBBcCqUxoi/cEYM9N+c/oIXnMzEwDYXFTj7EUr5fEqGtsoqW1ljpN2TcDx0FgOPG+//zxw+akTGGPKjDHb7fcbgXz66G0yHNLjwslKimTTwSpXrVIpj/HFoWoAzrB/eDqDo6GRYh9VHPu/I043sYhkArOALT0evt3e3/XZ3nZvesw7qGZJPS0al8jmouqTozEr5S82FFQSFxHMtPRYpy2z39AQkfUisqeX26DaJ4pIFPB34IfGmO4TJ54AsoGZQBnwu77mH2yzpJ4WZSfR3GFxuB2dUt7EGMM/Cqs4a1wSgQHOOZ4BAxiN3Bhzfl/PiUi5iKQZY8rsvVkr+pguGFtgvGiMeb3Hsst7TLMaeGcwxQ/UmdmJiMBnBVXMGeO8zTSlPNn+E41UNrazeMLgPmT74+juyVvAt+33vw28eeoEYjtk+wyQb4z5/SnPpfX4cQWwx8F6ehUXEcLcMfGs21fe/8RK+YgNBbbd+MXjPSs0fgNcICKFwAX2nxGRkSLS/U3IIuB64Nxevlp9QER2i0gesBT4kYP19OmiqanklzVwtFqvQ1H+4YM9J5iWHkOqk87P6ObQJW/GmGrgvF4eL8XWUQ1jzEZ6b/SMMeZ6R9Y/GBdNTeX+d/NZu/cEtywe66rVKuUWx+ta2Xmsjv+4eKLTl+3zZ4R2G50QwZS0GNbuPeHuUpQadu/vLgPg0ulp/Uw5eH4TGgAXT0sl92jtyaHPlPJV7+SVMXVkDGMSI52+bL8KjRWz0jEG3thx3N2lKDVsjtW0sPNYHZcMw1YG+FlojE6IYH5WAq/llmjnNeWzXt12jACxfUgOB78KDYAr5ozicFXzycuFlfIlFqvhtdwSFk9IZqQT2hX0xu9C45LpaUSEBPK3bSXuLkUpp/ussJKy+ja+NXf0sK3D70IjMjSIr81I482dpdS3dLq7HKWc6pWtx0iMDOG8ySnDtg6/Cw2AGxdm0dpp4ZVtR91dilJOc7yulQ/zy7lizihCgobvT9svQ2PKyBjmZyXw/KYjOqKX8hnPfX4YgBsWZg7revwyNAC+syiT43WtrM/X61GU92ts6+SvW49x6fQ0p/RrPR2/DY0LpqSSHhfO6n8c1q9fldd75ctjNLZ3cfPZWcO+Lr8NjcAA4XtLxpJ7pPbk6EZKeaOOLivPbjzMvKwEZoyKG/b1+W1oAFw1dzQpMaE88lGhu0tRashe3XaM0vo2fnBOtkvW59ehERYcyG1LstlyuIbNRbq1obxPW6eFxz4+yJwx8Sxx8mA7ffHr0AC4Zl4GydGhPLK+UI9tKK/z161HOdHQxp0XTHBai4L+DHvfE/t0xfbBdnaKyLbBzj+cwoID+cE52XxRVM2nBYMbsFgpd2rtsPD4p4eYn5XAwuxEl63XFX1Pui219zaZO8T5h82188eQmRjBr97N1/M2lNd46rMiKhvb+fGFE122lQEu6HsyzPM7RUhQAHcvm0RhRROv6jUpyguU1bfy5IZDXDI9lXlZrh0s21V9TwywTkRyReTWIczvUN+TgbhoairzMhP4/YcHaGrvcvrylXKm376/H4sx3LNsssvX7aq+J4uMMbOBZdjaMi4ebKGO9D0ZCBHhZ5dOpqqpg4c/LHD68pVyltwjtbyxs5Rbzs5idEKEy9fvkr4n9oGGMcZUiMgaYB7wGTCg+V0lZ3QcK+dn8H+birl8VrpTu1Ip5QxdFiv3vrmHEdGh/OCccW6pwRV9TyJFJLr7PnAh/+xv0u/8rnbXRZOIjwjh/63ZjcWqX8Eqz7L6H4fZW9rAL74+lchQh5oJDJkr+p6kABtFZBewFXjXGPPB6eZ3p9iIYO69bAp5JfX8+Ytid5ej1EmHq5p5eH0BF01NYdkwjf85EK7oe1IE5Axmfne7bEYar+WW8ODaA5w7KYWMRNfvNyrVkzGGe17PIyQogP9ZPs2ttfj9GaG9ERF+tWIaASLc+epO3U1RbveXzUfYXFTDzy6ZTEqMczumDZaGRh9GxUfwi+VT2Xaklic3HHJ3OcqPFZQ3cv+7+ZwzMZlvnTF8Y38OlIbGaayYlc6l09P4w4cF7Dle7+5ylB9q67Rwx8s7iA4L4sErclx65mdfNDROQ0T45YppJEaFcMdfd+hJX8rlfvP+fvafaOTBK3JIjg51dzmAhka/4iJCePhbsyiuauau1/L0SljlMmv3nuC5TcXcuDCTpZP6PFna5TQ0BuDM7ET+4+JJvLu7jGc2HnZ3OcoPHKxo5Mev7mLGqFjuXjbJ3eX8Cw2NAfre4rFcOCWFX7+/n62Ha9xdjvJhDW2d3PrnXEKDAnjyujmEBQe6u6R/oaExQCLCQ1flkJEQwQ9e3E5JbYu7S1I+yGo13PnKLo5Ut/D4tbOHrbWiIzQ0BiEmLJjVN8yhvcvCTc99SUObdmhTzvXgugOszy/nPy+dzIKxrhtYZzA0NAZp3IhonrxuDkWVzax6cTudOmiPcpI/bz7CE58e4pp5Gdw4zA2PHKGhMQSLxiXxqxXT+UdhFf/1xh79RkU57MN95fz3m3s4b9II7ls+1SPOx+iLey6T8wFXnTGaIzXNPP7JIRKjQvjpRZ51hFt5jx1Ha/m3l7czLT2WP66cRVCgZ3+Wa2g44CcXTqSmuYPHPzlEREgQq5a6Z3wD5b32HK/n289uZUR0GM98+wwiQjz/T9LzK/RgIsL9l0+npcPCg2sPEBkSyI2Lhr8tnvIN+080cP0zW4gOC+bFm+d7zBmf/dHQcFBggPDQlTm0dFj4+dv7CA0O5Jp5Ge4uS3m4gxWNXLt6CyFBAbx0y3y3DNs3VJ698+QlggMDeGzlLM6ZmMw9r+/Ws0bVaR040cg1q7cgIrx0ywLGJEa6u6RBGfZmSSIy0d4kqfvWICI/tD/3cxE53uO5Sxypx51CgwL50/VzuHhqKve9s4/HPtaObeqrth+t5ao/fYEAL98yn+zkKHeXNGjD3izJGHPA3iRpJjAHaAHW9JjkD93PG2PeO3V+bxIaFMhjK2fxjVnpPLSugN9+cECDQ530j8JKrnt6C3ERwfz9+wsZnxLt7pKGxNFjGsuBc+z3nwc+Be46zfTnAYeMMUccXK/HCgoM4KErcwgPCeTJDYcob2jjN9+cTmiQZ10/oFzr3bwyfvjKDrKTo3jhu/MYEe3e0bcc4apmSd2uBl4+5bHbRSRPRJ49XS/X4W6W5EwBAcL9l0/jJxdOYM2O41z/zFbqWjrcXZZyA2MMj31cyKqXtpMzKo5XvnemVwcGgPS3+Swi64HUXp76GfC8MSaux7S1xpi+mkCHAKXAVGNMuf2xFKAKWwe2+4A0Y8xN/RU9d+5cs23btv4m8whv7jzOT/+WR3p8OP934xlkJnnXQS81dG2dFu7+ex5v7Czl8pkj+c03Z3jcFat9EZHcU/oun+SSZkl2y4Dt3YFhX/bJ+yKyGninv3q8zfKZ6YyMC+fWF7ax/PHPefjqmSyd6DkDqqjhUdHYxvf+nMuOo3X89KKJ/OCcbI8+NXwwhr1ZUg/XcMquiT1ouq3gn02UfMoZmQm8sWoRabFh3PTcl/zhwwKsOsK5z9pcVM2lj24kv6yBJ66dzaql43wmMMA1zZIQkQj786+fMv8DIrJbRPKApcCPHKzHY41JjGTNDxaxYlY6j3xUyHee+5LaZj3O4UusVsPjnxxk5erNRIcGseYHi9za1Gi49HtMwxN50zGNUxljeGnrUX7x1j4So0L43ZU5LByX5O6ylIOqm9q589VdbCio5Os5I/nVN6YT5aa2ic5wumMaekaoi4kI184fw2vfP5Pw4EBWPr2F+9/ZR1unxd2lqSFau/cEFz38GV8UVfPLFdN45OqZXh0Y/dHQcJMZo+J4546zuG5BBk9vPMzyxz5nX2mDu8tSg1Df2smdr+7ke3/OZUR0GG/dvohr54/xqeMXvdHdEw/wyYEK/uO1PGqbO7h18VjuOG+813w1568+2V/BPa/vprKpnVVLx3H70nGEBPnOZ/Dpdk80NDxETXMHv3ovn9dyS8hIiOCXK6Zx9vhkd5elTnG8rpX/eXsva/eWM25EFL+/KocZo+LcXZbTaWh4kU2HqvjPNXsoqmrm6zkjueeSSaTFet6I1P6mo8vKs58f5pH1hRgMd5w3npvPGutTWxc9aWh4mbZOC098eognNhwiQOCWs8fyvSXZPn1wzVMZY1i3r5zffrCfospmLpySwr2XTWFUvPeMfzEUGhpe6lhNCw+uPcBbu0pJigrlzgsmcNXcUR4/hqSv2FZcw6/f30/ukVqykyP52aWTOXdSirvLcgkNDS+342gtv3w3n21HaslMjGDV0nFcPiudYA2PYbHneD0Pry9kfX45KTGh/Oj8CVwxx7/CWkPDB3RvJj/6USF7SxvISIhg1dJsVswa5bP71a629XANj39ykA0FlUSHBXHbkmxuWpRFeIj/fZOloeFDjDF8lF/Box8XkldSz8jYMK4/M5Nr5o0mLiLE3eV5HYvV8OmBCv60oYitxTUkRoZw01lZXH/mGGLCgt1dnttoaPggYwyfHqjkqc+K+KKomrDgAFbMGsV3FmUywUtHhHKl2uYOXtl2jL9sPkJJbStpsWHcungsV5+R4ZdbFqfS0PBx+WUNPL+pmDU7jtPeZeWMzHiumDOKS6anEe3Hn5ansloNW4tr+Nu2Et7OK6Wjy8r8rASuP3MMF01N1WNEPWho+Ima5g5e+fIYr+Ue41BlM2HBASyblsY3Zqdz5thEvzqQ11NBeSNrdhznzR3HKa1vIzIkkBWz07l+QSYTU3WrrDcaGn7GGMPOY3W8llvCW7tKaWzrIi4imPMnp3DR1FTOHp/k06epG2PYc7yBD/PLWbf3BPtPNBIYICwen8Tls9K5YEqKV3QycycNDT/W1mnh0wOVrNt7gvX55TS0dREREsjZ45M4a1wSC8clMTYp0usvsqpv6WTL4Wo2FFSyPr+c8oZ2AgTmjIln2bQ0LssZ6TUdzDyBQ8P9Ke8WFhzIxdNSuXhaKp0WK5uLqvlgzwk+PVDJ2r220RbTYsNYmJ3E/KwEckbHMW5EFIEBnh0iVU3t7Dhax+aiajYXVbOvrAFjICIkkCUTkjl/cgpLJ40gIVK/UXI2h7Y0RORK4OfAZGCeMabXj38RuRh4BAgEnjbGdI/wlQC8AmQCxcBVxpja/tarWxqOM8ZwtKaFzw9W8/nBKjYdqqK2pROw/eFNS49l5ug4JqVGk50cRfaIKLecxt5lsVJS28rBiib2lNaz53gDe47Xc6KhDYDQoABmZ8SzYGwiC8YmMDMjTttFOMGw7Z6IyGTACvwJ+ElvoSEigUABtuH+SoAvgWuMMftE5AGgxhjzGxG5G4g3xpyubwqgoTEcrFZDcXUzu0rq2HWsnp3H6thX2kCHxXpymtSYMLJHRJIeF05qbDhpsWGkxoSREhNGbEQw0WFBRIYEDWgrxRhDh8VKc7uFqqZ2KhvbqWhso7KxnRP17RRXN1Nc1czRmha67OOpisDYpEimp8cyLT2WGaPimDEq1qePz7jLsO2eGGPy7Ss43WTzgIPGmCL7tH/F1mRpH4NvtqSGSUCAMDY5irHJUayYNQqATouVI9UtHKxo4lBlE4cqmjhU1cynByqpbGqnr8+bqNAgIkMDCQoIsC8bAuz/R1o7LLR2WGjptGDpY3Dl8OBAxiRGMDE1moumpZKVFMnYpEgmp8UQqRftuZ0rfgPpwLEeP5cA8+33/6XZkoj0Oba/iNwK3AqQkaFd2V0hODCAcSOiGDfiq/1GOy1WKhrbOVHfRnlDG41tnTS2dZ28Nbd3YTEGqzFgwGoMBlsghIcEEhESSERIEBEhgSRFhTIiOpRk+y0qNMjrD8z6sn5D43TNkowxp2tZcHIRvTw26H0iY8xTwFNg2z0Z7PzKuYIDA0iPCyc9Tsf68DcONUsaoBJgdI+fR2HrtAYwmGZLSikP4IpTBL8ExotIlr0149XYmizB4JotKaU8gEOhISIrRKQEOBN4V0TW2h8/2SzJGNMF3A6sBfKBV40xe+2L6LXZklLKc+kZoUqpr9BmSUopp9HQUEoNioaGUmpQNDSUUoPilQdCRaQSODKASZOAqmEuZ7j5wmsA33gdvvAaYGCvY4wxptcWf14ZGgMlItv6OgLsLXzhNYBvvA5feA3g+OvQ3ROl1KBoaCilBsXXQ+MpdxfgBL7wGsA3XocvvAZw8HX49DENpZTz+fqWhlLKyTQ0lFKD4lOhISJXisheEbGKSJ9fKYnIxSJyQEQO2scm9RgikiAiH4pIof3f+D6mKxaR3SKyU0Q84uq9/t5XsXnU/nyeiMx2R539GcDrOEdE6u3v/U4RudcddZ6OiDwrIhUisqeP54f+uzDG+MwN26joE7GNNTq3j2kCgUPAWCAE2AVMcXftPep7ALjbfv9u4Ld9TFcMJLm73sG8r8AlwPvYRnNbAGxxd91DfB3nAO+4u9Z+XsdiYDawp4/nh/y78KktDWNMvjHmQD+TnRzo2BjTAXQPdOwplmMbZBn7v5e7r5RBGcj7uhx4wdhsBuLsI7Z5Ek///zEgxpjPgJrTTDLk34VPhcYA9TbQcbqbaunNvwy2DPQ12LIB1olIrn3QZXcbyPvq6e89DLzGM0Vkl4i8LyJTXVOaUw35d+F148F7ykDHjjjdaxjEYhYZY0rtI7h/KCL77Z8u7jKQ99Xt7/0ADKTG7diuzWgSkUuAN4Dxw12Ykw35d+F1oWGGd6BjlzjdaxCRAQ22bIwptf9bISJrsG1WuzM0BvK+uv29H4B+azTGNPS4/56I/K+IJBljvOlitiH/Lvxx9+R0Ax17gn4HWxaRSBGJ7r4PXAj0epTchQbyvr4F3GA/cr8AqO/eFfMg/b4OEUkVe2MWEZmH7e+o2uWVOmbovwt3H+V18hHjFdgStB0oB9baHx8JvHfKkeMCbEfJf+buuk95DYnAR0Ch/d+EU18DtiP7u+y3vZ7yGnp7X4HbgNvs9wV43P78bvr4hsvdtwG8jtvt7/suYDOw0N019/IaXgbKgE7738R3nfW70NPIlVKD4o+7J0opB2hoKKUGRUNDKTUoGhpKqUHR0FBKDYqGhlJqUDQ0lFKD8v8BMRrY65oV+LAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# RK4 method\n", "n=1000\n", "t=np.linspace(0,2*np.pi, n, endpoint=False)\n", "dt= t[1]-t[0]\n", "\n", "xcord = np.zeros(n)\n", "ycord = np.zeros(n)\n", "tcord = np.zeros(n)\n", "#print(dt)\n", "\n", "x0=1\n", "y0=0\n", "vx0=0\n", "vy0=1 # v = dx/dt = r0 /(r/vI) = vI\n", "\n", "def pag(x,y,t):\n", " r=np.sqrt(x**2+y**2)\n", " ax = -x/r**3\n", " ay = -y/r**3\n", " return ax,ay\n", "\n", "i=0\n", "for t0 in t:\n", " ax1,ay1 = pag(x0,y0,t0)\n", " vx1,vy1 = vx0,vy0\n", " \n", " ax2,ay2 = pag(x0+vx1*dt/2,y0+vy1*dt/2,t0+dt/2)\n", " vx2,vy2 = vx0+ax1*dt/2,vy0+ay1*dt/2\n", " \n", " ax3,ay3 = pag(x0+vx2*dt/2,y0+vy2*dt/2,t0+dt/2)\n", " vx3,vy3 = vx0+ax2*dt/2,vy0+ay2*dt/2\n", "\n", " ax4,ay4 = pag(x0+vx3*dt,y0+vy3*dt,t0+dt)\n", " vx4,vy4 = vx0+ax3*dt,vy0+ay3*dt\n", " \n", " vx0 = vx0+(ax1+2*ax2+2*ax3+ax4)*dt/6\n", " vy0 = vy0+(ay1+2*ay2+2*ay3+ay4)*dt/6\n", " x0 = x0+(vx1+2*vx2+2*vx3+vx4)*dt/6\n", " y0 = y0+(vy1+2*vy2+2*vy3+vy4)*dt/6\n", " \n", " xcord[i]=x0\n", " ycord[i]=y0\n", " tcord[i]=t0\n", " i=i+1\n", "\n", " \n", "plt.plot(xcord,ycord)\n", "ax=plt.gca()\n", "ax.set_aspect(\"equal\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Trijų kūnų problema (sunkesnė)\n", "Išspręskite Žemės judėjimo aplink Saulę judėjimo lygtį, su aplink Saulę besisukančiu Jupiteriu. Laikykite, kad Jupiteris juda apskritimu aplink Saulę ir nekreipkite dėmesio į Žemės įtaką jo orbitai.\n", "\n", "Raskite chaotinio judėjimo pavyzdžių. Nubrėžkite orbitų pavyzdžių. Patikrinkite, kad orbitos yra skaitmeniškai stabilios (žingsniai pakankamai maži)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. n-matės sferos tūris\n", "\n", "Raskitė n–matės vienetinės sferos tūrį ($n=2$ -- skritulio plots, $n=3$ sferos tūris, $ V_n = \\frac{\\pi^{n/2}}{(n/2)!} $)\n", "\n", "Pasinaudokite $r<1$ formule ir sumuokite tūrio elementus $[-1,1]^n$ kube.\n", "\n", "- Kiek taškų reikia padalinti kubą vienodais tūriais, kad pasiektumėte vieno procento tikslumą? \n", "- Kiek taškų reikia, jei taškus pasirenkat atsitiktinai? Pvz, su np.random.rand\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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.9.1" } }, "nbformat": 4, "nbformat_minor": 4 }