{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "数値計算試験問題\n", "
\n", "
\n", "
\n", "2021/12/17 実施\n", "
\n", "2022/12/16 予行演習実施\n", "
\n", "cc by Shigeto R. Nishitani 2021-2 \n", "
\n", "\n", "2022/12/16 予行演習:以下の問いに答えよ.答案はpdfとipynb形式でLUNAのd12へ全員が個別に提出せよ.pdfは2pageを一枚に集約して作成すること.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1 簡単な行列計算(25点)\n", "次のデータにフィットした二次関数を求める.\n", "``` python\n", "import numpy as np\n", "\n", "xdata = np.array([1,2,3,4])\n", "ydata = np.array([1,2,5,11])\n", "```\n", "\n", "最小二乗法の正規方程式(normal equations)から求められるデザイン行列$A$は,\n", "$$\n", "A=\\left( \\begin {array}{ccc} 1. & 1. & 1.\\\\\n", " 1. & 2. & 4. \\\\\n", " 1. & 3. & 9. \\\\\n", " 1. & 4. & 16.\n", "\\end {array} \\right)\n", "$$\n", "となる.$A^TA$の逆行列から\n", "\n", "$$\n", "a = (A^TA)^{-1} A^T y\n", "$$\n", "により最適パラメータ$a$を求め,データと同時に plot せよ." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/bob/opt/anaconda3/lib/python3.8/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.4\n", " warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 1. 1.]\n", " [ 1. 2. 4.]\n", " [ 1. 3. 9.]\n", " [ 1. 4. 16.]]\n", "[ 2.75 -2.95 1.25]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfdElEQVR4nO3de5zWY/7H8dfVWYZImaITOlha1FRCS1Pr7Cdsu6Io2R02h6iQoiyl8qt2rYTIod1qtPppMzpIjTOlSVQSRSWhomhMOs31++O6o6apmft43d/7fj8fj/vhPs73/fjq/sw11/c6GGstIiISPBV8BxARkciogIuIBJQKuIhIQKmAi4gElAq4iEhAVUrkwWrVqmUbNWoU0Wd/+uknDj300NgGigHlCo9yhUe5wpOquQoKCjZZa2vv94K1NmG3rKwsG6n8/PyIPxtPyhUe5QqPcoUnVXMBC20pNVVdKCIiAaUCLiISUCrgIiIBpQIuIhJQKuAiIgGlAi4iEi8TJ0KjRpzToQM0auQex1BCx4GLiKSNiRMhJweKijAAa9a4xwBdu8bkEGqBi4jEw8CBUFS073NFRe75GFEBFxGJh7VrAfiGTG7j73xHzX2ejwUVcBGReGjQAIB7GMJYerGZI/d5PhZUwEVE4mHoUBZXa8vT9OQWHqExq6B6dRg6NGaHUAEXEYkDe3VXbj/uRWpW2MI9DIGGDWHcuJhdwASNQhERiYv//hdeW16HRx+FD096kfbt28f8GGqBi4jE2I4dcMcdcNJJv44cjAe1wEVEYmzMGFi5EmbOhEpxrLJqgYuIxNCmTXD//XDBBe4WTyrgIiIxNHgwFBbCqFHxP5YKuIhIjHz8MTzxBNx4o+v/jjcVcBGRGOnbFzIy4L77EnM8XcQUEYmBWbPcbdQoqFUrMcdUC1xEJEq7dkGfPtC4Mdx8c+KOqxa4iEiUnngCli+HF1+EKlUSd1y1wEVEorB5sxt5kp0NnTol9tgq4CIiURgyBL7/HkaPBmMSe2wVcBGRCH32GTzyCPTsCaedlvjjl1nAjTFPG2M2GGOW7vVcTWPMHGPMZ6H/HhnfmCIiyefOO6FqVdcK96E8LfBngZITQvsDc621TYC5occiImkjPx+mTYO774Y6dfxkKLOAW2vfAL4v8XQn4LnQ/eeAy2IbS0Qkee3eDbff7pb4vv12fzmMtbbsNxnTCMiz1jYPPd5irT0idN8Am/c8LuWzOUAOQGZmZlZubm5EQQsLC8nIyIjos/GkXOFRrvAoV3gSlevll+swcuSJ3HvvMjp02Bj3XNnZ2QXW2lb7vWCtLfMGNAKW7vV4S4nXN5fn52RlZdlI5efnR/zZeFKu8ChXeJQrPInI9eOP1mZmWnvmmdYWF5fvM9HmAhbaUmpqpBN5vjXG1LXWfm2MqQtsiPDniIgEyvDh8O23MH164ocNlhTpMMLpQPfQ/e7Af2MTR0Qkea1e7dY66doV2rTxnaZ8wwgnA+8CzYwx64wx1wPDgXONMZ8Bvw89FhFJaf37Q4UKMGyY7yROmV0o1tqrDvBSxxhnERFJWu+8A88/D4MGQf36vtM4mokpIlKG4mI3XPCYY9zknWSh1QhFRMoweTIsWADPPguHHuo7za/UAhcROYiiItf3nZUF11zjO82+1AIXETmIYcNg3TqYNMldwEwmSRZHRCR5fPIJjBgB3brB737nO83+VMBFREphLfTq5fq8R470naZ06kIRESnFpEluxcGxYyEz03ea0qkFLiJSwubNbpPiNm0gJ8d3mgNTC1xEpISBA2HTJpg1CypW9J3mwNQCFxHZy4IF8PjjcMst0KKF7zQHpwIuIhKyaxfceCPUrQv33+87TdnUhSIiEjJ2LHzwAUyZAocf7jtN2dQCFxEB1q+He+6B88+Hzp19pykfFXAREdxiVTt2wJgx/jdqKC8VcBFJe6+84rpNBg6Exo19pyk/FXARSWvbtrkZl02bJtdSseWhi5giktaGD4dVq+DVV6FqVd9pwqMWuIikrU8/dQX86quhYwD3GFMBF5G0tGexqkMOcRsVB5G6UEQkLeXmwty58OijUKeO7zSRUQtcRNLOli1usapWreCGG3yniZxa4CKSdu65BzZsgLy85F6sqixqgYtIWlm40E2Zv+kmt89lkKmAi0ja2L3bLVaVmQkPPOA7TfTUhSIiaeOxx6CgwF3ArFHDd5roqQUuImnh66/dVPlzz4U//cl3mthQAReRtNCnD2zf7oYNBmWxqrKogItIypszx3Wb3H03NGniO03sqICLSEr7+Wc34qRJE7jrLt9pYksXMUUkpQ0ZAp995paMrVbNd5rYiqoFboy53RizzBiz1Bgz2RiTYqdHRILs/ffdYlXXXecuXqaaiAu4MeZY4FaglbW2OVAR6BKrYCIi0fj5Z+jRw21QPHq07zTxEW0XSiXgEGPMTqA6sD76SCIi0bvvPvj4Y5g1C444wnea+Ii4BW6t/QoYCawFvgZ+sNa+EqtgIiKReu89+N//hT//2W1SnKqMtTayDxpzJDAVuBLYAvwHeMFa++8S78sBcgAyMzOzcnNzIzpeYWEhGRkZEX02npQrPMoVHuUKT2FhIZUrH05OThY//1yRp59+n0MP3e07VtTnKzs7u8Ba22q/F6y1Ed2APwLj93p8LTD2YJ/JysqykcrPz4/4s/GkXOFRrvAoV3jy8/Ntv37WgrVz5vhO86tozxew0JZSU6MZhbIWaGuMqW6MMUBHYHkUP09EJCpLlx7OqFFuwarf/953mviLpg98PvACsAhYEvpZ42KUS0QkLEVFMGLEiTRsCA895DtNYkQ1CsVaOxgYHKMsIiIRGzgQ1q2rzrx5cNhhvtMkhqbSi0jgvfEGPPwwXH75OrKzfadJHBVwEQm0n35yMy2POw7+8pfPfcdJKK2FIiKB1r8/fP45vP46FBcX+46TUGqBi0hgvfYajBkDvXvD2Wf7TpN4KuAiEkiFha7rpHFjePBB32n8UBeKiATSnXfCmjXw5ptQvbrvNH6oBS4igfPqq26D4j594KyzfKfxRwVcRALlxx/h+uuhWTN44AHfafxSF4qIBEq/frBuHbz9NhxyiO80fqkFLiKBMXs2PPmkK+Jt2/pO458KuIgEwpYtruvkN7+Bv/3Nd5rkoC4UEQmEPn3gm2/gxRdTb3PiSKkFLiJJb8YMeOYZuOsuaN3ad5rkoQIuIklt82b4y1+geXMYNMh3muSiLhQRSVrWQk4ObNgA06dD1aq+EyUXFXARSVpPPAEvvAAjRkBWlu80yUddKCKSlD76CG67ze0q36+f7zTJSQVcRJJOYSFceSXUrAkTJkAFVapSqQtFRJLOzTfDihUwdy4cfbTvNMlLv9dEJKn861/w3HNw772k1fZokVABF5GksWIF/PWvbnOGe+/1nSb5qYCLSFL4+WfX712tGkyaBJXUwVsmnSIRSQr9+sGHH0JeHhx7rO80waAWuIh4N3UqPPqoW+/k4ot9pwkOFXAR8Wr1arfKYOvWMGyY7zTBogIuIt7s3AlXXeWmzOfmQpUqvhMFi/rARcSbe+6B996DKVPg+ON9pwketcBFxItZs+Chh+CGG+CPf/SdJphUwEUk4davh2uugd/+Fv7+d99pgksFXEQSavdu6NYNiorg+ee1MXE01AcuIgk1dCjk57sddn7zG99pgi2qFrgx5ghjzAvGmE+MMcuNMWfEKpiIpJ7XX3cbEnfrBt27+04TfNG2wB8GZllrOxtjqgDVY5BJRFLQxo1w9dVwwgkwdiwY4ztR8EVcwI0xNYCzgR4A1todwI7YxBKRVFJcDD16wKZN8PLLcNhhvhOlBmOtjeyDxpwGjAM+Bk4FCoDe1tqfSrwvB8gByMzMzMrNzY3oeIWFhWRkZET02XhSrvAoV3hSJdeUKfV47LHG3Hrrp1x++fqkyZUo0ebKzs4usNa22u8Fa21EN6AVsAs4PfT4YeCBg30mKyvLRio/Pz/iz8aTcoVHucKTCrnefdfaSpWsvfxya4uL45fJ2tQ4X6UBFtpSamo0FzHXAeustfNDj18AWkbx80QkxXz1FVxxBdSvD+PHq9871iIu4Nbab4AvjTHNQk91xHWniIiwbRtcdhls3QrTp8ORR/pOlHqiHYVyCzAxNALlc+C66COJSNBZC3/+MxQUwLRp0Ly570SpKaoCbq1djOsLFxH5xYgRbledoUPh0kt9p0ldmkovIjH10kswYAB06QJ33+07TWpTAReRmFm2zE3WadlSFy0TQQVcRGLiu+9cd0lGhuv3rq552XGnxaxEJGo7d7o1vb/6Cl57DerV850oPaiAi0jUbr/drTD43HPQtq3vNOlDXSgiEpUnnnA7yvfrB9de6ztNelEBF5GIvf463HwzXHghDB/uO036UQEXkYh88QX84Q/QuDFMngwVK/pOlH5UwEUkbEVFFbn0Urc92vTpUKOG70TpSRcxRSQsxcXw4IO/YflymDkTmjTxnSh9qYCLSFgGD4a3367Fww/Duef6TpPe1IUiIuX2/C1vMWQIXM9T3DKqEUyc6DtSWlMLXETKZdGQGVw3pj3teJOx9MKs3Qk5Oe7Frl39hktTaoGLSJm++QY63XcatdnIVP5AFXa6F4qKYOBAv+HSmFrgInJQ27e7XXW+312DtzmLo9m47xvWrvUTTNQCF5ED273bza58912YUKsvp/Hh/m9q0CDxwQRQAReRA7AWevWCKVNg5Ej4wz9+t/8Sg9Wru10bxAt1oYhIqQYOhHHj3KYMffsCdP3lBbt2LaZBA1e8dQHTG7XARWQ/o0bBsGFwww0lGthdu8Lq1bw+bx6sXq3i7ZkKuIjs45ln3MqCV17pVhnUrjrJSwVcRH7x4otuN/nzz4cJE7RAVbJTARcRAObOdRsRn346TJ0KVar4TiRlUQEXEd5/Hy67DJo2hbw8OPRQ34mkPFTARdLc8uVuQ4bateGVV6BmTd+JpLxUwEXS2Jo1bkXBypVhzhyoW9d3IgmHxoGLpKkNG1zx/ukntzXaCSf4TiThUgEXSUM//AAXXADr1sGrr8Ipp/hOJJFQARdJM9u2wf/8Dyxd6rZDO/NM34kkUirgImlk507405/grbfcRsQXXOA7kURDBVwkTRQXQ8+ebpjgY4+5mZYSbFGPQjHGVDTGfGCMyYtFIBGJPWvh9tvh3/+GIUPgxht9J5JYiMUwwt7A8hj8HBGJkwcegH/+0xXxAQN8p5FYiaqAG2PqARcDT8UmjojE2siRbif57t3dfS1OlTqibYH/A7gTKI4+iojEkrUwaBDccYe7cPnUU1BBU/dSirHWRvZBYy4BLrLW9jLGtAf6WWsvKeV9OUAOQGZmZlZubm5ExyssLCQjIyOiz8aTcoVHucITaa7iYhg7tjFTp9bjoou+pk+fFTFdWTDVzle8RZsrOzu7wFrbar8XrLUR3YBhwDpgNfANUAT8+2CfycrKspHKz8+P+LPxpFzhUa7wRJJr1y5rr7vOWrD2ttusLS5OjlyJkKq5gIW2lJoa8R9U1tq7rbX1rLWNgC7APGttt0h/nohEb8cOuOoqtynD4MEwerT6vFOZxoGLpIiiIujcGWbOdFui9enjO5HEW0wKuLX2NeC1WPwsEQnfjz/CJZe4GZZPPul21ZHUpxa4SMBt2uTW81682E2P1wzL9KECLhJg69e7JWE//xymTYOLL/adSBJJBVwkoL74An7/e7eu98yZ0L6970SSaCrgIgH08ceu5b1tm9uMuE0b34nEB83LEgmYRYvgnHPcZJ033lDxTmcq4CIB8tZbkJ3tdo1/801o3tx3IvFJBVwkIGbPhvPOcxsPv/kmNG7sO5H4pgIuEgBTp7pt0Jo1c90m9ev7TiTJQAVcJMk9/rhbTbB1a8jPh6OP9p1IkoVGoYgkqR07YNSopuTluYk6//mP6/sW2UMtcJEk9M037mJlXt4x9O8PL72k4i37UwEXSTILFkCrVm5q/KBByxg2jJiu5S2pQwVcJIk89xycfTZUrgzvvAPZ2Rt9R5IkpgIukgR27oTevaFHDzjrLFi4EE491XcqSXYq4CKebdwI55//667xs2fDUUf5TiVBoFEoIh4tXgyXXeYuWk6YANdc4zuRBIla4CKe5ObCmWfC7t1uiryKt4RLBVwkwXbvhrvucntXZmW5/u5W++83LlKmQBTw2bNhwYIjfccQidrmzW7ThYcegl693FKwmZm+U0lQBaIP/KGHYN68U/nySxgxAqpX951IJHzLlkGnTrB2rfatlNgIRAv85Zehc+cvGTPG/clZUOA7kUh4XnwRTj8dfvoJXntNxVtiIxAFvFo1uOmmVbz6KmzdCm3bwpAhsGuX72QiB7dtG/TtC1dc4dbuLihwFy5FYiEQBXyPjh1hyRLo3BnuvdfNWFu1yncqkdK9+y60aAGjR8Nf/+pa3scc4zuVpJJAFXCAI4+EyZNh0iS3L+Bpp8H48WCt72QizrZtcMcd0K6duz9nDowd6/6SFImlwBXwPa66yrXGW7d2/YmXX+525xbxaf58aNkSRo50/y6XLHE7x4vEQ2ALOLhdSV591f2JOmsW/Pa3kJfnO5Wko59/dmO7zzzTXaicPRueeAIOP9x3MkllgS7gABUquPUjFi6EOnXctlM33ACFhb6TSbp4/33X6n7oIejZ07W6zzvPdypJB4Ev4Hs0b+7WUb7jDjfGtkULeO8936kklW3fDnff7UZFbd3q/gp88kmoUcN3MkkXKVPAAapWda2g/Hy3HVW7djB4sFuqUySW9rS6hw93S8AuXepWFBRJpJQq4Huccw589BFcfTXcf79bX3nFCt+pJBVs3w4DB8IZZ8CWLW6S2fjxanWLHylZwMF9oSZMgClTYOVKd4Gzb1+3FoVIJAoK3KJTDz7oVg5ctgwuush3KklnERdwY0x9Y0y+MeZjY8wyY0zvWAaLlT/+0X3Rrr0W/v53aNzYLZyvbhUBYOJEaNSIczp0gEaN3OMSduxwE8dOPx2++86NdHrmGTjiiISnFdlHNC3wXUBfa+1JQFvgJmPMSbGJFVt168JTT8EHH7iLm717u4ue//2vJgCltYkTIScH1qzBWAtr1rjHoSJeXOzunniiW7qha1fXGLj4Ys+5RUIiLuDW2q+ttYtC97cCy4FjYxUsHk491c2Ky8tzww8vuww6dIBFi3wnEy8GDoSion2fKyrCDhjIjBnul323bq47btYst+HwkVrVWJKIsTFoghpjGgFvAM2ttT+WeC0HyAHIzMzMys3NjegYhYWFZGRkRJn0V7t2GV56qS7PPnscW7dW4rzzvuX66z+ndu0dXnPFinKV7ZwOHVzLey/vcAb9Gc6bnM0xx2yjZ88vyM7eQAVPV4uS6XztTbnCE22u7OzsAmvt/tt+WGujugEZQAFwRVnvzcrKspHKz8+P+LMHs3mztXfcYW2VKtZWr27t4MHWFhb6zxUt5SqHhg2tdb1odgkn20uZZsHazAob7KOPWrt9u++ASXa+9qJc4Yk2F7DQllJTo2pXGGMqA1OBidba/4vmZ/lyxBFu7Pgnn8All8Df/gZNmriLVLt3+04ncTV0KGuqNaM7z3IKH/Ea7RlaeTCrxs2lVy+oUsV3QJGDi2YUigHGA8uttaNjF8mP446D55+Ht9+GBg3clOhWrWDePN/JJB42boTb3u9K013LeJ4u9GE0n9c7hwHPNOXQ67v4jidSLtG0wM8CrgE6GGMWh26BHxV75pluHefJk92Y8Y4d4dJL3dK1Enxbt7q/so4/Hh55BK7pXpHP1lblkvxWHPXlYjfURCQgohmF8pa11lhrT7HWnha6zYhlOF+MgS5dXLfK8OFuIf6TT3ZTpfPy3PAyCZbt2934/xNOgPvuc/8vly1zw0vr1/edTiQyKTsTMxaqVXNLhK5a5abkL1niVjts0sQtYatZnclvyxYYM8aN5e7d283InT8fXnjBPScSZCrg5VC7tpuJt2YN5Oa6iUF9+0K9ejBqVFOWLPGdUPZmLbz1FnTv7rYwu+UWOPpoeOUVt358mza+E4rEhgp4GCpXhiuvdMVh0SLXzfLKK5mccgpkZ8PUqdpo2afvvnPLJZx8Mvzud24n+O7d3Rom8+fDuee67jGRVKECHqEWLdwqdFOmvMuIEfDFF26z5eOPh2HDYNMm3wnTg7Vu+eCrrnKt7T593MzJ8eNh/Xp47DG37KtIKlIBj1KNGru4807XTz5tGjRtCgMGuO6VHj1c609i79tvYcQId747dHBT3W+4wS0j/O67bhhoEk7IE4kpFfAYqVgROnVyfazLlsH117sLZa1auaGJzzzjio5ErrjY7TXZubP7Bdm/v7seMWGCa23/85/uIqVIulABj4OTToJHH4WvvoJ//MN1p/Ts6fbsbN0aBg1y271ppmfZrHWbcQwZ4oYAXnCBG9Z5661ubP4bb7i1uQ85xHdSkcSr5DtAKqtRww1du+UW+PBDmDHD3YYOhQcegKOOcuORL7rI/bdWLd+Jk8OXX8LcuW4W7Lx57hchuK6S4cPdKpJVq3qNKJIUVMAToEIFd9GzRQu3gun337shbTNmuL7bSZPc6Ig2bVwxv/BCyMrC2wp4ibZhg2tV7ynaK1e652vVckW7Y0e3y3ujRj5TiiQfFXAPatZ0QxC7dHH9ugUFrpjPnOlmCQ4e7MaeX3ihu513nvtMqvjhB9f1MW+eK9p7xtEffrjbz/Smm1zRPvnk9PklJhIJFXDPKlRw/eKtW7vCvXGju1A3c6abtj9hgnvPqae6mYPNmrlb06buluwjLax147MXL3YFe9q0lqxY4X5xVasG7dq5IYAdOri/OirpX6RIuenrkmRq13a7wHTr5i5yLljgWucLFsA777iZoHvvQXDssb8W9GbNoNn6fJpOHEy79e9Aw3quwz2OCzQVF7vRNatXu5mqa9bsf3/PpjcVK8KJJ1oGDHAt7LZtXREXkciogCexihXhjDPcbY9t21wf8YoV8Omn7r8rVrjCvmULQDaQTRW203jNSpp2X0Wz3GXUO/9kKlXil1vFiqXfP9BrRUWlF+e1a92mv3urWRMaNnS/UM4779f77dpBQcEHtG/fPkFnUCS1qYAHzCGHuLHOJcc7WwubGrRkxbrqrKAZn9KUFTTjk93NeDnvBHbmxS5DnTquKLdsCVdc4e43bOguMjZoAIcdFrtjiciBqYCnCGOg9leLqY2lHW/v89ouKvH9tzvZvdt1y+za5W4Hul/aa1WruiLdoIG6PUSShQp4KmnQwPVtlFCp4bEcfbSHPCISVxqklUqGDoXq1fd9rnp197yIpBwV8FTStSuMGwcNG2KNcX0e48ZpmzCRFKUCnmq6doXVq3l93jw3TETFWyRlqYCLiASUCriISECpgIuIBJQKuIhIQKmAi4gElLF7r4wU74MZsxHYf6ZJ+dQCknGrYOUKj3KFR7nCk6q5Glpra5d8MqEFPBrGmIXW2la+c5SkXOFRrvAoV3jSLZe6UEREAkoFXEQkoIJUwMf5DnAAyhUe5QqPcoUnrXIFpg9cRET2FaQWuIiI7EUFXEQkoJKugBtjLjDGrDDGrDTG9C/l9arGmOdDr883xjRKklw9jDEbjTGLQ7c/JyDT08aYDcaYpQd43Rhj/hnK/JExpmW8M5UzV3tjzA97natBCcpV3xiTb4z52BizzBjTu5T3JPyclTNXws+ZMaaaMWaBMebDUK6/lfKehH8fy5kr4d/HvY5d0RjzgTFmv40MY36+rLVJcwMqAquA44EqwIfASSXe0wt4PHS/C/B8kuTqAYxJ8Pk6G2gJLD3A6xcBMwEDtAXmJ0mu9kCeh39fdYGWofuHAZ+W8v8x4eesnLkSfs5C5yAjdL8yMB9oW+I9Pr6P5cmV8O/jXsfuA0wq7f9XrM9XsrXA2wArrbWfW2t3ALlApxLv6QQ8F7r/AtDRGGOSIFfCWWvfAL4/yFs6AROs8x5whDGmbhLk8sJa+7W1dlHo/lZgOXBsibcl/JyVM1fChc5BYehh5dCt5KiHhH8fy5nLC2NMPeBi4KkDvCWm5yvZCvixwJd7PV7H/v+Qf3mPtXYX8ANwVBLkAvhD6M/uF4wx9eOcqTzKm9uHM0J/As80xpyc6IOH/nRtgWu97c3rOTtILvBwzkLdAYuBDcAca+0Bz1cCv4/lyQV+vo//AO4Eig/wekzPV7IV8CB7CWhkrT0FmMOvv2Vlf4twazucCjwCTEvkwY0xGcBU4DZr7Y+JPPbBlJHLyzmz1u621p4G1APaGGOaJ+K4ZSlHroR/H40xlwAbrLUF8T7WHslWwL8C9v5NWS/0XKnvMcZUAmoA3/nOZa39zlq7PfTwKSArzpnKozznM+GstT/u+RPYWjsDqGyMqZWIYxtjKuOK5ERr7f+V8hYv56ysXD7PWeiYW4B84IISL/n4PpaZy9P38SzgUmPMalw3awdjzL9LvCem5yvZCvj7QBNjzHHGmCq4Tv7pJd4zHegeut8ZmGdDVwR85irRT3oprh/Tt+nAtaGRFW2BH6y1X/sOZYyps6ffzxjTBvfvMO5f+tAxxwPLrbWjD/C2hJ+z8uTycc6MMbWNMUeE7h8CnAt8UuJtCf8+lieXj++jtfZua209a20jXI2YZ63tVuJtMT1flSL9YDxYa3cZY24GZuNGfjxtrV1mjLkfWGitnY77h/4vY8xK3IWyLkmS61ZjzKXArlCuHvHOZYyZjBudUMsYsw4YjLugg7X2cWAGblTFSqAIuC7emcqZqzPwV2PMLmAb0CUBv4TBtZCuAZaE+k8BBgAN9srm45yVJ5ePc1YXeM4YUxH3C2OKtTbP9/exnLkS/n08kHieL02lFxEJqGTrQhERkXJSARcRCSgVcBGRgFIBFxEJKBVwEZGAUgEXEQkoFXARkYD6f7YeVlabU/11AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pprint import pprint\n", "import scipy.linalg as linalg\n", "\n", "xdata=np.array([1,2,3,4])\n", "ydata=np.array([1,2,5,11])\n", "\n", "def f(x, a0, a1, a2):\n", " return a0 + a1*x + a2*x**2\n", "\n", "def ff(x,i):\n", " return x**i\n", "\n", "Av = np.zeros([4,3])\n", "for i in range(0,3):\n", " for j in range(0,4):\n", " Av[j][i]=ff(xdata[j],i)\n", "\n", "\n", "print(Av)\n", "Ai = linalg.inv(np.dot(np.transpose(Av),Av))\n", "b = np.dot(np.transpose(Av),ydata)\n", "params = np.dot(Ai,b)\n", "print(params)\n", "plt.plot(xdata,ydata, 'o', color='r')\n", "\n", "x =np.linspace(0,4,20)\n", "y = f(x,params[0],params[1],params[2])\n", "plt.plot(x,y, color='b')\n", "\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2 ニュートンの差分商補間(25点)\n", "\n", "2を底とする対数関数($\\log[2](x)$)の$x=2$における値$F(2.0)$をニュートンの差分商補間を用いて求める.\n", "ニュートンの内挿公式は,\n", "$$\n", "\\begin{array}{rc}\n", "F (x )&=F (x _{0})+\n", "(x -x _{0})f _{1}\\lfloor x_0,x_1\\rfloor+\n", "(x -x _{0})(x -x _{1})\n", "f _{2}\\lfloor x_0,x_1,x_2\\rfloor + \\\\\n", "& \\cdots + \\prod_{i=0}^{n-1} (x-x_i) \\, \n", "f_n \\lfloor x_0,x_1,\\cdots,x_n \\rfloor\n", "\\end{array}\n", "$$\n", "である.ここで$f_i \\lfloor\\, \\rfloor$ は次のような関数を意味していて,\n", "$$\n", "\\begin{array}{rc}\n", "f _{1}\\lfloor x_0,x_1\\rfloor &=& \\frac{y_1-y_0}{x_1-x_0} \\\\\n", "f _{2}\\lfloor x_0,x_1,x_2\\rfloor &=& \\frac{f _{1}\\lfloor x_1,x_2\\rfloor-\n", "f _{1}\\lfloor x_0,x_1\\rfloor}{x_2-x_0} \\\\\n", "\\vdots & \\\\\n", "f _{n}\\lfloor x_0,x_1,\\cdots,x_n\\rfloor &=& \\frac{f_{n-1}\\lfloor x_1,x_2\\cdots,x_{n}\\rfloor-\n", "f _{n-1}\\lfloor x_0,x_1,\\cdots,x_{n-1}\\rfloor}{x_n-x_0} \n", "\\end{array}\n", "$$\n", "差分商と呼ばれる.$x_k=1.4, 1.8, 2.2, 2.6$をそれぞれ選ぶと,差分商補間のそれぞれの項は以下の通りとなる.\n", "\n", "$$\n", "\\begin{array}{ccl|lll}\n", "\\hline\n", "k & x_k & y_k=F_0( x_k) &f_1\\lfloor x_k,x_{k+1}\\rfloor & f_2\\lfloor x_k,x_{k+1},x_{k+2}\\rfloor & f_3\\lfloor x_k,x_{k+1},x_{k+2},x_{k+3}\\rfloor \\\\\n", "\\hline\n", "0 & 1.4 & 0.4854268272 & & &\\\\\n", "& & & 0.906425198 & &\\\\ \n", "1 & 1.8 & 0.8479969066 & & [ XXX ] &\\\\\n", "& & & 0.723766544 & & 0.0639712067 \\\\\n", "2 & 2.2 & 1.137503524 & & -0.1515578700 &\\\\ \n", "& & & 0.602520248 & &\\\\ \n", "3 & 2.6 & 1.378511623 & & &\\\\ \n", "\\hline\n", "\\end{array}\n", "$$\n", "それぞれの項は,例えば,\n", "\n", "$$\n", "f_1\\lfloor x_0,x_1 \\rfloor =\\frac{0.8479969066 - 0.4854268272}{1.8-1.4}=0.906425198\n", "$$\n", "で求められる.ニュートンの差分商の一次多項式の値はx=2.0で\n", "\n", "\\begin{align}\n", "F(x)&=F_0(1.4)+(x-x_0)f_1\\lfloor x_0,x_1\\rfloor \\\\\n", "& =0.4854268272+(2.0-1.4)\\times0.906425198\\\\\n", "& =1.029281946\n", "\\end{align}\n", "となる.\n", "\n", "## (1) 差分商補間の表中の開いている箇所[ XXX ]を埋めよ.\n", "## (2) ニュートンの二次多項式\n", "\n", "$$\n", "F (x )=F (x _{0})+(x -x _{0})f _{1}\\lfloor x_0,x_1\\rfloor+(x -x _{0})(x -x _{1})\n", "f _{2}\\lfloor x_0,x_1,x_2\\rfloor\n", "$$\n", "の値を求めよ.\n", "## (3) ニュートンの三次多項式の値を求めよ.\n", "(E.クライツィグ著「数値解析」(培風館,2003), p.31, 例4改)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1)\n", "-0.22832331749999996\n", "(2)\n", "1.029281946\n", "1.00188314784\n", "(3)\n", "1.0003478388792002\n" ] } ], "source": [ "print('(1)')\n", "print((0.723766544 - 0.906425198)/(2.2-1.4))\n", "print('(2)')\n", "print(0.4854268272+(2.0-1.4)*0.906425198)\n", "print(0.4854268272+(2.0-1.4)*0.906425198+(2-1.4)*(2-1.8)*(-0.2283233180))\n", "print('(3)')\n", "print(0.4854268272+(2.0-1.4)*0.906425198+(2-1.4)*(2-1.8)*(-0.2283233180)+(2-1.4)*(2-1.8)*(2-2.2)*0.0639712067)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 3 数値積分(25点)\n", "$$\n", "\\int_1^2 \\frac{1}{x} dx = \\log2\n", "$$\n", "の近似値をシンプソンの公式で求めよ.区間を2,4,8,16等分して片対数プロットで収束の様子を示せ.\n", "ただし$\\log_{e}(2)=0.6931471805599453$である.\n", "\n", "「大学教養数学」,児玉鹿三,技研社 1963, p.172." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.6931471805599453\n" ] } ], "source": [ "import numpy as np\n", "print(np.log(2))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def func(x):\n", " return 1.0/x\n", "\n", "def simpson(N):\n", " x0, xn =1.0, 2.0\n", "\n", " M = int(N/2)\n", " h = (xn-x0)/N\n", " Seven, Sodd = 0.0, 0.0\n", " for i in range(1, 2*M, 2): #rangeの終わりに注意\n", " xi = x0 + i*h\n", " Sodd += func(xi)\n", "# print(\"{0}\".format(i))\n", " for i in range(2, 2*M, 2):\n", " xi = x0 + i*h\n", " Seven += func(xi)\n", "# print(\"{0}\".format(i))\n", "\n", " return h*(func(x0)+4*Sodd+2*Seven+func(xn))/3" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeAklEQVR4nO3de7iVc97H8fe3XUmlHSJNok1JFKO6oqSTQYlGo6FMDmNMjVNUSs6FJMcQk07KqaREEhqjYkyTyqAUQg2l0YQOO+ZpjO/zx28/zzQp7b322vt3r7U+r+va19Vatdb+3Fx99up3/+7vbe6OiIhkvwqxA4iISPlQ4YuI5AgVvohIjlDhi4jkCBW+iEiOqBg7wI+pVauW169fP6XXbt26lWrVqqU3UCQ6luTJluMAHUsSlfY4lixZssHd99vx+UQXfv369Vm8eHFKr503bx7t27dPb6BIdCzJky3HATqWJCrtcZjZ33b2vJZ0RERyhApfRCRHqPBFRHKECl9EJEeo8EVEcoQKX0QkR6jwRURyRHYW/muvUXf69NgpREQSJTsL/5FHaDhqFAwZApr3LyICJPxK25SNHcu6v/+dOkOHwqZNcM89YBY7lYhIVNlZ+BUr8sHAgdQ57DAYORI2b4YxYyAvL3YyEZFosrPwASpUCGVfsybcfDNs2QKPPw6VK8dOJiISRfYWPoRlnKFDIT8fBgwIpT99OlStGjuZiEi5y86Ttjvq3x/GjoWXX4ZOncK6vohIjsmNwge46CKYPBkWLICOHWHDhtiJRETKVbkVvpk1NrPRZjbNzC4ur+/7X84+G557DpYvh7ZtYe3aKDFERGIoVuGb2QQzW29my3Z4vpOZfWBmH5nZ4B97D3df4e6/A84Cjk89cimdeiq89BKsWQNt2sDHH0eLIiJSnor7CX8i0Gn7J8wsD3gQ6AwcAfQ0syPMrKmZzdrha/+i13QFXgBmp+0IUtGuHbz6atiuecIJsGzZ7l8jIpLhzIt5JaqZ1QdmuXuTosetgCHufkrR42sA3H14Md7rBXfvsovf6w30Bqhdu3bzKVOmFCvfjgoLC6levfqP/pmqq1Zx9MCBVPjXv3h3xAi2HH54St+rrBXnWDJFthxLthwH6FiSqLTH0aFDhyXu3uIHv+HuxfoC6gPLtnvcHRi33eNzgVE/8vr2wP3Aw8ClxfmezZs391TNnTu3eH/w44/dCwrcq1d3L+5rylmxjyUDZMuxZMtxuOtYkqi0xwEs9p10arntw3f3ecC88vp+xXbIIfD663DyydC5M0ybBl12+o8PEZGMVppdOmuBets9PrDoucxTty7Mnw9HHglnnAEpLiOJiCRZaQp/EdDQzArMrDLQA5iZjlBmdrqZjdlUnhdI1aoVTuS2bg3nnBNm74iIZJHibsucDCwAGpnZGjP7jbt/B1wGvAysAKa6+3vpCOXuz7t77/z8/HS8XfHVqAEvvhiuxu3TB+66q3y/v4hIGSrWGr6799zF87OJvcUy3apWhWefhV69YODAMIbh5ps1XllEMl52D09LVeXKYQxDjRpw662h9EeODBM4RUQylAp/V/LywsC1/PxwA5XNm2HcOKio/2QikpkS2V5mdjpweoMGDWIHCev4+flw001hvPKTT8Iee8TNJSKSgkSuUUQ7abszZnDjjWFJ55lnoGtX2Lo1dioRkRJLZOEn0hVXwPjx8MorcMopsHFj7EQiIiWiwi+JCy+Ep56CN9+EDh1g/frYiUREik2FX1Ldu8PMmfDBB2Gm/mefxU4kIlIsiSz8KFfalkSnTuF2ievWhZn6K1fGTiQisluJLPxEnbTdlRNOgLlz4Ztvwq+XLo2dSETkRyWy8DNGs2bw2mthz367drBwYexEIiK7pMIvrcaN4U9/gn32gRNPDJ/6RUQSSIWfDgUFYaZ+/fphpv7MtAwNFRFJKxV+utSpE2bqH3UU/OIX4YpcEZEESWThJ36Xzq7suy/88Y9h506vXjB6dOxEIiL/L5GFnxG7dHZlr73CTP1TT4WLL4YRI2InEhEBElr4GW/PPWHGDOjRAwYPhmuvhXAjdxGRaBI5LTMrVKoEjz8eZuoPHx5m6j/wgGbqi0g0KvyylJcX1vHz8+HOO8NM/Uce0Ux9EYlCzVPWzMI6fn4+XH99mKk/ZQpUqRI7mYjkGK0vlAczuO46uP9+eO45OO00KCyMnUpEckwiCz9jt2XuzuWXw6RJ4Wrck0+Gr7+OnUhEckgiCz+jt2XuznnnwbRpsGQJtG8PX3wRO5GI5IhEFn7W69YNZs2Cjz4KkzY//TR2IhHJASr8WE46CebMCXfNatMGPvwwdiIRyXIq/JiOPz6s5//zn+GT/jvvxE4kIllMhR/bMceESZuVK4c1/QULYicSkSylwk+CRo3CTP1ateBnP4NXXomdSESykAo/KQ4+OHzSP/RQ6NIFnn02diIRyTKJLPys3Ye/OwccAPPmhWWe7t3hscdiJxKRLJLIws/qffi7s88+8Ic/hHvknncePPRQ7EQikiUSWfg5b6+94IUXoGtXuPRSDnriidiJRCQLqPCTqkqVcEXur37FIePGhbn6mqkvIqWgaZlJVqkSPPooa7dsoe6IEbBxIzz4YBi7LCJSQir8pKtQgZVXXkndxo3DmOUtW2DixPDDQESkBFT4mcAMbr8dataEa64JpT91qmbqi0iJaA0/kwweHJZ0nn8+3CR9y5bYiUQkg6jwM80ll4T9+a+9FgawffVV7EQikiFU+JmoVy+YPh3++tewX3/dutiJRCQDqPAz1c9/DrNnw6pVYdLm6tWxE4lIwiWy8HN2tEJJnXhiuCr3yy9D6b//fuxEIpJgiSz8nB6tUFKtWoX5O9u2Qdu2YZlHRGQnEln4UkJHHx3GK++5J3ToAG+8ETuRiCSQCj9bNGwYxivXrh1277z8cuxEIpIwKvxsctBBYbvmYYfB6aeHnTwiIkVU+Nmmdu1wn9wWLeCss2DSpNiJRCQhVPjZaO+9Yc4c6NgRLrgAHnggdiIRSQAVfraqXh1mzYJu3aBvX7j1Vo1XFslxKvxstsceYcjauefCDTfAwIEqfZEcpmmZ2a5ixTBOuUYNuPtu2LwZfv97zdQXyUEq/FxQoUJYx8/Ph9tuC6X/2GOaqS+SY1T4ucIMhg0LpX/11VBYCE8/HS7WEpGcoDX8XDNoEIweHQavde4cPu2LSE5Q4eeiPn3giSfCOIYTTwzD10Qk66nwc1XPnjBjBixdGoauff557EQiUsZU+Lns9NPhxRfh00/DeOVVq2InEpEylMjC1zz8ctShA7zyCnz9NbRpA8uXx04kImUkkYWvefjl7Nhjw9C1778PyztLlsROJCJlIJGFLxE0aRLGK1evHj71v/Za7EQikmYqfPmPBg3Czp26deGUU8L6vohkDRW+/LcDDwyf7hs3DjdKf/rp2IlEJE1U+PJD++0Hr74KLVtCjx4wYULsRCKSBip82bmaNcNM/ZNOgt/8BkaOjJ1IREpJhS+7VrUqPPccnHkm9OsHQ4ZovLJIBlPhy4/bYw+YMiXcOWvoUOjfX6UvkqE0LVN2r2JFGD8+zNQfOTIMXBszRjP1RTKMCl+Kp0KFUPY1a8LNN8OWLfD441C5cuxkIlJMKnwpPrOwrJOfDwMGhNKfPj2s9YtI4mkNX0quf38YOxZefhk6dQLNPBLJCCp8Sc1FF8HkybBgAXTsCBs2xE4kIruhwpfUnX02PPtsmLDZti2sXRs7kYj8CBW+lE6XLvDSS7BmTZip/8knsROJyC6o8KX02rULoxg2bQoz9Zcti51IRHZChS/p0aLFf0Yqt2sHixbFzSMiP6DCl/Q58sgwXjk/P5zInTcvdiIR2Y4KX9LrkEPCjVTq1YPOneGFF2InEpEiKnxJv7p1w/LOkUfCGWfAU0/FTiQiqPClrNSqFU7ktm4NPXuGC7VEJCqNVpCyU6NGuE1i9+7Qu3fYxdOiRexUIjmrXD/hm1k1M1tsZqeV5/eViKpWDRdn/fKXMHAgRw0cCH/5S+xUIjmpWIVvZhPMbL2ZLdvh+U5m9oGZfWRmg4vxVlcDU1MJKhmscuUwhuHOO6m+ciW0ahVO6L75ZuxkIjmluJ/wJwKdtn/CzPKAB4HOwBFATzM7wsyamtmsHb72N7OTgOXA+jTml0yRlwdXXcXCyZNh+PBQ9sceG67UXbw4djqRnGBezLsXmVl9YJa7Nyl63AoY4u6nFD2+BsDdh+/i9cOAaoQfDt8C3dz9+538ud5Ab4DatWs3nzJlSgkPKSgsLKR69eopvTZpsvFY8r75hrrPPEO9qVOptGULG1q1YvX551PYqFHsiMWSjf9PskG2HEtpj6NDhw5L3P2HJ8zcvVhfQH1g2XaPuwPjtnt8LjCqGO9zAXBacb5n8+bNPVVz585N+bVJk9XHsmmT+y23uNes6Q7uXbu6v/VWlGwlkdX/TzJYthxLaY8DWOw76dRy35bp7hPdfVZ5f19JqBo14PrrYfXqcHOV+fOhWTPo1g3eeSd2OpGsUprCXwvU2+7xgUXPiZRcfj7ceGMo/ptuCnv4f/rTsKVz6dLY6USyQmkKfxHQ0MwKzKwy0AOYmY5QZna6mY3ZpDsp5Z6aNWHIkFD8N9wAc+bAUUfBWWfBe+9FDieS2Yq7LXMysABoZGZrzOw37v4dcBnwMrACmOruafkb6e7Pu3vv/Pz8dLydZKK99w43S1+9Gq67LlzA1bQp9OgRbrgiIiVWrMJ3957uXsfdK7n7ge4+vuj52e5+mLsf6u7Dyjaq5KR99oFbbw3FP3gwzJoFTZrAOefA++/HTieSUTRLRzLDvvvCbbeF4h80CGbODMPZevWCDz+MnU4kI6jwJbPUqgW33w6rVsGAATBjBjRuDOedBytXxk4nkmiJLHydtJXd2m8/uOOOUPz9+sG0aaH4L7gAPv44djqRREpk4eukrRTb/vvDXXeFm6f37Rtm7zdqBBdeqBuqi+wgkYUvUmIHHAD33BNK/rLL4MknQ/FfdFFY9xcRFb5kmTp1YOTIUPwXXwyPPw4NG4Z5/H/7W+x0IlGp8CU7/eQncP/9YT2/Tx+YNCkU/+9+B59+GjudSBSJLHydtJW0qVsXRo2Cjz4KyzsTJkCDBnDJJbBmTex0IuUqkYWvk7aSdvXqwUMPheK/8EIYNw4OPTSs96/VCCjJDYksfJEyc9BBMHp02LN//vnw8MOh+Pv2hc8/j51OpEyp8CU3HXwwjBkTrtLt1St8+j/0ULjySli3LnY6kTKhwpfcVlAQlnc+/BB69gzr/YccAv37wxdfxE4nklYqfBEIJT9hQhjIdvbZcN994YfBVVfBet2GWbJDIgtfu3QkmgYNYOLEUPzdu8O994biHzQI/vGP2OlESiWRha9dOhJdw4bw6KNh9n63bmF8Q0FBGNG8YUPsdCIpSWThiyRGo0bhat3ly6Fr1zCwraCAgrFj4csvY6cTKREVvkhxHH54mM+zbBl06cJBkyeHT/zXXw9ffRU7nUixqPBFSuKII2DKFBaPHw+dOsGwYaH4b7wRvv46djqRH6XCF0nB1oICmDoV3n0XTjoJbrklFP+QIbBxY+x4IjulwhcpjaZNw81X3n4bOnaEoUND8d98M2iXmSRMIgtf2zIl4xx9NDzzDLz1FrRrBzfdFIr/1lth8+bY6USAhBa+tmVKxjrmGHj2WViyBNq0gRtuCMV/222wZUvsdJLjEln4IhmvWTOYORMWLYJWreC660Lx3347FBbGTic5SoUvUpZatIBZs2DhQmjZEq65JhT/HXfA1q2x00mOUeGLlIeWLWH2bFiwAJo3h6uvDsV/113wzTex00mOUOGLlKfjjoOXXoI33oCf/hQGDgzFf889Kn4pcyp8kRhat4Y5c+D118PWzgEDwjz+kSPh229jp5MspcIXialNG3jlFZg/Hxo3hn79QvHffz/885+x00mWSWThax++5Jy2beHVV2Hu3DCp84orQvGPGqXil7RJZOFrH77krPbtYd48+OMfw01ZLr88zOh/6CH4n/+JnU4yXCILXySnmYUxDa+9Bn/4Q7j/7qWXhk/+o0fDtm2xE0qGUuGLJJUZ/Oxn8Kc/wcsvQ926cPHFofjHjFHxS4mp8EWSzgxOPhn+/Gd48UU44ADo0yfcnGX8ePjXv2InlAyhwhfJFGZhBv9f/gIvvAC1asFFF4WbszzyCHz3XeyEknAqfJFMYwanngpvvgnPPw977w0XXhiKf9IkFb/skgpfJFOZwWmnhQFtzz0HNWrABReEu3I99piKX35AhS+S6czCDdaXLIEZM6BqVTjvPDjySHjiCfj3v2MnlIRQ4YtkCzM444xwE5bp02GPPaBXL2jSBCZPVvGLCl8k61SoAL/4Rbjt4tNPQ8WKcM45YWbPU0/B99/HTiiRJLLwNVpBJA0qVIDu3eGdd0LRm0GPHnDUUeEHgYo/5ySy8DVaQSSNKlSAs86Cd9/9z9LOWWeF8czTp6v4c0giC19EykBeXviEv2xZOJm7bRt0706L3/42nOx1j51QypgKXyTX5OWFNf333oPHHqPCtm1hzb9Zs7C9U8WftVT4IrkqLw969WLRxInhgq3CwrDLp0WLcEGXij/rqPBFcpzn5YV9+ytWhBENGzeGff0tW4YRDir+rKHCF5GgYsVwpe7774ehbF9+Ga7kPe64MLRNxZ/xVPgi8t8qVQqzeT74AMaOhS++CLN7WrUKY5pV/BlLhS8iO1epUpjG+eGH8PDDsG5dmNZ5/PHhxiwq/oyjwheRH1e5MvTuDStXwu9/D599Fubzn3BCuBWjij9jqPBFpHgqV4bf/Q4++ggefBBWrw535GrXLtx8XRJPhS8iJbPHHnDJJaH4H3gAPv443IO3fXuYPz92OvkRKnwRSU2VKnDZZaHw77svnORt3z6U/+uvx04nO6HCF5HSqVIF+vaFTz6Be++F5cuhbduw3PPGG7HTyXZU+CKSHnvuCVdeGYr/7rth6VJo0yac4F2wIHY6QYUvIulWtSr07x+K/847w1z+1q3Dls6FC2Ony2mJLHzNwxfJAtWqwVVXwapVMGJEuAXjcceFi7gWLYqdLiclsvA1D18ki1SrBoMGheIfPjx8ym/ZMoxtWLIkdrqcksjCF5EsVL06DB4c9u8PGwZ//nOYzNm1a7gPr5Q5Fb6IlK+99oJrrw3Ff8stYQtn8+ZhNPPbb0cOl91U+CISR40acP31ofiHDoV58+CYY8LNWN59N3a6rKTCF5G48vPhxhtD8d90U5jPc/TR4QbsS5fGTpdVVPgikgw1a8KQIaH4b7gB5syBo44KN1x/773I4bKDCl9EkmXvveHmm0PxX3dduPlK06bhBuwrVsROl9FU+CKSTPvsA7feGop/8GCYNQuOPDLcgP3992Ony0gqfBFJtn33hdtuC8U/aBDMnBmKv1evcHMWKTYVvohkhlq14PbbwwVcAwbAjBnQuHG4AfvKlbHTZQQVvohklv32gzvuCMXfrx9MmxaK/4ILqLJ2bex0iabCF5HMtP/+cNddYUhb377w1FMce9554Qbsn3wSO10iqfBFJLMdcADccw988glru3WDJ5+ERo3gt78N6/7y/1T4IpId6tTho8suC5/uL74YHnsMGjaEPn3g009jp0sEFb6IZJef/ATuvz/cerFPH5g4ERo0CD8EPvssdrqoVPgikp3q1oVRo8LN1i+6CMaPD8V/6aWwZk3sdFGo8EUku9WrBw89FIr/17+GsWPh0EPh8svh889jpytXKnwRyQ0HHQSjR4c9++efH359yCFwxRWwbl3sdOVChS8iueXgg2HMmHCVbq9e8OCDofj79YO//z12ujKlwheR3FRQAOPGheLv2RMeeCAU/4AB8MUXsdOVCRW+iOS2Qw6BCRPCQLazzoKRI8MPg4EDYf362OnSSoUvIgJhB8/EiaH4u3cPF3MVFMDVV8OGDbHTpUW5Fb6ZtTez181stJm1L6/vKyJSIg0bwqOPwvLl0K0b3Hkn1K8P11wDX34ZO12pFKvwzWyCma03s2U7PN/JzD4ws4/MbPBu3saBQqAKkJubYEUkczRqBI8/Hoq/a1cYMSIU/3XXwVdfxU6XkuJ+wp8IdNr+CTPLAx4EOgNHAD3N7Agza2pms3b42h943d07A1cDQ9N3CCIiZejww8N8nmXLoEsXGD48FP8NN8DXX8dOVyLm7sX7g2b1gVnu3qTocStgiLufUvT4GgB3H76b96kMPOnu3Xfx+72B3gC1a9duPmXKlOIdyQ4KCwupXr16Sq9NGh1L8mTLcYCOpaSqrVrFwZMmsf/8+XxXrRprzjyTNb/8Jd+l8fuW9jg6dOiwxN1b/OA33L1YX0B9YNl2j7sD47Z7fC4w6kde/wvgYeApoH1xvmfz5s09VXPnzk35tUmjY0mebDkOdx1Lyt591/3MM93BPT/ffcgQ940b0/LWpT0OYLHvpFPL7aStuz/j7n3c/Wx3n1de31dEpEw0bRpuvvL229CxIwwZEpZ6brkFNm+OHG7nSlP4a4F62z0+sOg5EZHccfTR8Mwz8NZb0K4d3HhjKP5hw2DLltjp/ktpCn8R0NDMCorW5XsAM9MRysxON7MxmzZtSsfbiYiUvWOOgWefhSVLoE0buP76UPzDhyem+Iu7LXMysABoZGZrzOw37v4dcBnwMrACmOru76UjlLs/7+698/Pz0/F2IiLlp1kzmDkTFi2CVq3g2mvDBVwjRkBhYdRoxSp8d+/p7nXcvZK7H+ju44uen+3uh7n7oe4+rGyjiohkkBYtYNYsWLgQWraEwYND8d95J2zdGiWSRiuIiJSlli1h9mxYsACaN4dBg8L8nrvvhm++KdcoKnwRkfJw3HHw0kvwxhvhRO9VV4Xiv/de+PbbcomQyMLXSVsRyVqtW8OcOfD669CkCfTvH4r/vvvKvPgTWfg6aSsiWa9NG3jlFZg/Hxo3hiuvDLdefOABKmzbVibfMpGFLyKSM9q2hVdfhblzw6TOvn059le/Cks/aabCFxFJgvbtYd48ePVVthYUhPJPMxW+iEhSmEGHDrx7xx2w//5pf/tEFr5O2oqIpF8iC18nbUVE0i+RhS8iIumnwhcRyREqfBGRHKHCFxHJEYksfO3SERFJv0QWvnbpiIikn4X73SaTmf0D+FuKL68FbEhjnJh0LMmTLccBOpYkKu1xHOzu++34ZKILvzTMbLG7t4idIx10LMmTLccBOpYkKqvjSOSSjoiIpJ8KX0QkR2Rz4Y+JHSCNdCzJky3HATqWJCqT48jaNXwREflv2fwJX0REtqPCFxHJEVlX+GZWz8zmmtlyM3vPzK6InSkVZlbFzN40s3eKjmNo7EylZWZ5ZvZXM5sVO0tpmNlqM1tqZm+b2eLYeUrDzGqa2TQze9/MVphZq9iZSsrMGhX9v/i/r81mdmXsXKkys35Ff+eXmdlkM6uStvfOtjV8M6sD1HH3t8xsL2AJcIa7L48crUTMzIBq7l5oZpWAPwFXuPtfIkdLmZn1B1oANdz9tNh5UmVmq4EW7p7xF/iY2STgdXcfZ2aVgaruvjFyrJSZWR6wFjjW3VO9aDMaM6tL+Lt+hLt/a2ZTgdnuPjEd7591n/DdfZ27v1X06y3ACqBu3FQl50Fh0cNKRV8Z+9PZzA4EugDjYmeRwMzygbbAeAB335bJZV/kRODjTCz77VQE9jSzikBV4PN0vXHWFf72zKw+cAywMHKUlBQtgbwNrAf+4O4ZeRxFRgKDgO8j50gHB+aY2RIz6x07TCkUAP8AHilaahtnZtVihyqlHsDk2CFS5e5rgbuAT4F1wCZ3n5Ou98/awjez6sB04Ep33xw7Tyrc/d/u/lPgQKClmTWJHCklZnYasN7dl8TOkiZt3L0Z0Bm41Mzaxg6UoopAM+D37n4MsBUYHDdS6oqWpLoCT8fOkioz2xv4OeGH8U+AambWK13vn5WFX7TmPR14wt2fiZ2ntIr+mT0X6BQ5SqqOB7oWrX1PATqa2eNxI6Wu6FMY7r4emAG0jJsoZWuANdv9y3Ea4QdApuoMvOXuX8QOUgo/A1a5+z/c/V/AM0DrdL151hV+0cnO8cAKd78ndp5Umdl+Zlaz6Nd7AicB70cNlSJ3v8bdD3T3+oR/cr/q7mn71FKezKxa0WYAipY/TgaWxU2VGnf/O/CZmTUqeupEIKM2N+ygJxm8nFPkU+A4M6ta1GUnEs5DpkXFdL1RghwPnAssLVr/BrjW3WfHi5SSOsCkol0HFYCp7p7R2xmzRG1gRvi7SEXgSXd/KW6kUrkceKJoOeQT4NeR86Sk6IfvSUCf2FlKw90Xmtk04C3gO+CvpHHMQtZtyxQRkZ3LuiUdERHZORW+iEiOUOGLiOQIFb6ISI5Q4YuI5AgVvohIjlDhi4jkiP8F4Uq7avgNX1UAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x, y = [], []\n", "for i in range(1,4):\n", " x.append(2**i)\n", " y.append(abs(simpson(2**i)-0.6931471805599453))\n", "plt.plot(x, y, color = 'r')\n", "plt.yscale('log')\n", "plt.grid()\n", "plt.show()" ] }, { "attachments": { "image-3.png": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAMbGlDQ1BJQ0MgUHJvZmlsZQAASImVVwdYU8kWnltSSWgBBKSE3gTpBJASQgsgvQg2QhJIKDEmBBU7uqjg2kUUK7oqothWQOzYlUWx98WCirIu6mJD5U1IQNd95Xvn++beP2fO/Kfcmdx7AND6wJNK81FtAAokhbLEiBDmqPQMJqkDaAAqoAAyMOPx5VJ2fHwMgDJw/7u8uwEQ5f2qs5Lrn/P/VXQFQjkfAGQMxFkCOb8A4uMA4Gv5UlkhAESl3mpSoVSJZ0GsJ4MBQrxCiXNUeLsSZ6nw4X6b5EQOxJcBINN4PFkOAJr3oJ5ZxM+BPJqfIXaVCMQSALSGQRzIF/EEECtjH1ZQMEGJKyG2h/ZSiGE8gJX1HWfO3/izBvl5vJxBrMqrX8ihYrk0nzfl/yzN/5aCfMWAD1s4aCJZZKIyf1jDW3kTopWYBnGXJCs2TllriD+IBaq6A4BSRYrIFJU9asKXc2D9gAHErgJeaDTEJhCHS/JjY9T6rGxxOBdiuFvQyeJCbjLEhhDPF8rDktQ2G2UTEtW+0PpsGYet1p/jyfr9Kn09UOSlsNX8b0RCrpof0ywWJadBTIXYukicGguxJsQu8rykaLXNiGIRJ3bARqZIVMZvDXGiUBIRouLHirJl4Ylq+7IC+UC+2EaRmBurxvsKRcmRqvpgp/i8/vhhLthloYSdMsAjlI+KGchFIAwNU+WOPRdKUpLUPB+khSGJqrU4VZofr7bHLYX5EUq9JcSe8qIk9Vo8tRBuThU/ni0tjE9WxYkX5/Ki4lXx4EtADOCAUMAECjiywASQC8StXQ1d8JdqJhzwgAzkACFwVmsGVqT1z0jgNQkUgz8gEgL54LqQ/lkhKIL6L4Na1dUZZPfPFvWvyANPIS4A0SAf/lb0r5IMeksFT6BG/A/vPDj4MN58OJTz/14/oP2mYUNNjFqjGPDI1BqwJIYRQ4mRxHCiA26MB+L+eAy8BsPhjrNw34E8vtkTnhLaCI8I1wnthNvjxSWyH6IcCdohf7i6Flnf1wK3hZxeeAgeANkhM26AGwNn3BP6YeNB0LMX1HLUcSurwvyB+28ZfPc01HYUVwpKGUIJptj/uFLTUdNrkEVZ6+/ro4o1a7DenMGZH/1zvqu+AN6jf7TE5mP7sbPYCew8dhhrAEzsGNaItWBHlHhwdz3p310D3hL748mDPOJ/+OOpfSorKXetde10/ayaKxROLlQePM4E6RSZOEdUyGTDt4OQyZXwXYYx3V3d3QBQvmtUf19vE/rfIYhByzfdnN8BCDjW19d36Jsu6hgAe33g8T/4TWfPAkBHA4BzB/kKWZFKhysvBPgvoQVPmhEwA1bAHubjDryBPwgGYSAKxIFkkA7GwSqL4D6XgUlgGpgNSkE5WAJWgjVgA9gMtoNdYB9oAIfBCXAGXASXwXVwF+6eDvASdIN3oBdBEBJCRxiIEWKO2CBOiDvCQgKRMCQGSUTSkUwkB5EgCmQaMgcpR5Yha5BNSA2yFzmInEDOI23IbeQh0om8QT6hGEpD9VBT1BYdjrJQNhqNJqNj0Rx0IlqMzkUXoZVoNboTrUdPoBfR62g7+hLtwQCmgRlgFpgzxsI4WByWgWVjMmwGVoZVYNVYHdYEn/NVrB3rwj7iRJyBM3FnuIMj8RScj0/EZ+AL8TX4drweP4VfxR/i3fhXAp1gQnAi+BG4hFGEHMIkQimhgrCVcIBwGp6lDsI7IpFoQLQj+sCzmE7MJU4lLiSuI+4mHie2ER8Te0gkkhHJiRRAiiPxSIWkUtJq0k7SMdIVUgfpA1mDbE52J4eTM8gScgm5gryDfJR8hfyM3EvRpthQ/ChxFAFlCmUxZQuliXKJ0kHppepQ7agB1GRqLnU2tZJaRz1NvUd9q6GhYanhq5GgIdaYpVGpsUfjnMZDjY80XZojjUMbQ1PQFtG20Y7TbtPe0ul0W3owPYNeSF9Er6GfpD+gf9BkaLpocjUFmjM1qzTrNa9ovtKiaNlosbXGaRVrVWjt17qk1aVN0bbV5mjztGdoV2kf1L6p3aPD0HHTidMp0Fmos0PnvM5zXZKurW6YrkB3ru5m3ZO6jxkYw4rBYfAZcxhbGKcZHXpEPTs9rl6uXrneLr1WvW59XX1P/VT9yfpV+kf02w0wA1sDrkG+wWKDfQY3DD4NMR3CHiIcsmBI3ZArQ94bDjUMNhQalhnuNrxu+MmIaRRmlGe01KjB6L4xbuxonGA8yXi98WnjrqF6Q/2H8oeWDd039I4JauJokmgy1WSzSYtJj6mZaYSp1HS16UnTLjMDs2CzXLMVZkfNOs0Z5oHmYvMV5sfMXzD1mWxmPrOSeYrZbWFiEWmhsNhk0WrRa2lnmWJZYrnb8r4V1YpllW21wqrZqtva3Hqk9TTrWus7NhQblo3IZpXNWZv3tna2abbzbBtsn9sZ2nHtiu1q7e7Z0+2D7CfaV9tfcyA6sBzyHNY5XHZEHb0cRY5VjpecUCdvJ7HTOqe2YYRhvsMkw6qH3XSmObOdi5xrnR+6GLjEuJS4NLi8Gm49PGP40uFnh3919XLNd93ietdN1y3KrcStye2Nu6M7373K/ZoH3SPcY6ZHo8drTydPoed6z1teDK+RXvO8mr2+ePt4y7zrvDt9rH0yfdb63GTpseJZC1nnfAm+Ib4zfQ/7fvTz9iv02+f3p7+zf57/Dv/nI+xGCEdsGfE4wDKAF7ApoD2QGZgZuDGwPcgiiBdUHfQo2CpYELw1+BnbgZ3L3sl+FeIaIgs5EPKe48eZzjkeioVGhJaFtobphqWErQl7EG4ZnhNeG94d4RUxNeJ4JCEyOnJp5E2uKZfPreF2R/lETY86FU2LTopeE/0oxjFGFtM0Eh0ZNXL5yHuxNrGS2IY4EMeNWx53P94ufmL8oQRiQnxCVcLTRLfEaYlnkxhJ45N2JL1LDklenHw3xT5FkdKcqpU6JrUm9X1aaNqytPZRw0dNH3Ux3ThdnN6YQcpIzdia0TM6bPTK0R1jvMaUjrkx1m7s5LHnxxmPyx93ZLzWeN74/ZmEzLTMHZmfeXG8al5PFjdrbVY3n8NfxX8pCBasEHQKA4TLhM+yA7KXZT/PCchZntMpChJViLrEHPEa8evcyNwNue/z4vK25fXlp+XvLiAXZBYclOhK8iSnJphNmDyhTeokLZW2T/SbuHJityxatlWOyMfKGwv14Ed9i8Je8ZPiYVFgUVXRh0mpk/ZP1pksmdwyxXHKginPisOLf5mKT+VPbZ5mMW32tIfT2dM3zUBmZM1onmk1c+7MjlkRs7bPps7Om/1biWvJspK/5qTNaZprOnfW3Mc/RfxUW6pZKiu9Oc9/3ob5+Hzx/NYFHgtWL/haJii7UO5aXlH+eSF/4YWf3X6u/LlvUfai1sXei9cvIS6RLLmxNGjp9mU6y4qXPV4+cnn9CuaKshV/rRy/8nyFZ8WGVdRVilXtlTGVjautVy9Z/XmNaM31qpCq3WtN1i5Y+36dYN2V9cHr6zaYbijf8GmjeOOtTRGb6qttqys2EzcXbX66JXXL2V9Yv9RsNd5avvXLNsm29u2J20/V+NTU7DDZsbgWrVXUdu4cs/PyrtBdjXXOdZt2G+wu3wP2KPa82Ju598a+6H3N+1n76361+XXtAcaBsnqkfkp9d4Ooob0xvbHtYNTB5ib/pgOHXA5tO2xxuOqI/pHFR6lH5x7tO1Z8rOe49HjXiZwTj5vHN989OerktVMJp1pPR58+dyb8zMmz7LPHzgWcO3ze7/zBC6wLDRe9L9a3eLUc+M3rtwOt3q31l3wuNV72vdzUNqLt6JWgKyeuhl49c4177eL12OttN1Ju3Lo55mb7LcGt57fzb7++U3Sn9+6se4R7Zfe171c8MHlQ/bvD77vbvduPPAx92PIo6dHdx/zHL5/In3zumPuU/rTimfmzmufuzw93hndefjH6RcdL6cvertI/dP5Y+8r+1a9/Bv/Z0j2qu+O17HXfm4Vvjd5u+8vzr+ae+J4H7wre9b4v+2D0YftH1sezn9I+Peud9Jn0ufKLw5emr9Ff7/UV9PVJeTJe/6cABgeanQ3Am20A0NMBYMC+jTpa1Qv2C6LqX/sR+E9Y1S/2izcAdfD7PaELft3cBGDPFth+QX4t2KvG0wFI9gWoh8fgUIs828NdxUWDfQrhQV/fW9izkZYD8GVJX19vdV/fl80wWNg7HpeoelClEGHPsJH7JasgC/wbUfWn3+X44x0oI/AEP97/BfIdkKvVzisOAAAAOGVYSWZNTQAqAAAACAABh2kABAAAAAEAAAAaAAAAAAACoAIABAAAAAEAAAGKoAMABAAAAAEAAAEGAAAAAEPC2WIAAEAASURBVHgB7V0HvNRU9j6gIEhHEWn6AAHBggqKCCKgKDYQXRUVFVbEjmVde+G/rmXtXXRZBCtWFAVxlWZXyiJKU0DQJwpSlKIID/L/vmTyyJuXmZeZyUySmXN+vzPpN/d+yeTce9oVUVIEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUUgSghUilJlvdZ1l112MYqKiryeXua8jRs3So0aNcrsi+pGvrQlX9rB90jbEr5/kz4T65nMnDlzFdYahO8JZalGHTp0MNKlKVOmpHtp6K7Ll7bkSzv4gmhbQvc30WcSeyT4HM9I9EmunOiA7lcEFAFFQBFQBIiACgp9DxQBRUARUASSIqCCIik8elARUAQUAUVgR4VAEVAEFAG/EdiyZYsUFxfLpk2b/C7a9/Lq1Kkj8+fP973cIAr00pZq1apJ06ZNpUqVKp6rqILCM1R6oiKgCHhFgEKiVq1aQu/DSpXC7Vy5fv16s65e2xbm8ypqC+zWsnr1alOIN2/e3HNTglY9jURNV4K/TlBjvmEPgxeB54APAispAopAyBHgSAJu6qEXEiGH0ffqUWjzuaQ60gtaUIwCEr2ToHEsjrWK8RAsn0hyrh5SBBSBECEQ9pFEiKDKaVXSeS5BC4oPgNCaJCj1xbFnwAb4M3BdcCNwVmgNavL2243kt9+yUrwWqggoAopAJBEIu42iCVD9wYFsMda57yfHPnuVIw6yqX+bOnUqV1Oi2bPryH33HSgjR26W2277WvbZZ11K14ft5A0bNkg6OGg7sodAvjwTIpSsLTSqUl8eBdq6dWtGdT3qqKPk/fffl2XLlsnnn38up512mtnsWbNmyYsvvij33HNPzmDw2haqnqL2bSgCiolsFONxrKsD5UlY7+DYdl1NNzJ72zbDePzxGcZeexlGnTqGMX9++KJIU6lRvkQB50s7+OwKpS3z5s1L5VUN9Nx169b5cn8+2+OPP96XstItxGtb3J4PPqaRjczmCKKZQxo0xfpyx7avq3TOaNt2PXoHAtcxkXPOEUFnQ0kRUAQihsD06dNl//33N422zOW0zz77yNdfl+2PDhw4UC688EI55phjpHXr1lA7v222kr3tQYMGyX777ScHHnigQACY++fOnSuHHHKIHHDAAWbZ3377rbm/Zs2a5vK6666TDz/80Dz+wAMPmD32E044wTy2Bnrtk046ybzu0EMPlTlz6JsjMmzYMPnrX/8q3bt3lxYtWsjDDz9s7g/bT9hVT+MA2KXgMeBOYFoP3NRO2O0f7bmnyCOPiFx5pcC/WmTfff0rW0tSBAoRAXwHyxE1NBdfLPL77yLHHVfusOA7bvKqVSJ/+UvZ4xVplg8++GDp06eP3HTTTfLHH3/IgAED8D8u/0deunSpvPPOO7Jy5Urp0aOHLFq0SB577DHzZl999ZUsWLBAjj76aPnmm29k+PDhcvnll8tZZ50lmzdvRieybC/yrrvuknvvvbdU4DhVO7feeqspdN544w2ZPHkyOqHnyOzZs8378B4URlTVtWnTRi666KKUYhzKIpOdraAFxYtoVnfwrmCOHm4F21Egw7E+AcxXiO6xeJ1kEDgndPrpgh6ACGJTlBQBRSCCCNxyyy1CgcEAs0Q9ddoTKleuLK1atTJ79Pxof/TRR3LZZZeZLd57771lT/QcKSg6d+4st99+u2kDPfnkk81rvMLCMl977TXz9J49e5qxDL/FvGagrpKddtrJ5N12201WrFhhBsR5LTsX5wUtKM6ooJH0drqkgnOycphqKAoJBJjKDzCnY1SopAgoAmkikGwEsPPOAjVN4oJ3RTcy2fFEV1LdQ4M7o8SpTrrjjjtk/HiaPaW0Nx/vKspt2AdcizzzzDOlU6dOZhlUV40YMUL40fdCbmXa96aQsGmHHXaQkpISezM0y8qhqUlIK9Kvn2AIC/9c93cnpLXWaikCisCQIUPgvXibqSq69tprzdEA1T22yocIvfLKK7Jt2zZZvHixLFmyxFT9dOvWTZ5//nkTQI4kvv/+e3M/j9OOMHToUFOtZdsZbKQZiZ7I08tZJlVSu0L61a5d27409MugRxShB+iUUwTGJoGRSgTvj5IioAhEAIFnnnlGdtxxR+EogLaEww47zLQNxI8AaBM49thjZRUMIbRBUE11MQwnNHLTmM0yRo0aZaqFXnrpJXnuuedM+8Huu+8uVG05icZznt++fXvYVgaaNgn7OI3WNJDznJ0xhBo9erR9SJdBIZCueyxd0uji5qSNGw2jbl3DGDDAuTca6/FtiUaty9cyX9rBlhVKW9zcL8s/2WD3nHvuuQZGFIZXl9Jga+vt7l7b4vZ88L2ekeibraqnRMjE9lN/So8LOCuY3hkVnK6HFQFFQBHIOwRUUHh4pP37MwpVZOJEDyfrKYqAIhAJBKhS+ku8320kap77SqqNwgPm3buLTJuGEPGuHk7WUxQBRUARyDMEVFB4eKDwWFNDtgec9BRFQBHITwRU9eTxuWKuD7niCpEPPvB4gZ6mCCgCikCeIKAjCo8PskYNkX//2wrAUzdZj6DpaYqAIpAXCOiIwuNjZJQ2sgkjKlOD7zxCpqcpAopABQgw3uLVV1+t4Kzyh2fMmGEG/vEIA/g++eST8if5uEcFRQpgMhEkUs7LvHkpXKSnKgKKgCLgMwIdO3YszV+lgsJncDMtrndvqwSmIVdSBBSB8CLArLBM6Dd48GAzaywzvnJyoS5dupjJ/L744guz8lxy4iGmE2f09sKFC839binFma6cCfwYec1MtIzUdtJ8pJpmGnKbWAdGYpNmzpwpRxxxhCAY2Exr/tNPP9mnlS4nTZpk1oMR4Uw9/ueff5rHmDKddeN9WT7ThFA4MIU578GIcqY1Z/pzjiyaN29u5rfixQjAk6KiotLt0puluKI2ihQAa4aZMfCeIWlXChfpqYpAoSNAL5BYSm3foMBHUR58MGlxTBnOXE5PPfWUmUX2hRdeMDPDjhs3zkwQyJTfFCYTESBVr149U5DccMMNZpZXt5TiEyZMkMaNG5cmFrSzv9qVaNu2rZl+3M4JRUHC7LRMSshstG+++aY0aNDAFDA33ngjZtIcaV9qJi2kGorCgnNjMA35E088YaYTOR2prFkWM+Hyw1+9evXS6ygEmG6Ec2JcffXVphDh3BZMfsj5L8aMGSOnIA9RFU6wkwGpoEgRPGQLVlIEFIEIIMCeNXvnJE5cdOSRRwoztnIfe+IkfuyZ2+m7774zj/GjTnJLKc7r+DFmgkH25g8//HDzXOcPBcPLL78snMSIH3cyRymcNKlXr17mqcw91ahRI+dl5jmsL4UECelFzHkxWGeeSyFB8pJIkKOou+++2xQUTz/9NJxw4IWTIamgSBNAjiqQ/0tJEVAEKkKggp5/RZene9yZvptzTtjbXLdTed98883mB/+tt94yhQd746REKcWpQuLI4vrrrzcnNIpPDMje/6mnniqcr4JCifNccAIkCqpPP/3ULNvtB5mc3HabKc/tdOSuJ7jspHqNgnAaooQplNwmbHK5LOkuNWYnhaf8QagpoQMUuf/+8sd0jyKgCEQLAY4oqE4iMaWHTbb6yJlSfPny5WbmV86Wx5HFrFmz7NNLly1bthTOKcH05hQaJGao/eWXX0oFBUcttIE4iSowftypLiM9++yzpk2D+3lf2ilItE/YQs7cgR+39OZUXZ1xxhlmxlr7vEyWQQsKmodpPSI617k0pA72vQX+EkxkB4EDJcZTMFEgZjNUUgQUgYgjcM011whTgLMX7pzalCoj9sRpIOasd/zwcmRgz5nNme44zaobUUAwHTnVUKSqVauaLrBUWdEgbRudndcyvTnVRByNUMXFUQ9tD7yWdaGNg9dSfcVJmJx04oknytixY8uUS+P92rVrTWHhPDeK60iMIYvBLcBVwRQG7cBOugEb/4rtaIDlGjDPTUp+phl3S+578cWGUauWYZSUuB0Nz758SWmdL+3gm1EobXFLYx2ef0bZmnhNzV32qnBu2W1h+nSMfBJW0u354KMayjTj9CPjSGIJeDN4DLgv2ElU3NUCVwLXBFNQlIADJXiqYQiIIU7Z0WOgddKbKwKKgCJABDj6oDGd9he/iB/goAizPAhVT4NjFTgby07gS2PbXFBIjAPvHVun0m882I2GYCdZGjZs2IFuYekQ59ilq1kyWr68GqZXPBS5n76Rvn2XJzs10GNe2hJoBT3ePF/aweYWSlvq1Kkje+21l8cnHOxpVDnRrpAP5LUttIXEu/f26NFjJjDo6IZDkH47bkIq3vR/DCo9G9wT3BL8HvhD8DpwPD2FHWRp2rSpYXsvxJ9U0TYDWSq6lg4KcH6Qfv1aI/V464qKDOy4l7YEVrkUbpwv7WCTC6UtDD5jhytVj50UXgvfTqWBmAbhfCAvbYE+ypzylUGGXilIY3YxKtnMUdGmWI/vntN4/TqYAoRqqu/AHF0ESvB6k/vu0/kpAn0IevNQI0Dj7GqkXOZHSSk8CPB58Lnw+aRCQY4o6O/VCtwc/CMY88jJmWAnfY+NI8EcRTQEtwHTphE4YbRq5nzac08GwQReHa2AIhAqBDCql+LiYtMtNFQVc6kMvYhS/XC6FBOKXV7awrby+aRCQQoKGqVpj3gXvAOY8ew0D18IJg0H3wYeBf4KTFXVteBV4MCJqWJo1H4d451+/QKvjlZAEQgVAkwZwUjjKBDVgamoYcLcpmy1JUhBQbwnxNiJPQWETVRFHW1vhGl50EH0jxYE0aigCNNz0booAoqA/wgEaaPwvzU5LHGnnQS9EEHEZA5vqrdSBBQBRSAABFRQZAA6UsKb3k/btmVQiF6qCCgCikDIEVBBkcEDoqBg4N0332RQiF6qCCgCikDIEVBBkcEDOgZRHpwalfNUKCkCioAikK8IBG3MjjSuTCkfl1Y+0u3RyisCioAi4IaAjijcUElhH+Y4l2eeSeECPVURUAQUgYghoIIiwweGbMJy0UU6PWqGMOrlioAiEGIEVFBk+HBo0P79d0HO+gwL0ssVAUVAEQgpAiooMnwwFBQkqqCUFAFFQBHIRwRUUGT4VDkXOrOSq6DIEEi9XBFQBEKLgAqKDB9NZSDIdB6zmQxdSRFQBBSBPERA3WN9eKgvviiyyy4+FKRFKAKKgCIQQgRUUPjwUBo39qEQLUIRUAQUgZAioKonHx7M2rUiV14p8sEHPhSmRSgCioAiEDIEVFD48ECqVxd55BHM08qJWpUUAUVAEcgzBIIWFL2B50Iwpzm9LgG23bGfpmJOajQNHDrirIJt2oh8+WXoqqYVUgQUAUUgYwSCtFFwVrvHwL3AxWDO7DAOPA9sU12sPA6mQPkevBs4lHTAASIffRTKqmmlFAFFQBHICIEgRxSHoOYcSXAO7M3gMeC+YCediQ1MNmoKCe5fyZ8wUvv2qCRE2Zo1Yayd1kkRUAQUgfQR8DKi2BfFf53+LRJe2QRHfnAc5aiik2ObqwhnkyrgqeBa4IfAz4DdaAh2ks1J3Tl3bDq0YcMGSe/aelK//t4yduwcadlyYzq39v2a9Nvie1UyKjBf2kEQtC0ZvQpZuVifiT+wUqHyBfhiMFVBftGpKGiEo7CzsQ6TcBl6FFufgWuAdwV/C24NTkodOnQw0qUpU6akdem2bWldltWL0m1LViuVRuH50g42XduSxguQ5Uv0mVgA46OaMBGRF9VTVxRwFpjT87CgF8C0K2RKHEE4p/xpiu3lcYXynIlgdtFXgT8AQ8kTPqpUKXx10hopAoqAIuAHAl4EBe/DnvxN4GvBR4AfBjNf6sngdInG61bg5uCq4P5gGrOd9CY2DgfvCN4ZTNXUfHAo6c47RU45JZRV00opAoqAIpA2AvwAV0T744RB4OPB74FPBM8CMx75UzCNzelQCS66FPwumB5QI8FzwReCScPBFAocUcwBbwOPAGfDXoJiM6dffxV5+22RLVtgWKFlRUkRUAQUgTxAwIugoJ3g3+AbwH842kw1EUcZmdAEXEx2EgWEk+7BBjn0RM+nzfDf4twU++0X+upqBRUBRUAR8ISAF9UTRwzPgp1C4vJY6dyvFENg332tlbkcFykpAoqAIpAnCHgRFOe4tHWgy76C38Xo7B2gRPs6tMqxgn9ECoAioAikgUAy1dMZKI8BbzQ2O43MjGdYDVaKQ2CnnWDdh3l/t9DGj8dVWDcVAUVAEfCAQDJB8Qmu/wm8K/g+R1nrsU7jspILAi+/7LJTdykCioAiEGEEkgmKZWgXuXOE2xdI1Q1DhMzZ75QUAUVAEYg6Ask+ZR/FGscRxDoH29tRb3tW6j95skidOiL/+19WitdCFQFFQBHIOQLJBEXXWG1ok6jtYHs755WNwg2bNBFZD1GqBu0oPC2toyKgCHhBIJmgsK8/FCsUDjbVxEone0OXZRFo2VKERm0VFGVx0S1FQBGILgJeBMUTaN4GRxN/xzr3KbkgsCOsPm3bqqBwgUZ3KQKKQEQR8CIoKqFtMM2WElNpJDOCl55YqCsMvNMRRaE+fW23IpB/CHj54C9Bs4eC7VEE041zn1ICBPr1E9lzTySngkhVz6cEIOluRUARiAwCXkYUTNJ3GPjHGNM+MQSslAABBt39858qJBLAo7sVAUUgYgh4GVGsRJv6R6xdgVd3HRyK//xTpEGDwKuiFVAEFAFFICMEvIwomuIOY8EUGCvAr4G5TykBAgy2awqEbrstwQm6WxFQBBSBCCHgRVA8jfaMAzcGI0pA3gJzn1ICBDjbXbt2atBOAI/uVgQUgYgh4EVQUHlCwVAS41FYqkIFICQj9XxKho4eUwQUgSgh4EVQrEKDBoCRQNtkrq8G+0G9UchC8CLwdUkKPBjHtoL/kuScUB2ioPjlF+jrqLBTUgQUAUUgwgh4ERR/RftOA/8MZjZZfqy5L1Oi4HkMfCwYihphWnMu44nn/Qv8bvyBMG/vs49VO42nCPNT0ropAoqAFwS8eD19j4L6eCksxXMOwfkcSdgxGWOw3hc8D+yky7BBAzpHFZGhDh1EHn9cpHXryFRZK6oIKAKKgCsCjLpORI/gAPx3EhKD8DIhjkyoehocK+RsLBmjcWlsmwsaz18A9wT/B/w2+FWwGzG2w4zvaNiwYYcxYyh3UqcNGzZIzZo1U78whFfkS1vypR18RbQt4fuj6DOxnkmPHj1mYq2j2xNKNqKY4XaBj/vchFS8YHoQ97sWTPtERfQUTiDDNbWp0b17d66mTFOnTpV0r42/2XffiSxdKtKjR/yR3Gz72Zbc1Nj9LvnSDrZO2+L+jIPcq8+kYvSTCYrRcZfXwPbGuH2ZbBbj4maOAppifbljm6uUbvbQYFesHwcuAb8BDj39C5aVl14SWbNGhC6zSoqAIqAIRBGByh4q3Rnn0G4wP3Zueyyhfc+YpqOEVuDm4KpgRn8zXsNJPFYUY6qcLgZHQkignkKD9q+/wgOALgBKKSGwYoXIjNiYljmzOM/H0UdjyIgx46ZNKRWlJysCikCGCHgRFFT/HAO2XWK/xHq3DO/LyzkyuBRMbyYKoZfBc8EXxhiLaJPt+TSXrVLyhACj2h+BdYzzetx4o3XJli2V5TiMJYsxBr3gAvQu0L14N1I+cJ6aricpAqFFwIugYOV/iGuBF5tB3CWumxOwtzUYnwW5PXbGcCzJ8TQQOziqiAypoEjtUZWg6zBoEFIVw02iG7oijz5qXb/TTtvk3/9GLwICd9IkTLdYG14QcIOYNSu18vVsRUARSA+BZDYKu0QKicPANDRTRYS/cakaCqtKiRDYbTeRXXaxPnCJztH9FgJUL1FIPPecyLBhIjffXD77Lu08PXuKzIRvxvjxIgcdpOgpAopALhDwIiioCnoITFfVYvB/wZeAlSpAgB+2N2BR4dwUSskR4Dzj9BBjenZb5ZToimrVRE45xTr6xRf0JBK55ppEZ+t+RUARyBQBL4ICnzs5K9MbFer1XbsWastTa3edOiJTpiBHzA6pXffMMwjvf0xkjz3gDUF3CCVFQBHwHQEvNopPcFeOIs4D1/W9Bnle4A9Q3D2E8diqVXne0DSbR6+w8/Bm/YhpsTjfeKpuxA88INKlC6I2B8M1b16aldDLFAFFICkCXgQFXVhvAsPZU2g+fBs8AKzkAYElS0SuuMLSq3s4veBOocpo1CgE0CxPr+lVqsBdDv5yNWqInH66NVlUeiXpVYqAIpAIAS+CgtdCEyxXgQ8BrwGPBit5QMD2fNLkgOXB+vxzMb2ZrsKbdfDB5Y973dO4scjIkdb8H08/7fUqPU8RUAS8IuDFRlEbhfUDUwPcEjwWTIGh5AGBXRFPTu8njaUoCxbjJTiaIDa33FL2WDpbxx9veUIdw4gfJUVAEfAVAS+CggF28N2Rf4A/9fXuBVIYRxUqKMo+7LehwPzgA5EnnhCpVavssXS3GJRHot2DsRaVvY6Xrcv0VxFQBBIg4OWv1ALXXglWIZEAxIp2U1AsWIBAFPSilSwEDj3UipegIdtP+vZbkb32EqE3lJIioAj4g4AXQaGftwyxHjbMyveUqkdPhrcN9eUNGojceqsIjdF+ElN/MMXHddeJrFvnZ8laliJQuAh4ERSFi45PLWd09s47+1RYHhRDL7Bp07LTEKqbHn5YhEkF77gjO/fQUhWBQkNABUUOnvhWZMai4fb113Nws5Df4qOPrLiSOXOyV1F6UA2AAzcFRrput9mrnZasCEQPgWTG7EfQnGRqp6HRa24wNWa0MXMYrVwpcvLJwdQhLHe9/34r/5Xfton49lHdx0kOX3tN5LLL4o/qtiKgCKSCQDJBEZsNQBD3Ku3AmILHpFPxi7RsSqkgoJ5PVvT1uHEif/tb9lVxtFXQgYBLJUVAEcgMgWSCwg6qG4hb9ABvid1qOJZM6aGUAgIUFEyVzSypheq2abefc0rkgmwhwRkG69fPxR31HopAfiLgxUaBuFdxerrXxDb3+UGYVUAWgheB4adSjs7CnjkxZs4pzq4XSaKg+P13kWXLIll9XyrNWerOP1+kRQtfivNUyFiEhzZtiheMb5iSIqAIpIWAF0FxF0r+H3hUjGdh6Yc/CTT38hj4WDBVW2fElliU0ndYOwK8P/g28FPgSBIFBXu1hTwtKoXEk0/m9vExXoOjOM5frqQIKALpIeBFUDyNojuBmbqD3Blsq6WwmjYdgivZz1sC3gyG6VH6gp3EUcTa2I7PsETfMJrUGagxg+xhh0Wz/pnWmlHYf/yRaSmpX9+okZWddjTeWGbyVVIEFIHUEfAiKFgqe/+/gPnRbg3uBs6UoIgoM8VqMba5LxGdhwPvJDoY9v0MtivUgDt6e3FmOk5KFATRNZlR8ffeG8Td9Z6KQPQRSGbMtlvHQfvp4LlgDOJNwt9O0EfMiPDpLEcs1416YCcFRVe3g7F9Q7AkS3FxMWY9m8rVlGnDhg1pX1vRzV58sRl05TUxzef8ik715Xg225JKBV97rYls3doKqTW+ALYw1KRIfrTjqKPayIgRDTDX9qdSvbpfU76n2BCc7kdbUr9rdq7Il7bkSzv4lINsy0Lcf6csvGpUYb3rKPd6rJPjifaJxWCOZDxRhw4djHRpypQp6V5a4XVXXWUY1asbxtatFZ7qywnZbEsqFezY0TAOPDCVK8qe60c7vv/eMBYvLltuEFt+tCWIervdM1/aki/t4DPKpC34uNohEeW+s15UT7QhVCl3ZeY7pqOIVuDm4Krg/mB42ZehPbD1Ovhs8DdljkRwgwZt6um/+y6ClU+zyoxlmIHX72w+wQCpWbPt3laanDHAB6G3jiQCXlRP1BXMBk8C/+lo5VDHejqrJbjoUjBHFbSBjARTvXUhmDQcfAsYmZLkcTCJ13Q01yL4Q0FBYspx28ff2pO/v2+9ZcWNnHFG8G2ke/Jpp4lwzgqN1g7+eWgNooOAF0HBXn58T9+vFk5AQWQnUUDYNBgr5LygdnQCBnG2uz59rPXI/P72G5yk/yfy/fciXOcE1w0bWqla2TDmKXGhq6+22rr77i4Hc7yLiRk5VwXTiFx0kdWEHFdBb6cIRBIBL4JidCRbFsJKc4Ie9mbr1Alh5dyqtHq1yPPPW4mqZiJrCwMS3IgNOvJIKxMfp5qrSk2iRfT0atPG3gp+yfQhzLfFQLxTTw2+PloDRSAKCHgRFLQj3Almf7iao1EtHOu66hGBiRM9nhjkaWvXitxzj5XmlfqaDh0ErlpWEAjDquvWhRIQWkCmZp03T+DKhHEhBoZMj8tRBocR6LLf83gNMyJ6OMaIYXEN5kiOar/77hP5y1/CU68gH3cq92ZfYelSa2D58ce7wMtGpF49kb33tpI9plKWnhsdBLwIiqfRnFvBD4DppjoIjH6iUroI2MbUsHw8y7SDH3vqZX75Be4F/a0ZgPbfv8wppRvUJx10kDWSoOD473/xluA1+fvfxUDQwupK/5JFbc+BkAjP60IN2ZWYr/FSWMc++USkS5fS1uhKAgSoaaTKbs89rUSLtq1NZL/SKx55xMKUg1D2Gfr2taajLT1BVyKNgBevp+poIQ3Z/LcvAw8D9wQrpYHAO++I7LabCKfsDBXRHWvgQJFTTrGSI82aJfLCC0iekkBIxFeeNgtOWv3eeyKYdOKP3VvIXT8PlGd/6B66xrKZd2KMzF6wUmIEmG5m6FDrdbCDJVu3ttKwTMIXYfjwmfL555ZgOOkkq5x34ZpyzjlIBtdYhMkfOeBUij4CXgTFJjST5/HTRi+lfmB86pTSQWDXXa1UHvR8Cg0hQFEOPxyJWUbDzwyOZp99JnLAAelXD930e/p+JOfLv2X3VV+JHHggMoWNCs2k4TVqWAMlzjyoVB6BP+HbeNttlp8C1Yb98I/nIJPE/sCQIVakfZs26+WQQ5Cs7VhLmPA4B6EcqdG7jPOW77uv5RodRPoW1kfJHwS8CIorcCv4iwj6FtIBPAB8LlgpDQTatrUuCo2g+AbhKUxExSUni/i///NlIuvXxlaWhYcPlspfzRHhlHODBomcdZbIxo1poJadS6hlo6ZMqSwCfAXYX+jdW2T+fOuDTw2jF2IKfb5OI0daubWghTRznFWrZl1tq129lKXnhAcBL4JiOqoLk5Wg22naJ6CbEHQ5ldJBoGZNkaIiK5Yinet9vYZ+ut26IToGXcgPPxQ58URfit+yBSl/j7CS8Zk5vt9/X+T22zH11UuWUSAkudbfflvkxhtFqFdX2p60kbmxaGd49dXM4n04embWXpZFMxWTMlLgZGu+dH2G2UPAi6DI3t0LtGQaAwMfUdBIwkx9tO4ytWv79r49jSpVRGjcPNced/IeN9wgMn685TLDEQYFU8BEV1mqRKheKWSy53RnZmMO+OjURnWSX2T7MtA/Yv16eMT0QK6e65EyerNfd9Byso2ACopsI+xSPt0yfeq8u5TuYRetlAzooB5gyhTfrbo0YLqGXFCXQesn/SmPOkrk5Zc9VDZ7p1Bgs0oUaptoiStA4ujvzDMtb2gKCkcIjO9ocDQxe7Y10rzrLsvbWieU8h3mrBSogiIrsCYvdOBASxOT/KwsHV23zuouMvc3dQJ0Y/GR6EZJRykaQ12J0Xc0ltMKSssnv9IBEkcVK1ZYDl4BViOQW1PjyKBDymuGzTz2mC/mqaRtoeqVU+LSPsScZ3ffnfR0PRgSBODDUCE1wBnng4vAzvP/im2lNBGgyoND75xGabObT99F2iaoBqIKyGeieyRVGb16JSmYIwrGXNC4Tf9LBu7dcUcg0W8MKOcIL6fPIQk0uTzEfFdvviny6KMil1ySyztbnlQcYdCOQWK/heuFOp+8hUJ4f50f/kS1xKskVCjDIin4BChligDVHLVri9x0EyIZb820tBSupzM8vwwPPmipnlK41OupNBDT7bRTpwquqF5d5JVXrCgt6iGoDmNXkwaOHBL156xGpIjBjbQML1myPfcWR4p8sag7IrPrzun9GNDASDlOWG4bC2KNpWcTzVQc2AVBrBaJ1WY9uP3cc5Zm0jqiv2FBwIug2BmVvTYsFc6HetBVcI89cmzQZhpXSiWOKNiLzwJxJEFtFtM90X5dIfGkxx+3PmisG92P6BnF7H05JhpxOchhzEDoaM0aq3IMUKCNh4p+N0swgxwoRNyIQyYYZYz99pcZVQ+Tg4Z2laYtiyAkICkDpp12svoLfC2pkXzjDbOqAdcq3LenbcnuUzEdDX1TOCo76yx0ErJAXgQF+ohyHBifACW/EGAgUs48n5YutaKemLMpi4mXaHrgN+2EE1JAib1cdm0Zrn7xxSJHHy1CoUb1VA6JsDBFFRPkZhJr6FuVf/zRMpwwtoUCgmpDCtCOHS1Bz7By5t1iN5xYMeMkBQUdFPgVYd4NjtKo1uPIgy8bVI6bR78gB29CYx9BTTnKYFeekp3ODXR3CoD4Clx4IRKC7GclBjj0UJFnnxWxo70DqFKobslHSkFAv5PJky0TX1HRdjfjUaMsOxvTrG3a5KWHlp3mrUexeEsFA0ThOhnj3PBSWGe4c84Sdv31hrHjjobx55/Ovf6um7NdbdliGIcdZhi1axvGkiX+3iCutD/+MIyJEw3jt9/iDnjdfOUVw6ha1TD23dcwfvyx9KpMZu0qLaSClbVrDaNmTcM4++wKTszwcNK28GUYM8Ywevc2jMqV+X0wjAMOMIybbjKMTz81DD7LDOjttw1jBykxrj76S2PrI48axumnG0b9+tZ9dtjBMLp1M5CjyzCWLfN0l6Rt8VRC+ZOKiw3j4IMNg7MilpSUP56NPdloh5/15GOyegCG0aSJ9dgefnj7HZw4ZdIWfNFnhPernoWaRUFQPPec9fC//nr7A/d7zXxpbrnFutELL/hdfHbKmzTJ+mLvuadhLFxo3iOTlz+VSl5+uSW8f/ghlatSO9e1LZRS//qX9RXgF6FZM0s4fPNNaoUnOXvePKuvcNBBhrFxo+NEfmU++sgw2HPZf//tX6TOnQ3jgQcMIwkYrm1xFJ3uKjscP/1kXc1OR9odD48V8L0dFOhr1hjG8uVWh4fSj3Pxkn/+2TDWr084HzL7R7fdZhitW1tFsAlvvmkYTzxh/R22bUveqEzakkxQeFE98VPeB9yNK6CpYKqj/KDeKOQhMMdLI8B3gZ1EBSqPU/X1O3ggeBY48tS1q5U+Ipv5hurMmSNCAzYj37I8xRxtq08+aeUEokYjbaIqZOpUy4WXqV2ZRTFHdMUVlrcuPXYZUZx1Yjp3GvJpp2G+brpgPfWUFdzho/sPtVann445AqpZ83CUMQHRTkScyfQ8Y2ADrfv0mWWaXTKP0Y+W7mEZPVxviLKe9kRXf/2rpTWjD4bPntzeKuM8i77f1AGRmR+Naj1bvffzz9CzQNFiOxU4r0u0TocOJB7DqE5+q9JAFqxpIHN+hgrWaCCXtmoom15sCq+QPaRP5z0sl7A4Z4RExWZjvxdBgTdZ6Ef5fKwCl2OJz5xcF9tOd4E3VB4D9wIDdZkOhkJW5oFtOhYrrWLcCcsnwFxGnqha5ocpawQddVv+8anHzkGsAr1tmaXDl7mxaUv5+GPLXoEw3rrDhol07541qOyCqfflt/Crr6yBftb+lwxgYNAChTg/Pox4o4EkS8YRyhw6lNHOTSeKpLTXXlbYNEOnmf+LQoMOBnxZyRQazPjHLMM5IKaDp4yiJzeTGdOcknWifefLLy2DFTtbCxdaWDC03Em0C9meZbQd0WGA+2ymIOBLxAdgv0x89vSccPD6ZWtk1ju/SKMdvpUzqn8iNTetkkrfQrpf4rgZpScnfucD5JIvK7l5c4vp3ebJg8RRps+rQErQ0lLiB577MqXOKOBdRyHXY53spCexcYZjx0KsN3Jsu66mpXrieHzYMGP23XcnH9v5eJSq4OnTfSzQWdR55xnbqOf+/HPn3qyt9+tnGNQWVTQ0TqkCHIfvt5+xtUoVw3j55ZQuTffk339P90oP123dasylvaGoyFLx0Bbx5ZceLkz/FN+0VwsWWDoRPA9TYV6pkrGW61SWO+xJ6dc08ZVLlxrGgQcaBm5p3H67z+8YdFxz7rjDMP7xD8M4+WTDaN58u/qN/YWGDQ3jiCMM4/zzDYPfhjfeMIy5cw1j3brEFa7gCF4D4/XXzc9N6ZkTJhjGpk2xTZ6wcqVhzJxpGGPHWhhffbVhnHaaYVAl2LixBYbVn7Hqy/9Iy5bGLzyeJuFjmtBGsaPrl7b8zrrYBX8WkyA2faEmKOUHR0nFWI8fLbidw31w5yhHQ7CHjFFhMbQXU7nqmSrBt7MTVABNMOadmoVANLeK3HzzPrJsWQ2kY/7C7XDa++pNny7t//MfWYxeXzFnqEsRi1RvvHVrJbiWdkGnfyUSvqEX6iPtiFFR22uvlfrQnXwL75/lnBEnB7R2bRXEupSgk4bPog9UFy6tLeBa1Q690/XotS/BxE5rOXKim1iWns+cOXWgOTpArrlmAZyaVmTeCupLwTtj3vQGqPMudMOBT6tx+eXyG1yWVuHYGvi3/s5er92DzvyuZgl33FEZ0eNt5N5760m7dtPhoLUl5ZJ3gHqvFvCvvWCBuayFZTWMEuBsZdLvTZvKBjybDVABbmjVylzfDLVQOaIfKjlFovv4tGkNECuyJ6LSa2Jg8Dsy7c5A2Ms24eDj009dCqQnGpkuYQ6qBPfoaqhDNai8qkH9ZS6xvhUjllS/fY5iM1plj34ZeBR4NPg7cH9wpoQBpWmXsMs5GyuP2Bux5XgsqeayaRJW8O9KTmmNKCiF77rLks5ffZWmTE7tshtvhBcKnE1KexKpXe5+Nns6e+xhGG3aGNPefdf9HJ/30hbKzg2dlrJB0+hKdeKJ1k0w6vN32FK+xrNnG8ZOOxnGiy+WP5byHnorHH+8VXcYqefRaMweY5bpl18s2/hee2XPGDxlyhTDmD/f6o3TU83u4dIYzx74q69axluf2srRqm1bpw3eXnctnhZwvpgPPmgYZ51lWYft+nGJ3rfRv79h3H+/MeuhhzIaIbjeP27nF1+Yf0kTorZtDYPOLBk6scXdwdo0n4nrkYp34quacESR/Iu7/SjVPTRoszu3+/bdGa11xtXvOkq4HutkJz2JjTMcOxZivZFj23U1bUGxapVRQvfMCy6oGFUfzqAjEt9ZX7UPF11kDUs//tjI5KVJpXl86evU2e6lkcq1Xs4128F/1cCBFmADBsS57ngpxfs5W/Edb9fOMPbZJ4NvOtUxgwdbbq4Eh15NcOfJxTNh/Smb+CrPmuW93ameWa4t1BE9+aRhUA9Jd2z7w0x1zplnWiqUyZMt4ZGhjpKaonr1DGP8C79ajaRbMXteffpsV+3Z96eqpm9fw/jnPw2DnafVq8s0tVw7yhxNf4PezrYw4+vQqZPVmeLzyRZl0hZ8TBMKimSqp71x4QLwQbGvMVVDJFhNTM7U+2g6yqGhGtYY+RHMUcqZYCeNwwbMWTIGTLUUrEyuaifs9oHggrQSWU0bMdqHc2VmOeiLQXckxkJ5nXHUuiLBL1UBT8DeT08VpgLNkkoj/u5M2USvGsZ7ZY1YOGfDadnSCtBjvipmlqMxz2ei7ZHpVWhjHjs2Rbst82gzwx7DZRn4xnBjFpZN97a49t9/v5XKizmcOLlgzogeGpz+jsy2QwVq6lOoU+G7SGu0TfxvQb1jphihQZhMFQt1MHTLYrg2y2AEOpm40pi8apW5/PviYrls3RKpe+Yau0TLmEvXKOaPOf98S13DhFI58NTaXgkrJQlfVXrOFRVZgXG0NTMgNaqU7K99FRqFJy5448sROgvSs9ze1HaU4HQKgXfBNJADWsEnUy4Ek4aDJ4CPAy8CQ9luTpyERfaoGDkcGjEPBZ80U4tmkfhO01HBlwhtelEMHixCrxV60+SI2G2kOjqrQsJuC2/Ejy51+/yKc8npW7OQs52OPcOGWVlwTz7Zg8qdHzW6FnF6OOqvKTnpBkbBlmNihC4ztTDQPTBifgl2Vsgkvih0J2UOek6bR1682HLH5XwotNVURPyzMHNggwZSDYKlyqGnybtLWsiIKS3kt11ayk3P7S3djq5WUSlZO05Zxgh/Cmp6y3KmPzqP2f+RrN04JAW7Ie+2LyTV5fejQ9ojO3PoxghVeqY4Qx7TLjH5hePHGwZH7BnT0KHWUH/atNKiMhmGlhZSwQrV0IzT8hjMW0Fp7odd27FokWG0b2+1mSqeDLxQ3O9qGKNHW1q8pJ5pW6FHoDGDOm9+E/jufPZZoiIN17YkPDvcB3xtCw11UP2aQWkMtJwzx7J/LF5s6W8YwEasXYjPh+/gjBnWwVS1Wn61g/GJfAWOOsowqGFLtR4uTUt5VyZtwQd8RiYf8VkuF7vtczktmF0ZCwp+/fjEGRIZBfrwQ+uLdumlZWqbyUtTpqAkG7RZUh2dDcOcfduE7eDH5dprrba3aGHlD7Ev8mHJNtEr1JX4FaCUp98m3xW6inK7gq9Dwra43sT7Tt6WNtuRI71fk+mZ2WpLOvXa6pAhQ4YYBrzDDfYlvFA67SDe7A/QXMYOBYkB1znyRrdu6PKbTlvsYvC1TigonPER8V/13bEDY3uB0lCo6YSyz+TuWEKJmMdEF0wGtTzwQNYbuXSpNVzl/BRpES9k+GpRkWVXSauQ9C7iF5LzTzCgOCeqp/hqUo/NyGaqLqiW4HR1zCRHlYYPxDZxniUSvYxNYtTamDGW8p/RX1SZPPOMFZx1HLSkVI8FQIzfe/55EQZ7FyLRrkRiFDpj05iunKpdaiinTqUkNw9n/MMMBPBsNuf+ZvJCRozbcXjM7M7st4VG56LBU8DQvJlLrpPHgaG1DS9lPKKgiL3vPqunmESNYEviTJYvvWTdJm3vFAbi8H/w/vvlqpFJ76JcYS472NvmrYcPdzno4y5P7eDo4s47DaNGDcvn+NxzkwwHUqvc3/5mGEe2W25s/eft2wOy9t7bMEaNMozNm1MqzFNbUirRcvqhhxM9nSoY0KRYcvLTs9GW5Hf0fpRplvjc6HDGd5QxdSQ+rvhHlqwd9FxyeiV26WKVx6SFjz2WFY2nVdE0f5O1paIi8UVPOKLw8rU/xctJYTrHF0FBnTf97+jql0Wimz1f5GefTeMmHOcy+ppjbRfK5KVxKa7cLgblsu5ZTkprpNQO+iFecYVhVK9uVa5nT8tpPZ3McitWGMaIEUbxgccbW2QHq7wePayw2q0OXUc5ZBLvSKktiYspPcLXtFUrK1iXsRO5JL/bko26M+EC/1v2O8rAasbI0Iw5aJBh3HyzYQwdutBgXkYSQy+ozaRKiSYwBjwz1ulXeOGSGA/hW7S7VaSvv5k8k2SCAoPrhDQARzCAkyLwVeB4gm0/j4n5WphohpM/00OjbdusNJYeglRxpOz5xJwxVDnRrTCgiYeplqFnTRY8VNPHmn6IVBnS3YQuKKNGiQzAq0yQ6YbCXEX0S6Z3GF1WqS+g+yW9xqhX4ETOs2CCo0snXXAhHhrD7fOZ3a6Sp3cYLBMntDZVG+lX0N8rOXUHp5uYhFBUe1pRf+8Q7dLoacvHbxM1ypwClvOOMN8k50tHNJ75HtM7dwb61PRa4t+qXTtLm8lEDVWrWiXkKGmDXd3QLPHvSUg1YkfwTypQ4htFhSQ/xE8/nRUQ+AJSl5qyoKDrJS/i3KMBTfjMOYbIoSROhMQJkehOywSDEydaBhU+z0SzwNkNIZ5UQDMbHVxvK7VvL00mVZJpvaxEr1e5dZvsa3O8pA6eHy92OJQqRoBhFWSb+CqMG/cxcmeiAwFix4ehLwGZmuxqhW6ZTFA8GastHMMLlOCvbQbuMA00/eMrTL2ZHk7s4LIn45mY2ZIBgUzVmpN0muVrxgnYGDpQVFT+WKj20Mp5+OEWU7hy9MA00UuXWpZfpoWmUZyBXsj1Yw6PGKBlW0djjUEcpikUaTRGaqOgE3WaCU1puKY8UyGR/hvHgWb9+ltKpxW1pxdNv8T8vDKZoLBbjO60/BNMvxx0y6Q9+Arwc+D8JwbdMdqZH2Yus0Ac6lID4on4dR40yFKbPPigp0uycRKjftk558fKc92zUZFUy+QQDnNHm5zitXz8mD4gcCHBqSsYBEjsqXail4+SIpBNBGJOZUlvQeUCul1yArgY3Br8d3BhEEcRTAcwYoRvbpfxwLED61l7xPQQVLCya+uW2TK+8Cxt//e/VtBtpIREhlhwag9GPdMFk7rtIIgW9fPOQ26dBda80iokgngKhXdPL4KiSgyW47B8Ebym4GCinptj0ltvzUrT2UO8+Wb4Hk+poHgKiGHDLN15jiaOcasR/cZp7w2tfcKt0j7u69/fmoCPaaNzTRx9cvI5ho8wfkVJEcgFAl4EBfwqzOSAHbGcBIbiXjaBC4foAkELF5Oacfozn4m9QtrLk876ycA6Zt+ja0uWVGBem/X++9aZhSooOAseZXYO4jHLPBIK52uusWbh44R4SopArhDwIiiuQ2XgVygUFFCQy0ZwX3BhEf+hdJm9jnD4SzSoUW1OG3VC4n3ppjtqVE4zkbrVh2onar2c3iNu5+XrPjpDIXek3HhjVvoNCWHjTKkcUTBfpXrlJIRJD2QBAS+Cgqqns8EvgV8FQ0Mqq8GFRfwy0t2SmWXpkuozMR00e4zUQZej994Tefhha1QTgm487fpMv82sGYVI/Eg/CZ9AZsqmjz5DWrJJtIcsXWo5YtHjiv0VJUUglwh4ERTw9TBzPsFHVMgHgbmv8IjqJwbe8d+6yV/tG3vnTLVPt9MyhGkOTTdYRv9QMR0Cwmyx0q1bCCoSYBXoOU3/BqaW5iPKFjEOsE8fyx5Bz14lRSAIBLwIioNRsXPBk2MM30zhvsIjGrTZs6dPIn1DfSSOKNhTZM+xlOgKy3kN+DV65RXL17/0YDArjAQmBPT8KXQ6AX6AjHnMViwJExHyHna0sB0dXOi4a/tzj4AXQUHfjpaOqrXAegD+Ho4aBLnKyCvOasPUHj4athk49euv5tz121t3ww0iH35oTYjDEUUIiFkx6JkbF48WgpoFUwXG6bGnz8HmtGn+1YG+C0xizDKZnJbrSopAUAh4ERR/R+XouDkVzL8CRxZ/A2dC9XExFO/ybWwJbW85aoY9vO98MPptAn1PSIjRZkwMw2nEfNIH8MNb5uNLDyuOWphTgHkaQkDUxU+dKtKrVwgqE6Iq0L2ZnmDs/X/xhT8Voyc28zcxcwyd3ZQUgSAR8CIo8Lqac1ujzyTkNmB+wDOh63CxXS6X3I6nEuygQIJRQNDflkvA4ehWU0FNa+bs2ZaBGxXzg6jz5pQG5iiC0ddHHGG5ufhRuA9lME8e1SEhsKf70Br/iqCfA/0NmF7qmGMyG1nYzgwM3Rk/Hjrfc/2rp5akCKSLgBdBwQQB/EgPA98CvgjMfZkQB9KjYwVweZJLYT9h36zYfijpzZFFE5fzgtnFCXIYsc0Z1F97zZc6/PYbkpe+M1+29UXZTMlK1yLmIQoJ8WNIT6fu3UNSoRBVg9H1kzHWZsgNR1wvvZRa5SggOP03+wYcudWubQX1pVaKnq0IZAeBSh6KfRnn8EP9XOzcM7CkqujU2HY6i19xUV3HhWux7qZ+sk8pwsoH4H3B68BuNAQ7yUiz0LDDGM5ClgZtgB6hpse8FJWgdjrgyiulJozbsx55RDYydXUGtPSdNdL77oulVs0S+erJh2UTU2ZnQKm0xctt7r67jSxfXl0efHC2l9N9O8fvdvhWMZeC1q/fUf7xj3Zwm10m7dv/Zro7O2Me3NqyaFENefzxvRDEVw/X/Ir8k3OR0mWLS+nh2uXWlnDV0Ftt8qUdbG0mbenRo8dMFMF4ubToS5er3PbFnwatrXztwhxNUFA4iYIiEdXEATbg5EQnxO/3ZeIir9OJcKKcJk0MY7fdrMngvV4Xf97s2cbW3RoaP0lD499XzYs/mtZ2JpOYJLphNufGTnTPbLQj0b382O+cZW7oUMPo398wxo41DM66NnnylNJZ6Djp0FFHcSxhzZHFGdNKSvyoQW7KiNpzSYRKvrSD7cukLfiOzoj/ltrbXlRP/8PJtBHY1AkrH9sbSZZwDzJHABwFOPlNbK8AY5BuEpcrY+vxiyrYQb3O8+DX4w+GYpu9fuoc2G3s2dOKmku1Ygzg69pVKlfZUc5uPFneK26bagk5O59R5ErJEXCOIOgVxUh2RnLzVTnyyCNMj2eWwIEr1Ux33CGyaJHlt1CoQYzJEdWjQSPgRVBQMHwCXhrjT7GEJlW+As8Bp0PjcNG5sQu5pPCIJ6rF/gOm19P98QdDtc2Zh+iiwjgLfPBNf0bbKpmsovSBZPAeJscRThcHl5l2f2kn1HeHjTiHE3McKaWGAGMkf/7Z8hZ76CFGci8rTeZHgfIBFKqcjC/ARMCpNUjPLkgEvPQPe2cBGfx9hLaP88Dfg217B5XyI8D0/ekCPhtMgWQrxW/A+gRw+IjJmugbya8pXVXo3jpsmDWzTHxtGW77PAZJjMUoLrbmZuQXBfM28mMSNqLMGwfR3rFj2GoWjfqw/0AjNXnq1KVwBiiKRsW1lopADAEvgmJZFtBajTKPdCl3OfZRSJA+AnNUER3iZAVTp4owzoIz4nGOZo4UOMMa9Q5M+zFvnnUOne87YbBGgRGXD4MfZqawDouahxPCfQ9xzp6vkiKgCBQeAl4EReGhkkmLqWSmOomzyzCkll3xN6FZ4yQO7Fpy3koG0DGKigLEqdDGfRmnQE0WI32ZsDYMRLdYksZPWDjoryJQaAiooMjWE6elklHVZBKHCB4sldA+mVNbfv65dVkYfmmM5exuZCVFQBEoPARUUOTqmXsQEnZVqJHyM2+QXW66SzpzkZUUAUWgMBFQQRHC505BQVs4U46HwQOKmjQlRUARKFwEKhdu08PbcgoKUhjUT0yjvW6dVR/9VQQUgcJEQAVFCJ87p7ykMbuoKPjK0e5+8snB10NroAgoAsEhoKqn4LBPeGfmAQxDPAVDPObMEbn77oRV1QOKgCJQAAjoiCKkD5mT23EO7ZKS4Co4caJ172OPDa4OemdFQBEIHgEVFME/A9cavPwyJirvYMXnuZ6Qg50TEAPfrJkIg86VFAFFoHARUEER0md/yCFWxThZUBDEifs4axtHE3ExgUFUR++pCCgCASKgNooAwU92a05twRnTPvpI5IILkp2ZnWMMIqeQ4lJJEVAEChsBFRQhff7sxTMF1IcfBlNB3l9VTsFgr3dVBMKGgKqewvZEHPVhKqhly6yEfI7dOVm94QYrd2FObqY3UQQUgVAjoIIixI+H8QvMs0QVVC6JwunOOzGt4Mxc3lXvpQgoAmFFQFVPYX0yqFfTphbnuopMdkvq08da6q8ioAgUNgJBjSjqA3Ymr/42tqyX5DEgb7dwOta3k5yTt4f+h5Y/+GBum/fGG5Z9ghnRlRQBRUARCEpQXAfoJ4H5KeKS24mIKenmJzqY7/uperrySkwqvjI3LV2NKaU4PedJJ+XmfnoXRUARCD8CQQmKvoBmdAweLhN9lqB8kePBI2LnFtzCnvwuV95P330nssceKigK7kXTBisCSRColORYNg/9isLrOm6wFutu6qdXsR9mVakFvhp8AjgRDcEBsjRs2LDDmDFjEp2XdP8GTFFak5MOhYRKSirBVtBFevVagZEFNXXeKd22cCpWUlgC7dJth9WKcP1qW8L1PFgbfSbWM+nRowfdVzpaW7n7RVyvfO3CHE1QUDiJgiKeKBQej+3sjqVnG0WHDh2MdGnKlCnpXpq16/r0MYwWLVIvPtW2bNliGCUlqd8n21ek2o5s1yeT8rUtmaCXnWv1mVi44hs7I/a9LbfIpurpKNxtXxemT80KcCMwiUs3DXxuZtlfAAALFUlEQVQX7KffzVIwhwc9wc+BC444VzWn3F5B1LJI77zD0ZgI56BQUgQUAUXARiCbgsK+h9tyHHaeGzvAJYVHPF2PHbRRFIH7gyeDB4ALjgYNEqGRmR/xbBITEW7bBg8D9XbKJsxatiIQOQSCEhR3AaleYCrdueQ2qTF4grmmP6UI7Lxz9nMu/fGHCN1iTzlFpGrV0lvriiKgCCgCEpSgQP9YjgSz78rlGjBpOfg4c63sz1Rs0mZRsPTqqyKdO4twnopsEFOKw44v/Tl2U1IEFAFFwIFAUILCUQVd9YIAPZA++0zkk0+8nJ36OXQSo2qre/fUr9UrFAFFIL8RUEERkedLgzanSKV6KBt0+eXW9Ks7MA5eSRFQBBQBBwIqKBxghHm1FiJJesGaM3asiB3n4Gd9u3YVOf10P0vUshQBRSBfEFBBEaEn2a+flXZ89mx/K/3AAyJz5vhbppamCCgC+YOACooIPUtmcz31VBE/1UOLF4tcdZXIW29FCAitqiKgCOQUAU0znlO4M7vZrruKMNbBTxo5UqQyugsDB/pZqpalCCgC+YSAjigi+DQXLRL54YfMK15SIjJqlEjv3iJNmmRenpagCCgC+YmACoqIPdf165EXBYlR7r0384pT3bQckSuDB2delpagCCgC+YuACoqIPVt6P514osgLL4hs3pxZ5Skk9tvPKi+zkvRqRUARyGcEVFBE8Omei+xYq1aJjB+fWeUvuUSEHlQ7qqUqMyD1akUgzxFQQRHBB0ybQrNmIg8/nH7lv/zSisegIVtJEVAEFIFkCOhnIhk6IT3GEcDQoVZKjx9/TL2S8+aJHHigyKOPpn6tXqEIKAKFh4AKiog+8wsuEPn++/S8lYYNE6lRQ+SMMyLaeK22IqAI5BQBFRQ5hdu/m9Go3aCBpT5i1levNH26yCuviFxxhQjjMpQUAUVAEagIATVjVoRQiI8z5xOjtWlneNNt6qe4um/dKnLRRZhSEHMKXn113EHdVAQUAUUgAQI6okgATBR2M/V4F0wYOw7zBZIrouJiK26CuZ3q1KnobD2uCCgCioCFQFCCoj5u/x7429iynlWdcr91sQdT9sgC8Hwwpu5RciLAPE377y8yZAgmHl/pPFJ+fc89rfmwTzut/DHdowgoAopAIgSCEhTXoUKTwJzhjktuu9FD2DkRvDe4PZjCQsmBAKctffZZkXXrrMA5N3vF3LkAGAhzPux6EMkciSgpAoqAIuAVgaAERV9UcHSsklye5FLh2tjXDfyf2DHGIf8aW9eFAwGOKDhD3c8/izDFh020STzzjMhhhwFsoMzjSoqAIqAIpIpAUH1LfvCpVrJpLVbi1U8HYN9TYHj9m6OJmVheDt4IdiMoX4SMKT0bdhjDL2catAFd8po1a6ZxZfCXbNpUWapV2yYlJZVk0KCDZe3aKrJxYxVp23ad3HLLXNl99z+Dr2QaNYjyM4lvrrYlHpHgt/WZWM+gR48e/MZ2zPUTeR83/NqFOZqIHxlQUMQTK4z8ptIpdoBqqNti60kXHTp0MNKlKVOmpHtpaK5budIwzjnHMPr0KTbeeMMwtm4NTdXSqkg+PBO74doWG4nwLPWZWM8CH9UZiT6s2XSPPSrRTbF/BRhOmvJTbOlmhoWPjpA/B5No1E5kyzBP0B8LAcZXUNU0deq30r275g/X90IRUAQyQyAoGwWdOZHaziQu34ytOxfUqHPWhTaxnUdiSTWUkiKgCCgCikAOEQhKUNyFNvYC0z2WS26TGoMnmGvWz2VYPA+eA6bN4g6wkiKgCCgCikAOEcim6ilZM1bjIEcI8bQcO45z7JyN9ZwbVxz311VFQBFQBAoegaBGFAUPvAKgCCgCikBUEFBBEZUnpfVUBBQBRSAgBFRQBAS83lYRUAQUgaggoIIiKk9K66kIKAKKQEAIqKAICHi9rSKgCCgCUUEgqBQe2cbnF9xgWZo32RXXrUrz2rBdli9tyZd28P3QtoTtX6LPxH4iyC8tCNdV8oJAwjB2LxeH7Jx8aUu+tIOvh7YlZH8SfSYVPxBVPVWMkZ6hCCgCikBBI6CCoqAfvzZeEVAEFIGKEdih4lMK8gym280Xype25Es7+F5pW8L379JnEr5nojVSBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQKCYHeaOxC8CJw1CZIGok6c/Knr8E21cfKe2CmcucyfqpZ7AolNUOtpoDng+eCOf0tKWrtqYY6fwH+Esx2/B+YFLV2WLW2fmnT/B/47djOqLZlKer/FZjZqW135ai2pS7a8Cp4AZj/mc7gqLYFVQ838Q+wGNwCXBXMP3c7cFSoGyp6ENgpKO7Gti3wuPwXOArEmQ/ZFlIt8DdgPouotYfBrDXBpCpgztR4KDhq7UCVS+kqrL0AtgVFVNuyFG3YFeykqLZlNBoxONYQfrsoOKLallgzwrugFH7XUb3rsU6OEhWhsk5BsRDbjWIN4JLbUSTOfsjJraLcnp1R/1lgzv8e1XY0Rd0ngXuCbUER1bYsRRviBUUU21Ib7fgOzE6Jk3xvi8ZRWPA2weIHB9Kcq5v7okwNUXnOSU7icjdzLVo/RajugWD2xqPYHo5Uqd5YCab6L6rtQNXlQfA14G3ciFEUnwmrboD/C6ZL7BAwKYptoQaE6YqeBlMlOAJcA+x7W1RQAFVQvETmPr5MSsEhQLXNa+ArwOuCq0ZGd96Kqw8Aszd+CHhfcBTpBFSawi5fYg26oC1Ubx4LvgRM1W0UiTOUsh1PgNmh2gi21c1Y9Y9UUFhYcgRBI6pN/GMvtzciulyBejeK1Z1L/tGjQtTpU0hwvvTXY5WOcnt+RRumgukwEcV28MPaB7wUPAbcE/wcOIptQbVL/9v8T4wFU4hHsS38bpE/B5NeBVNw+N4WFRSEV2Q6uBW4OZgGof7gceAoE+t/bqwBXFLXHwXi6O4/4Png+x0Vjlp7mIWzbqz+1bE8CrwAHLV2sAnXg9l5KgLzvzEZPAAcxbZQNVMLTOL60WDa9qLYlp9Rb6rM24BJR4LngaPYFtY/EnQcakkPm8XgGyNR4+2VfBGrtENsAbOHcR54F/Ak8LexZX0so0BdUUmq/eaAZ8eYzyZq7dkfdabemO3gh+gWMClq7bBqvf23O1ZtY3YU29IC9f8yxnRbtv/rUWwLqm+qNmdgyffsDXA9cFTbgqorKQKKgCKgCCgCioAioAgoAoqAIqAIKAKKgCKgCCgCioAioAgoAoqAIqAIKAKKgCKgCCgCioAioAgoAoqAIqAIKAKKQDgQYLzExbGqNMaSgU9KioAioAgoAopAKQJFWGOshJIioAgoAoqAIuCKANNa/AFmgOArYFtoDMQ6g6DeAn8HvhR8FZiBeJ+B64NJLcETwcyl9CF4b7CSIqAIKAKKQB4hUIS22MLBuT4Q+xeBmTKCKT1+A18IJj0AvsJcsyLomUKG1AnMFBlKikBoENgxNDXRiigC+YnAFDRrfYwpKDi6IH0FZpqPmuDDwByJ2LSTvaJLRSAMCKigCMNT0DrkMwJ/Ohq3Dev2Ntf5/2Nizl/BB4CVFIFQIqDZY0P5WLRSEUOAIwY7I2mqVV+HC2i/ODV2IbPnto+t60IRCAUCKihC8Ri0EhFHYDXq/zGYdop70mjLWbiGGX+Z1ZQZTfuClRQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRUARUAQUAUVAEVAEFAFFQBFQBBQBRSDHCPw/0233vjR4+xMAAAAASUVORK5CYII=" } }, "cell_type": "markdown", "metadata": {}, "source": [ "# 4 スムースな静止(25点)\n", "\n", "バネと質点(mass)の運動において,地面から摩擦力(friction)が速度(velocity)に比例して働いているとする.\n", "以下のコードに「速度に比例する時間に依存しない一定の摩擦項」を加えると,質点が減衰(damping)する様子が図の通り再現される.\n", "1. friction=0.1でこの図を作成せよ.\n", "1. さらに,質点が原点を超えることなくできるだけ早く減衰するにはfrictionはどの程度の値が最適か.\n", "小数点以下一桁程度で答えよ(厳密に導かなくていいよ).\n", "1. 摩擦力と質点の振る舞いを定性的に解説せよ.\n", "\n", "これは,ロボットアームなどの静止をダンパー制御するときの振る舞いとなる.\n", "\n", "![image-3.png](attachment:image-3.png)\n", "\n", "``` python\n", "import matplotlib.pyplot as plt\n", "\n", "def my_plot(xx, vv, tt):\n", " plt.plot(tt, xx, color = 'b', linestyle='--',label=\"x-position\")\n", " plt.plot(tt, vv, color = 'r', label=\"mass velocity\")\n", " plt.legend()\n", " plt.xlabel('time')\n", " plt.ylabel('position and velocity')\n", " plt.grid()\n", " plt.show()\n", "def euler3(x0,v0):\n", " v1 = v0 +(- k * x0 ) * dt\n", " x1 = x0 + v0 * dt\n", " return [x1, v1]\n", "\n", "friction = 0.1\n", "t, dt, k=0.0, 0.01, 0.1\n", "tt,xx,vv=[0.0],[1.0],[0.0]\n", "for i in range(0,6000):\n", " t += dt\n", " x, v = euler3(xx[-1],vv[-1])\n", " tt.append(t)\n", " xx.append(x)\n", " vv.append(v)\n", "\n", "my_plot(xx, vv, tt)\n", "```" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvaElEQVR4nO3deXhU1fnA8e+bnSWAsgmiEG0AlbDvKqsiglutQhVEtIqKoq3VgtW6Vq21rdal2taKe0H0JyLSukBSsVpZFNnBCAGiyI5JgECW9/fHGXDIOklmcmd5P89zn8xd5s57ksm8c8859xxRVYwxxsSuOK8DMMYY4y1LBMYYE+MsERhjTIyzRGCMMTHOEoExxsS4BK8DqKkWLVpohw4davXcffv20ahRo+AG5BErS/iJlnKAlSVc1aUsS5cu3amqLSvaF3GJoEOHDixZsqRWz83KymLIkCHBDcgjVpbwEy3lACtLuKpLWURkU2X7rGrIGGNinCUCY4yJcZYIjDEmxkVcG4ExxntFRUXk5uZSWFjodSjVatq0KWvWrPE6jKAIpCwpKSm0a9eOxMTEgM9ricAYU2O5ubmkpqbSoUMHRMTrcKqUn59Pamqq12EERXVlUVV27dpFbm4uaWlpAZ83ZFVDIvK8iGwXkZWV7BcReUJEskVkuYj0DFUsxpjgKiwspHnz5mGfBGKNiNC8efMaX6mFso3gBWBkFfvPBdJ9yyTgmRDGYowJMksC4ak2f5eQJQJV/QjYXcUhFwIvqfM/oJmItAlVPACLFh3LygqvT4wxJnZ52UZwPLDFbz3Xt21r2QNFZBLuqoHWrVuTlZVVqxd86aUM7rhD+eUv1zFq1He1Oke4KCgoqPXvIdxES1mipRxQfVmaNm1Kfn5+/QVUByUlJXWK9ayzzuLDDz9k06ZNfPbZZ4wZMwaAzz//nH/+8588+uijwQq1WoGWpbCwsGbvRVUN2QJ0AFZWsm8ucIbf+nygd3Xn7NWrl9bWW299rCNGqMbFqX70Ua1PExYyMzO9DiFooqUs0VIO1erLsnr16voJJAjy8vKCcp7MzEwdPXp0UM5VW4GWpaK/D7BEK/lc9fI+gm+AE/zW2/m2hUyzZkXMmgVpaXD11XDoUChfzRgTKosXL6Zr164UFhayb98+TjvtNFaWqfedOHEi119/PYMHD6Zjx47MnTsXcN+Wr7rqKjIyMujRoweZmZkArFq1ir59+9K9e3e6du3KV199BUDjxo0BmDZtGgsXLqR79+489thjZGVlcd555wGwe/duLrroIrp27Ur//v1Zvnw5APfeey9XX301Q4YM4aSTTuKJJ56ol99PTXlZNTQHuElEZgD9gO9VtVy1ULA1aQJPPgmjRsELL8CkSaF+RWOiX0XD34wZA5Mnw/797v+trIkT3bJzJ1xyydH7qqvV6NOnDxdccAF33XUXBw4cYPz48XTp0qXccTk5OWRmZrJ9+3aGDh1KdnY2Tz/9NCLCihUrWLt2LSNGjGD9+vU8++yz3HLLLYwbN45Dhw5RUlJy1Ll+97vf8Yc//OFIQvGvernnnnvo0aMHs2fPZsGCBUyYMIFly5YBsHbtWjIzM8nPz6dTp07ccMMNNerjXx9ClghE5J/AEKCFiOQC9wCJAKr6LDAPGAVkA/uBq0IVS1nnngsffFDxm9cYExnuvvtu+vTpQ0pKSqXftMeMGUNcXBzp6emcdNJJrF27lo8//pgpU6YA0LlzZ9q3b8/69esZMGAADz74ILm5uVx88cWkp6cHHMvHH3/Mm2++CcCwYcPYtWsXeXl5AIwePZrk5GSSk5Np1aoV27Zto127dnUsfXCFLBGo6mXV7FfgxlC9fnXOOsurVzYm+lT1Db5hw6r3t2hR/RVARXbt2kVBQQFFRUUUFhby0EMP8e677wIc+TZetitlVV0rL7/8cvr168e7777LqFGj+Otf/8qwYcNqHlgZycnJRx7Hx8dTXFxc53MGW0yPNfSPf8D553sdhTGmNq677joeeOABxo0bx9SpU3nwwQdZtmzZkSQAMGvWLEpLS/n666/ZsGEDnTp14swzz+TVV18FYP369WzevJlOnTqxYcMGTjrpJG6++WYuvPDCI/X8h6WmplbaY8f/nFlZWbRo0YImTZqEpuAhENNDTBQUwNy5sHIlVFC9aIwJUy+99BKJiYlcfvnllJSUMHDgQBYsWFDuG/yJJ57I0KFDKSgo4NlnnyUlJYXJkydzww03kJGRQUJCAi+88ALJycm8/vrrvPzyyyQmJnLcccfx61//+qhzde3alfj4eLp168bEiRPp0aPHkX2HG4W7du1Kw4YNefHFF+vl9xA0lXUnCtelLt1Hy3aJ275dNSFB9bbban1Kz8RSV8VIES3lUI2O7qNXXnmlzpo1K2jdR8NBNHYf9VzLlnDeefDyyxCG1XbGGFMvYjoRAIwfD9u2wcKFXkdijAmmF154gUvK9ks1FYr5RDByJEyYAM2aeR2JMcZ4I6YbiwEaNYJIa9cxxphgivkrgsPWrIHNm72Owhhj6p8lAiAvD7p1c0NPGGNMrLFEgBt/aMgQ8N2UaIwxQTFx4kTeeOONGj9vyZIl3HzzzYC7Qe2TTz4JdmhHsUTgc+65rnpo0yavIzHGxLrevXsfGT/JEkE9GumbVPPf//Y2DmNM9XJycujcuTMTJ06kY8eOjBs3jg8//JDTTz+d9PR0Fi1aBMCiRYsYPnw4PXr0YODAgaxbtw6oeMjpffv2MXr0aLp160aXLl2YOXPmUa+5du1a+vbte1QMGRkZACxdupTBgwfTq1cvzjnnHLZuLT+Q8vz58+nRowcZGRlcffXVHDx4EHBDag8cOJBu3brRt29f8vPzjwxxnZOTw7PPPstjjz1G9+7d+eSTT0hLS6OoqAiAvLy8o9ZrK+Z7DR3WuTO0b+8SwXXXeR2NMRHk5z8Hv/F9gqJ7d3j88SoPyc7OZtasWTz//PP06dOH1157jY8//pg5c+bw0EMPMXv2bDp37sx7773HMcccw4cffsivf/1r3nzzzQqHnJ43bx5t27Y9MnDd999/f9Trde7cmUOHDrFx40bS0tKYOXMmY8eOpaioiClTpvD222/TsmVLZs6cyZ133snzzz9/5LmFhYVMnDiR+fPn07FjRyZMmMAzzzzD5MmTGTt2LDNnzqRPnz7k5eXRoEGDI8/r0KED119/PY0bN+a2224jPz+fIUOG8O6773LRRRcxY8YMLr744joPa21XBD4i8NZbMH2615EYYwKRlpZGRkYGcXFxnHbaaQwfPhwRISMjg5ycHMB9mE+YMIEuXbrwi1/8glWrVgEwYMAAHnroIR555BE2bdpEgwYNyMjI4IMPPmDq1KksXLiQpk2blnvNMWPGHLlSOJwI1q1bx8qVKzn77LPp3r07v/3tb8nNzT3qeevWrSMtLY2OHTsCcOWVV/LRRx+xbt062rRpQ58+fQBo0qQJCQlVfz+/5pprmO77oJo+fTpXXVX3EfztisCP3xhSxphAVfPNPVT8h3eOi4s7sh4XF3dkqOff/OY3nHnmmbzzzjvk5OQwxDcJSWVDTn/++efMmzePu+66i+HDh3P33Xcf9Zpjx47l0ksv5eKLL0ZESE9PZ8WKFZx22ml8+umn9VLu008/nZycHLKysigpKalwQp6asiuCMh5/HP7+d6+jMMYEw/fff0/btm0BN+TEYRUNOf3tt9/SsGFDxo8fz+23387nn39e7nwnn3wy8fHxPPDAA4wdOxaATp06sWPHjiOJoKio6MiVx2GdOnUiJyeH7OxsAF5++WUGDx5Mp06d2Lp1K4sXLwYgPz+/3HwFFQ1/PWHCBC6//PKgXA2AJYJyZs+Gv/7V6yiMMcHwq1/9invvvZcePXoc9QH7+uuv06VLF7p3787KlSuZMGECK1asONKAfN9993HXXXdVeM6xY8fyyiuvMGbMGACSkpJ44403mDp1Kt26dTvSqOsvJSWF6dOnc+mllx6pzrr++utJSkpi5syZTJkyhW7dunH22WdTWFh41HPPP/983nrrraPOO27cOPbs2cNll1U5/1fgKhuWNFyXYA5DXZH77lMVUd29u9YvUy9iacjjSBEt5VCNjmGoD4vGYahnzZql48ePr/Q4G4a6joYNA1X4z3+8jsQYY8qbMmUK06ZN4ze/+U3QzmmNxWX07evmWF2wAC66yOtojDHmaE+GYCwcuyIoIynJTWy/f7/XkRgT3lxtgwk3tfm72BVBBWbPdvcVGGMqlpKSwq5du2jevDli/yxhQ1XZtWsXKSkpNXqeJYIKHH5fq1pCMKYi7dq1Izc3lx07dngdSrUKCwtr/MEYrgIpS0pKCu3atavReS0RVGLUKDfkxDPPeB2JMeEnMTGRtLQ0r8MISFZWFj2i5G7RUJXF2ggqkZwM77/vdRTGGBN6lggqMXQobNhgw1IbY6JfSBOBiIwUkXUiki0i0yrYf6KIZIrIFyKyXERGhTKemhg82P1cuNDbOIwxJtRClghEJB54GjgXOBW4TEROLXPYXcDrqtoD+Cnwl1DFU1NdukCzZvDRR15HYowxoRXKxuK+QLaqbgAQkRnAhcBqv2MUaOJ73BT4NoTx1Eh8PNx2G5xwgteRGGNMaIUyERwPbPFbzwX6lTnmXuB9EZkCNALOCmE8NXbnnV5HYIwxoSehujtQRC4BRqrqNb71K4B+qnqT3zG3+mL4o4gMAP4BdFHV0jLnmgRMAmjdunWvGTNm1CqmgoICGjduXKPn7N6dSElJHC1bHqzVa4ZKbcoSrqKlLNFSDrCyhKu6lGXo0KFLVbV3hTsrG42urgswAHjPb/0O4I4yx6wCTvBb3wC0quq8oR591F9RkWrDhqpTptT6JUMmlka6jBTRUg5VK0u4qktZ8Gj00cVAuoikiUgSrjF4TpljNgPDAUTkFCAFCJtbFRMSYMAAazA2xkS3kCUCVS0GbgLeA9bgegetEpH7ReQC32G/BK4VkS+BfwITfZkrbAwaBMuXw549XkdijDGhEdIhJlR1HjCvzLa7/R6vBk4PZQx1NWiQG3Pov/+F887zOhpjjAk+u7O4Gv36QWKiVQ8ZY6KXDTpXjQYN4I03oGtXryMxxpjQsEQQgAsuqP4YY4yJVFY1FICCAnjuOVixwutIjDEm+CwRBKC0FK67zlURGWNMtLFEEIAmTaBHD2swNsZEJ0sEARo0CP73PzgYXiNNGGNMnVkiCNCgQVBYCEuWeB2JMcYEV7WJQEQy6iOQcHfGGW4i+y+/9DoSY4wJrkC6j/5FRJKBF4BXVfX70IYUnlq0gJ074dhjvY7EGGOCq9orAlU9ExgHnAAsFZHXROTskEcWhiwJGGOiUUBtBKr6FW5ayanAYOAJEVkrIheHMrhws26du7ls2TKvIzHGmOAJpI2gq4g8hhtBdBhwvqqe4nv8WIjjCyupqfDOO7BggdeRGGNM8ARyRfAk8DnQTVVvVNXPAVT1W9xVQsxo2xZ+9CO7n8AYE10CSQRvqerLqnrg8AYRuQVAVV8OWWRhatAgWLjQ3W1sjDHRIJBEMKGCbRODHEfEGDQIdu+G1au9jsQYY4Kj0u6jInIZcDmQJiL+U0ymArtDHVi4GjTITV+Zn+91JMYYExxV3UfwCbAVaAH80W97PrA8lEGFs7Q0+OQTr6MwxpjgqTQRqOomYBMwoP7CiRyFhZCc7O42NsaYSFZpG4GIfOz7mS8ieX5Lvojk1V+I4Wf2bGjaFLKzvY7EGGPqrqorgjN8P1PrL5zI0LkzHDrkupGmp3sdjTHG1E0gN5T1F5FUv/VUEekX2rDCW6dO0LKl3U9gjIkOgXQffQYo8Fvf59sWs0Rc7yFLBMaYaBBIIhBV1cMrqlqKTXrPoEGQkwObN3sdiTHG1E0gH+gbRORmfrgKmAxsCF1IkeHcc2H/fkhJ8ToSY4ypm0CuCK4HBgLf+JZ+wKRQBhUJ0tNh2jRo1crrSIwxpm6qvSJQ1e3AT+shloizZ4+buvLsmJydwRgTLQLpNdRORN4Ske2+5U0RaRfIyUVkpIisE5FsEZlWyTFjRGS1iKwSkddqWgAvPfccjBgB27Z5HYkxxtReIFVD04E5QFvf8o5vW5VEJB54GjgXOBW4TEROLXNMOnAHcLqqngb8vCbBe23QIPdz4UJv4zDGmLoIJBG0VNXpqlrsW14AWgbwvL5AtqpuUNVDwAzgwjLHXAs8rap74Eg1VMTo2RMaNrRupMaYyBZIr6FdIjIe+Kdv/TJgVwDPOx7Y4reei2to9tcRQET+C8QD96rqv8ueSEQm4Wugbt26NVlZWQG8fHkFBQW1fm5lTjmlK/PmJXHxxUuCet7qhKIsXomWskRLOcDKEq5CVhZVrXIB2uOqhnYA24HZwIkBPO8S4Dm/9SuAp8ocMxd4C0gE0nCJo1lV5+3Vq5fWVmZmZq2fW5n771cVUd29O+inrlIoyuKVaClLtJRD1coSrupSFmCJVvK5GkivoU3ABbXIMd8AJ/itt/Nt85cLfKaqRcBGEVkPpAOLa/F6nrj6avjJT6BZM68jMcaY2qlqYponAa1sv6reXM25FwPpIpKGSwA/xU104282rqppuoi0wFUVRdTNascf7xZjjIlUVV0R1KnSW1WLReQm4D1c/f/zqrpKRO7HXaLM8e0bISKrgRLgdlUNpP0hrCxYAPPnw4MPeh2JMcbUXFXDUL/ovy4iDVV1f01OrqrzgHlltt3t91iBW31LxFq0CB56CG6+GVq39joaY4ypmUBuKBvg+8a+1rfeTUT+EvLIIsjw4e7nggXexmGMMbURyH0EjwPn4OsyqqpfAoNCGFPE6dnTNRZbIjDGRKJAEgGquqXMppIQxBKx4uNhyBDXTmCMMZEmkESwRUQGAioiiSJyG7AmxHFFnGHDIC7ODURnjDGRJNBhqG/E3Sn8DdDdt278TJ7sJrM/5hivIzHGmJoJZIgJUdVxIY8kwsXHex2BMcbUTiBXBP8VkfdF5Gci0izUAUWyp56CjAzQSm/DM8aY8FNtIlDVjsBdwGnA5yIy1zcInSmjUSNYuRJWr/Y6EmOMCVygvYYWqeqtuKGldwMvVvOUmDRsmPtpvYeMMZEkkBvKmojIlSLyL+ATYCsuIZgy2reHk0+GDz/0OhJjjAlcII3FX+IGh7tfVT8NbTiR7+yz4ZVX4NAhSEryOhpjjKleIIngJN+YQCYAY8ZAgwawb58lAmNMZAhkPgJLAjUwdKhbjDEmUgTUWGxqpqgIltTvzJXGGFNrlghC4PHHoU8f+PZbryMxxpjqhXKGspg1YgT86lfw/vswcaLX0RhjTNWquiJYAiwFUoCewFe+pTtgzaBV6NoVjjsO/v1vryMxxpjqVTtDmYjcAJyhqsW+9WeBhfUTXmQSgZEj4e23oaTExiEyxoS3QNoIjgGa+K039m0zVRg50g1JvXix15EYY0zVArmP4HfAFyKSCQhudrJ7QxlUNBg5EjIz3exlxhgTzgK5j2C6b3iJfr5NU1X1u9CGFfmaNnWzlhljTLgLtPtoPLAD2AN0FBGbszgAOTlw++2wdavXkRhjTOWqvSIQkUeAscAqoNS3WYGPQhhXVMjPhz/8ATp2hGuv9ToaY4ypWCBtBBcBnVT1YIhjiTpdukCHDjBnjiUCY0z4CqRqaAOQGOpAopEIXHihG5Z63z6vozHGmIoFkgj2A8tE5K8i8sThJdSBRYsLLoDCQpujwBgTvgJJBHOAB3CT0iz1W6olIiNFZJ2IZIvItCqO+4mIqIj0DuS8keTMM1310LZtXkdijDEVC6T7aK2mpRSReOBp4GwgF1gsInNUdXWZ41KBW4DPavM64S4xETZscNVExhgTjgKZqjJdRN4QkdUisuHwEsC5+wLZqrpBVQ8BM4ALKzjuAeARoLBGkUeQw0mgqMjbOIwxpiKB9BqaDtwDPAYMBa4isCql44Etfuu5/HBTGgAi0hM4QVXfFZHbKzuRiEwCJgG0bt2arKysAF6+vIKCglo/ty5KSoRrr+1Fv367ue66QHJo9bwqSyhES1mipRxgZQlXISuLqla5AEt9P1eU3VbN8y4BnvNbvwJ4ym89DsgCOvjWs4De1Z23V69eWluZmZm1fm5djRihevLJqqWlwTmfl2UJtmgpS7SUQ9XKEq7qUhZgiVbyuRrIN/uDIhIHfCUiN4nIj3EDz1XnG+AEv/V2vm2HpQJdgCwRyQH6A3OiscEY4NJL4euvYdkyryMxxpijBZIIbgEaAjcDvYDxwJUBPG8xkC4iaSKSBPwU1wMJAFX9XlVbqGoHVe0A/A+4QFWjcpLHH//YDUf9+uteR2KMMUerNhGo6mJVLVDVXFW9SlV/oqr/C+B5xcBNwHvAGuB1VV0lIveLyAV1Dz2yNG8Ow4fDrFmglc77Zowx9S+QxuJaU9V5wLwy2+6u5NghoYwlHNx+O+zcCaWlNlmNMSZ8hDQRmKOddZbXERhjTHmBDkNtgmTLFvjLX6x6yBgTPgIZhrolcC3Qwf94Vb06dGFFr/fegxtvhP79bfYyY0x4COSK4G2gKfAh8K7fYmrhxz92w0689prXkRhjjBNIG0FDVZ0a8khiRPPmMHo0vPoq/O53kGCtNMYYjwVyRTBXREaFPJIYMmECfPcdzJ/vdSTGGBP4DWVzRaRQRPJ9S16oA4tmo0ZBq1awfLnXkRhjTGDDUKfWRyCxJDnZTWzfoIHXkRhjTID3EfjuBB7kW81S1bmhCyk2HE4CRUWu8dgYY7wSyHwEv8NVD632LbeIyMOhDiwWXHMNnHuu11EYY2JdIG0Eo4CzVfV5VX0eGAmMDm1YsaFDB9dg/PXXXkdijIllgd5Z3MzvcdMQxBGTrrrKjTn09797HYkxJpYFkggeBr4QkRdE5EXcxPUPhjas2HD88XDeeTB9Ohw65HU0xphYFcgw1P/ETRrzf8CbwABVnRnqwGLFddfB9u0we7bXkRhjYlWliUBEOvt+9gTa4OYczgXa+raZIBgxAp56CoYN8zoSY0ysqqr76K24CeP/WME+BeyjKwji490gdMYY45VKE4GqTvI9PFdVC/33iUhKSKOKQa+9Bnv3wuTJXkdijIk1gTQWfxLgNlMHc+fCtGmQZ4N3GGPqWVVtBMeJSC+ggYj0EJGevmUIbjJ7E0S33gr5+fD8815HYoyJNVW1EZwDTATaAX/y254P/DqEMcWk3r3hzDPhz3+Gm26y4amNMfWnqjaCF4EXReQnqvpmPcYUs2691U1cM3s2XHKJ19EYY2JFpYlARMar6itABxG5tex+Vf1TBU8zdXD++e4GsxRrijfG1KOqKiAa+X42ro9AjOtK+s47XkdhjIk1VVUN/dX38776C8cAFBTAnDlw+eVeR2KMiQWBDEP9exFpIiKJIjJfRHaIyPj6CC5WPf88jBsHH3/sdSTGmFgQyH0EI1Q1DzgPyAF+BNweyqBi3TXXQMuW8MADXkdijIkFgSSCw9VHo4FZqvp9oCcXkZEisk5EskVkWgX7bxWR1SKy3He10T7Qc0ezhg3httvg/ffhs8+8jsYYE+0CSQRzRWQt0AuYLyItgcJqnoOIxANPA+cCpwKXicipZQ77Auitql2BN4Df1yT4aDZ5MrRoAXfcAapeR2OMiWaBDEM9DRiI+8AuAvYBFwZw7r5AtqpuUNVDwIyyz1PVTFXd71v9H+7mNQM0bgx33+2SQH6+19EYY6KZaDVfN0UkEbiBHyav/w/wrC8pVPW8S4CRqnqNb/0KoJ+q3lTJ8U8B36nqbyvYNwk3EiqtW7fuNWPGjCpjrkxBQQGNG0dOb9jSUhBxS1mRVpaqREtZoqUcYGUJV3Upy9ChQ5eqau+K9gUykMEzQCLwF9/6Fb5t19Qqmgr4eiH1BgZXtF9V/wb8DaB37946ZMiQWr1OVlYWtX2ul3JzYeNGNwTFYZFalopES1mipRxgZQlXoSpLIImgj6p281tfICJfBvC8b4AT/Nbb+bYdRUTOAu4EBqvqwQDOG3OuuALWr4e1ayE11etojDHRJpDG4hIROfnwioicBJQE8LzFQLqIpIlIEvBTYI7/ASLSA/grcIGqbg887Njy8MPw7bdwn93aZ4wJgUASwe1Apohkich/gAXAL6t7kqoWAzcB7wFrgNdVdZWI3C8iF/gOexQ3hMUsEVkmInMqOV1M698frr0WHn8cVqzwOhpjTLSptmpIVeeLSDrQybdpXaBVOKo6D5hXZtvdfo/PqkGsMe3hh+H//s91K/3oI6+jMcZEk2oTgW9aysnAGbi5iheKyLNlp680odW8OTz6KCxcCAcOeB2NMSaaBNJY/BJuMponfeuXAy8Dl4YqKFOxq65yizHGBFMgiaCLqvrfEZwpIqtDFZCp3pdfwmOPpTNoEMQF0spjjDFVCORj5HMR6X94RUT6AUtCF5KpzpIlMGfO8Tz5ZPXHGmNMdQJJBL2AT0QkR0RygE+BPiKyQkSWhzQ6U6Grr4b+/Xfxq1/B0qVeR2OMiXSBVA2NDHkUpkZEYNq0tUyZcjqXXuqSwTHHeB2VMSZSBTLo3KaqlvoI0pTXtGkRr7/uhp/4vY3Zaoypg0CuCEyY6t8fPvgABgzwOhJjTCSzPicRbvBgSEqCXbtg3rzqjzfGmLIsEUSJqVPhoosgM9PrSIwxkcYSQZT4wx8gPR0uvNCmtzTG1IwlgijRrJmb47hVKzjnHOtWaowJnCWCKHL88bBggetK+otf2FzHxpjAWK+hKHPiiS4ZJCRUPMWlMcaUZVcEUSgtDU44AUpKYNw4ePFFryMyxoQzSwRRrLAQtm2DiRPhllugqMjriIwx4cgSQRRr1Aj+9S/4+c/hiSdg+HD4ptys0caYWGeJIMolJsJjj8Grr7pRS885B0pLvY7KGBNOrLE4Rlx+OfTu7aqK4uJctdGePdCmjdeRGWO8ZlcEMaRjRzjzTPf4T3+CTp3cjWiHDnkblzHGW5YIYtSYMW6cottvhy5d4JVXXC8jY0zssUQQo370I3jnHZg7Fxo0gCuugGuv9ToqY4wXLBHEuNGj4Ysv4I034IYb3LaNG+Huu62HkTGxwhKBIS4OfvIT6NPHrS9YAL/9rbtL+ayzYPp0+P57b2M0xoSOJQJTzs9+BtnZcNddkJPj5khu3/6HG9L27rVxjIyJJtZ91FTopJPgvvvg3nth8WJYtszdkwAwYgRs3w6nnw4DB7olI8ONb2SMiTwhvSIQkZEisk5EskVkWgX7k0Vkpm//ZyLSIZTxmJoTgb59YdIkt67qrhh693aT4Nx0E/Ts6YaxOOypp9yQ2Bs3QnGxJ2EbY2ogZN/hRCQeeBo4G8gFFovIHFVd7XfYz4A9qvojEfkp8AgwNlQxmboTgeuuc4sqbN4Mn3wCLVq4/Tt2wJQpPxwfHw/t2rlqpmuugbw8eOklaN36h6WgIIHSUtdWYYypf6G8mO8LZKvqBgARmQFcCPgngguBe32P3wCeEhFRDUENdGkpDTdtCvppY5mIazto3/6HbS1bumSwciV8/bW7Kti40X3gg3vsnyicM3juOXelsWwZjB8PTZu6pUkTN2bSTTdBjx7w1VcwcyakpBy9nH22e41vv1G+WFpKQrySEFdKfJySGF9K1wwltVEpu3Yq32xx++MoJU4U0VJObFdKcpKyd6+b//nwEN4ibmnXzlV97d2j5OUdvR+gbVtI2rWLvHVb2Vfww9v38P7Wrd3j77+HAwd8O31vc5Effj9797q7vo9QJS7OTTgE7m7wgweP/u3FxyktW7rHu3eXv0EwMRGaN3ePd+4sf5Xmv3/HDnc/ScGaPXzXaAsAyclujgtwd6aX/e9MSXETIwF89x3lVLe/YUP3dy4tdecvq1Ejt7+kxFVJlpWaCo0buzasnTvL79fNO2HLFg4dcn/bspo1c12oDx50v7+yjjnGlaGwsOL9xx7r9h844P4+ZbVo4eYV37ev4k4XLVu6v0FBgfuiVFarVu69V1QE8UfePMElofjMBRCRS4CRqnqNb/0KoJ+q3uR3zErfMbm+9a99x1Tw53R69+6tS5YsqXlA995L6W9/S9zOnT+8KyNYVlYWQ4YM8ToM96lw4IB7h5dd8vJg/36337fovv0c2HOAg3sOcCjvACV5+8nfuZd2xyXQKKmYA3lFfPdNMRQXE1dchJQWE19aRPOmxaTEFXHoQDGH9hWRQDGC78P8yAe6tWCb6Lb+F7+g45/+VKvnishSVe1d0b6IaN4TkUnAJIDWrVuTlZVV43M0admSniUlrP7jH9k+fHiQI6x/BQUFtfo9BCp+/35SvvuO5G3bSN69m8Q9e0jyLYcfJ+7dS0JBAXEB3pJcmphISXIy8cnJpCQnk5SURElyMsnN4jmkyRwsjkdTE0g9LRlNSEDj4yn1/dwdH4/Gx6MJCZTEJVBCPMWl8ZSUxlFcGkeDRqUkJELhoQTyChJdelBfqlChRasikpKVvIIkduxMppR4ShVU4lEROqTtJylZ2bU7me3bk4/EfDi3pKcXkJiobN+ezPadKb6dP+w/5ZQ8iosL2b27GTt3JqEi4JeXMrp+jwC53zRg966kH86PEBenZGS4r4KbNzdkz57Eo35vCYlw2mlu/8aNDcnLO3p/UnIpnTvnA/D1140pKDj637pBgxI6dnT7v/oqlQMH4o/a36hRMSefXADAunWpHDwYT1FxEYkJ7nVSU4tIS9sHwOrVTSguProOr2nTQ7Rvvx+AlSubUlp69IxIxx57iHbt3P7ly5tRVosWB2nb9gClpcLKVU3L7W/VqpDjWhdSVBzHmjVNyu1vc9wBWrY8yMGDcaxbX35/yxZ7adMGDhxI4KvsxuX2n3jCfpo1O8S+fQl8vaH8/g7t99GkSRF5eYnkbGpYbv9Jafto3LiYvXsT2byl/P70H+2jQYNidu9OJveblHL7O3UsIDm5hB07k9m6tfz+Uzrnk5hYysGD8exJb8+3Ifi/D+UVwQDgXlU9x7d+B4CqPux3zHu+Yz4VkQTgO6BlVVVDtb4iKC3lUIsWJI0YATNm1Pz5YSYoVwS7d8OaNbB6NaxdCxs2wKZNrs9oRde4TZr8ULHfqpVbmjVzdTiHf/ovTZq46/4GDdy1c3x8+XMGqyxhIFrKAVaWcFWXsnh1RbAYSBeRNOAb4KfA5WWOmQNcCXwKXAIsCEn7AEBcHDsHDqTtvHmuMjA5ufrnRAtVV2G/ZInrC/rFF+7D379CtkED12e0fXvo3x86dPihAaBNG/eh36CBZ0UwxoROyBKBqhaLyE3Ae0A88LyqrhKR+4ElqjoH+AfwsohkA7txySJkdp5xBm3ffdfdOnvuuaF8KW8dPAiLFkFWFnz0kUsAe/e6fSkp0K2bG1vi1FPhlFPc0r69ddsxJkaFtI1AVecB88psu9vvcSFwaShj8Le3Z09XXTFjRnQlAlXXTWfuXJg/3/XnPHDAdUfp2hXGjnUd//v0cR/+iYnVn9MYEzMiorE4WEqTktyH4muvwdNPuz5nEUqKi+HDD2HOHDeMaE6O29G1q7v7a8gQGDTI9W0zxpgqxFQiAODKK+Hvf4c333SPI4kqfPopvPYaA155xXXRTElxI8PdcQecd57r0G6MMTUQe4lg4EA3GP/06ZGTCLZsgeeeg5dfdndkpaSwt39/Wt1yixv4p2H5LmvGGBOo2GsdFHG3sP7nP7B8udfRVK60FP71L7jgAteD54EHXAJ74QXYto3V99wDF11kScAYU2exlwjA1aE3bAiPPeZ1JOXt3w9PPOE+9EeNgs8+g6lTXR//9993VzFNyt80Y4wxtRWbieDYY90g+6++Ct9+63U0Tn4+/P73kJYGt9wCxx/vejdt2QIPPeSuCowxJgRiMxEA3Hqra3x94AFv49izB+6/3/XjnzoVund31VYLF7oeTklJ1Z7CGGPqInYTQVoaXH+960G0bl39v/6OHXDnne6b/j33wBlnuGqg995z3T6NMaaexG4iAPjNb9ywCT//ef3Nvbh1K/zyly4BPPwwnHOOG/Jhzhw3A4wxxtSz2E4ErVq5+vd//xtefDG0r7V5sxtUPy0N/vxnN1v8qlXw+uuuOsgYYzwS24kA4MYb4cwz4eab3UicwZad7abmOvlk+Nvf4IorXFXUSy+5MX6MMcZjlgji4lzvoQYN4PzzK54CqTa+/BIuuww6dYJXXnHtEdnZrk3i5JOD8xrGGBMElggATjgB3nrLdSUdPBhyc2t3HlU32ufo0a66Z+5c1ztp40Z48kk48cSghm2MMcFgieCwgQNdW0Furpsc9+23A29A3rPHfdBnZLhEsmiR65a6eTM8+qgbz98YY8KUJQJ/gwa5D/G2bd3wDUOGuGqjsjNeFxXB55/DU0+5sX5atXJtDA0auKqfTZvgrrt+mPHbGGPCWOwNOledU05xyeC55+CRR2D8eLe9VStITXXj/G/bBofn6U1Pd9U/Y8dCz57exW2MMbVkiaAiycmuN9ENN7hhnz/9FNavh3373LDPbdpAly7Qr5+7H0Ck2lMaY0y4skRQlbg4OP10txhjTJSyNgJjjIlxlgiMMSbGWSIwxpgYZ4nAGGNinCUCY4yJcZYIjDEmxlkiMMaYGGeJwBhjYpxofc3MFSQisgPYVMuntwB2BjEcL1lZwk+0lAOsLOGqLmVpr6otK9oRcYmgLkRkiar29jqOYLCyhJ9oKQdYWcJVqMpiVUPGGBPjLBEYY0yMi7VE8DevAwgiK0v4iZZygJUlXIWkLDHVRmCMMaa8WLsiMMYYU4YlAmOMiXExkwhEZKSIrBORbBGZ5nU8NSEiz4vIdhFZ6bftWBH5QES+8v0M+wmSReQEEckUkdUiskpEbvFtj8SypIjIIhH50leW+3zb00TkM9/7bKaIJHkdayBEJF5EvhCRub71SC1HjoisEJFlIrLEty3i3l8AItJMRN4QkbUiskZEBoSqLDGRCEQkHngaOBc4FbhMRE71NqoaeQEYWWbbNGC+qqYD833r4a4Y+KWqngr0B270/R0isSwHgWGq2g3oDowUkf7AI8BjqvojYA/wM+9CrJFbgDV+65FaDoChqtrdr799JL6/AP4M/FtVOwPdcH+f0JRFVaN+AQYA7/mt3wHc4XVcNSxDB2Cl3/o6oI3vcRtgndcx1qJMbwNnR3pZgIbA50A/3F2fCb7tR73vwnUB2vk+VIYBcwGJxHL4Ys0BWpTZFnHvL6ApsBFfh55QlyUmrgiA44Etfuu5vm2RrLWqbvU9/g5o7WUwNSUiHYAewGdEaFl81SnLgO3AB8DXwF5VLfYdEinvs8eBXwGlvvXmRGY5ABR4X0SWisgk37ZIfH+lATuA6b4qu+dEpBEhKkusJIKopu7rQcT0AxaRxsCbwM9VNc9/XySVRVVLVLU77ht1X6CztxHVnIicB2xX1aVexxIkZ6hqT1w18I0iMsh/ZwS9vxKAnsAzqtoD2EeZaqBgliVWEsE3wAl+6+182yLZNhFpA+D7ud3jeAIiIom4JPCqqv6fb3NEluUwVd0LZOKqUJqJSIJvVyS8z04HLhCRHGAGrnroz0ReOQBQ1W98P7cDb+ESdCS+v3KBXFX9zLf+Bi4xhKQssZIIFgPpvp4QScBPgTkex1RXc4ArfY+vxNW3hzUREeAfwBpV/ZPfrkgsS0sRaeZ73ADX1rEGlxAu8R0W9mVR1TtUtZ2qdsD9XyxQ1XFEWDkARKSRiKQefgyMAFYSge8vVf0O2CIinXybhgOrCVVZvG4UqcfGl1HAelw97p1ex1PD2P8JbAWKcN8Ufoarx50PfAV8CBzrdZwBlOMM3KXscmCZbxkVoWXpCnzhK8tK4G7f9pOARUA2MAtI9jrWGpRpCDA3Usvhi/lL37Lq8P95JL6/fHF3B5b43mOzgWNCVRYbYsIYY2JcrFQNGWOMqYQlAmOMiXGWCIwxJsZZIjDGmBhnicAYY2KcJQJjquAbAXKy73FbEXnD65iMCTbrPmpMFXxjIs1V1S5ex2JMqCRUf4gxMe13wMm+weW+Ak5R1S4iMhG4CGgEpAN/AJKAK3BDVI9S1d0icjJuCPSWwH7gWlVdW9+FMKYqVjVkTNWmAV+rG1zu9jL7ugAXA32AB4H96gYI+xSY4Dvmb8AUVe0F3Ab8pT6CNqYm7IrAmNrLVNV8IF9Evgfe8W1fAXT1jbI6EJjlhlkCILn+wzSmapYIjKm9g36PS/3WS3H/W3G4cf2713NcxtSIVQ0ZU7V8ILU2T1Q318JGEbkU3OirItItmMEZEwyWCIypgqruAv4rIiuBR2txinHAz0Tk8IiYFwYzPmOCwbqPGmNMjLMrAmOMiXGWCIwxJsZZIjDGmBhnicAYY2KcJQJjjIlxlgiMMSbGWSIwxpgY9/90hwBOZIcXAwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "def my_plot(xx, vv, tt):\n", " plt.plot(tt, xx, color = 'b', linestyle='--',label=\"x-position\")\n", " plt.plot(tt, vv, color = 'r', label=\"mass velocity\")\n", " plt.legend()\n", " plt.xlabel('time')\n", " plt.ylabel('position and velocity')\n", " plt.grid()\n", " plt.show()\n", "def euler3(x0,v0):\n", " v1 = v0 +(- k * x0 - friction*v0) * dt\n", " x1 = x0 + v0 * dt\n", " return [x1, v1]\n", "\n", "friction = 0.6\n", "t, dt, k=0.0, 0.01, 0.1\n", "tt,xx,vv=[0.0],[1.0],[0.0]\n", "for i in range(0,6000):\n", " t += dt\n", " x, v = euler3(xx[-1],vv[-1])\n", " tt.append(t)\n", " xx.append(x)\n", " vv.append(v)\n", "\n", "my_plot(xx, vv, tt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": true }, "vscode": { "interpreter": { "hash": "f3f87633aac09da3bda522f97956bee375b5501d1579e6458804e567301cb62a" } } }, "nbformat": 4, "nbformat_minor": 4 }