{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br />\n",
    "\n",
    "<div style=\"text-align: center;\">\n",
    "<font size=\"7\">数値計算試験問題</font>\n",
    "</div>\n",
    "<br />\n",
    "<div style=\"text-align: right;\">\n",
    "<font size=\"4\">2022/12/23 実施</font>\n",
    "<br />\n",
    "<font size=\"4\">cc by Shigeto R. Nishitani 2022 </font>\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAAAXNSR0IArs4c6QAAAERlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAABcqADAAQAAAABAAAA/AAAAACR+ioEAAAvqUlEQVR4Ae2dB3xUVfbHfwmEZuiETgpIQDok9A6KdEFBRSyshT+6oIJYEF0FQdZCseC6rN3FAoKIVBFC7116DR2kKQlNQuZ/zpvMOsSUSTKTee/N73w+hzcz7777zvme4czNfbcAFBIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARLIEoEgD0vHS7kE0euiSaKxohQSIAESIAELEYgXW0tZyF6aSgIkQAIBQyA4YDyloyRAAiRgUwKedq0cFP/PizpE/y06STS19JcPVFGgQIGY8PDw1Odt8T45ORnBwfb9/aN/1v6aMn7Wjd+ePXvOiPVh2fHA00ReQSo/JlpadIHoINGlomlKdHS0Y/fu3Wmes/qHixcvRps2bazuRrr207900VjiBONniTClaWRQUNAGOZGt54+eNi01iav8Kvq9aCN9QyEBEiABEvA/AU8S+U1iZuEUU/V1B9FtKe95IAESIAES8DOBvB7cv4yU0Va4ipb/SnSevqGQAAmQAAn4n4AnifyAmFnX/6bSAhIgAW8QuHbtGo4ePYorV654o7pcraNo0aLYuXNnrt7T2zeTwSCoWLEiQkJCvFa1J4ncazdjRSRAAv4noEm8cOHCiIyMhDxg879BWbAgISHBsD0Ll5iqqMPhwNmzZ40f0qioKK/Z5kkfudduxopIgAT8T0Bb4iVLlrRcEvc/uZxboD+cyt7bfw0xkec8NqyBBCxHwGotccsBzsBgX7BnIs8AOE+RAAmQgBUIMJFbIUq0kQQClEDnzp3x22+/Zej9P/7xD/z8888ZlknvpE6g6tq1a3qnLfM5H3ZaJlQ0lAQCh4A+FFSdM2dOpk6PHDky0zJ2L8AWud0jTP9IwKQExo0bh1q1ahk6YcIExMfHo1q1anjwwQeNz44cOWKMrDlzRpcgAV577TU0aNAALVq0QJ8+ffD2228bn/fr1w/fffed8VpH4rzyyitGudq1a2PXrl3G52vXrkXTpk1Rv359NGvWDHZbQoQtciPM/IcEApPA008Dmzd71/d69QDJyxnKhg0b8Omnn2LNmjVGy7tx48Zo3bo19u7di88//xxNmjS54fp169Zh2rRpWLlypS7KZyTqmJiYG8q43pQqVQobN27EBx98YCT7jz76CNWrV8eyZcuQN29eoxvmxRdfNOpzXWP1IxO51SNI+0nAggSWL1+Onj174qabdNUP4M477zQSbURExF+SuJ5fsWIF7rjjDiOJ6xj4bt266cdpitalool++vTpxuvff/8dDz30kPFDoaNGdFKUnYSJ3E7RpC8kkEUCmbWcs1hdjou7EntOKsqfP79xeZ48eZCUpBuaAS+//DLatm2L77//3ujCsdsKpuwjN8LMf0iABHKTQMuWLTFjxgxcunQJFy9eNBKsfpaeNG/eHD/++KMxkSYxMRGzZs1Kr2ian2uLvEIFXY0b+Oyzz4yjnf5hIrdTNOkLCViEgD601IeUjRo1gvaPP/rooyhevHi61jds2BDdu3c3Hlh26tQJ+iBT113xVJ577jkMGzbMeNjpaqV7em3AltONJewqcXFxdnXN8Iv+WTu8nsRvx44dlnRS1llxXLhwwSEteIf0fzvkgakl/VCj04qB/GCsz+6PBvvIs0uO15EACeQqgf79+2Pbtm34448/jAeX2qqnOAkwkfObQAIkYAkCX331Fay++qGvQLOP3FdkWS8JkAAJ5BIBJvJcAs3bkAAJkICvCDCR+4os6yUBEiCBXCLARJ5LoHkbEiABEvAVASZyX5FlvSRAAukS0IWrMhNdSEsnDPladILQwIEDM7yNLner67xkVXQRL9eiX1m9NivlmcizQotlSYAEvELAk6SYnUR+/fp1r9iXupLsJvLU9fjqPRO5r8iyXhIggXQJhIaGGuc0Qeq6J7169TJWKOzbt6+xGuK7776L48ePG+uj6BopKj/99BPat29vrHzYu3dv6FR9FW31Pv/888bnU6dONep76qmnUE+WYdRlcnUJW5Vz586hR48eqFOnjrEw19atW43P3f/RZQB0pqkud3vrrbfi1KlTxtosH374IcaPH2/Uqasonj59GnfddRd0xqmqLuqlohsrd+jQATVr1jRmq8rcH/fqffaa48h9hpYVk4AFCPhrHVs3NJs2bcL27dtRvnx56JoqmhSffPJJ6HrlMlMVuiytdk+MGjUKM2fORNmyZfHGG28Y53V3IBXd0FiXrlXRpKtdMptlfd6lS5fi4YcfNiYS6TrlmqB1jZdFixYZ655rGXfRtc5Xr15tbEyty9+++eabGDt2LAYMGAD98Rk6dKhR/L777sPgwYONtdEPHz6M22+/HTt37sSIESOMz9Su2bNn4+OPP3av3mevmch9hpYVkwAJeEJA11upWLGiUVRb0fGywYQmVHfR5CrT2o3WbnBwsDG7UzeKcMk999zjemkcdeMJlVatWkGm9RvbxenSubqmuUq7du2M1rOec5ejR49C6zpx4oRxj6ioKPfT/3utW8upPS7RevQvBP3hcC2d26VLlwzXj3Fd640jE7k3KLIOErAqAROsY+tadlYRui89645Uuyhuu+02TJo0CboeeWpJvfxt6p3qU79Pfb3r/aBBgzBkyBBjgS7t9nn11Vddp244JicnGy133eTCDMI+cjNEgTaQAAn8hYAmbJ2Sr6I7BmmXy/79+433uvTtnj17jNdp/fPtt98aH2srXFdJVNVlcidPnmx8rklau2yKFClyw+Xuy93qTkUucbdFP9N+8Pfee8912ujG0Tf6F4AuJaAyd+5cnD9/3njt63+YyH1NmPWTAAlki4AuktWxY0fjgWdYWJixjrj2d+vDSu1Wce3HmVbl2lLW/nDt23b1U2vrWreY0+tfeOEFY0u51NdqGX2QqrsLaaJ3ie5IpJtSaNePPuzUh7Hr16836qpRo4bRL69ltR9eu1f0Yad2sYSHh7uqsN6Ry9jqQpXWFE+WQbWmZ06r6V/aS6haJaa6jG1mInt/OmSPz8yK+fW8t5exZYvcer+TtJgESIAEbiDAh5034OAbEiABqxPQ/u9AE7bIAy3i9JcEhID0K5CDnwj4gj0TuZ+CyduSgL8I6INAnYHoi4TiL5+scl9lruy9PWyRXStW+QbQThLwEgGdfKMTX3SaudXkypUrXk+Cuc1Ak7hrApS37s1E7i2SrIcELEIgJCQE6c1YNLsL2v+twwopNxJg18qNPPiOBEiABCxHICuJPI94t0l0VmZeFtYZV7IimUyjyqwoz5MACZAACeSQQFYS+VNyr50e3+/QIUBmZjGZe0yMBUmABEggWwQ87SPXpcm6iI4WHZLZnZKQF6dRDEGXHLjywDBcfPlLJFeMQP4alVG+ZxMUaNkQKFQos2p4ngRIgARIwAMCQR6U0SLfiY4R1WXHdEHerqKpRZrfUEU0Cse8i6a4Jgk9BEkoGXQO4Y5DKC3pXUU/P1S8Jn6NbYKgu5vgapUIyALAxjmz/6NLVboWxTe7rdmxj/5lh5p5rmH8zBOLrFoiG2hskGtis3qdp+U1aX+QUriNHDPtI4/R+QYujYhwJCU5HPHxDse8yWccn/We5fii0jDHiqDmjusIMsodK1rdsfeBEY7rR4/LMEtzC9fqMHd8MrOO8cuMkLnP2zl+klvXp+TZLB886SNvLrV2F40X/Ua0neh/RTMX7T4ZPVrWGAYipNF9+30l8dCULnjg8Ouof3E55n9yHJ81+gD7E8vi5i9fwfWK4dhRrw/Oz16Zed0sQQIkQAIkYBDwJJEPk5LaRx4peq/oItH7RTMWzdyyCDxkD760pGBBoNPfyqLfmsfR6GIcZo/fgxkVBqLCljko3rU5dldoh5PfLU/rUn5GAiRAAiTgRsCTRO5W3LOXCdHR0n6PTzeJp64lf355kvp0VfQ+Oh7H1x7DlKbjUfT4DpTt3RLbK3bA8ZnZ/osj9a34ngRIgARsRyCriXyxEEjrQafXwNzSMBR3r3waSbsPYHqztxF2bDPK3tEI6+s+gsT9p7x2H1ZEAiRAAnYhkNVEnmt+V4wuhDtXPIM/tu/D3BpDUWerDGGsGo01fcYj+Y+kXLODNyIBEiABsxMwbSJ3gatYowi6bH8TO7/9BduLNkPjb4Zgd6lmODRU9suLjARkR23OInXR4pEESCAQCZg+kbuCUvfuamh8Zg4W9f8GYQkHUW7sM1h2qBKuOcQFziJ1YeKRBEggAAlYJpFrbILzBKHdv+/B9bIVsBaN0BLLsUemHx2CbHB66RIwfHgAhpAukwAJBDoBSyVyV7DKnNqKFliBlWiC8jiBEjiHJWgFh7bMKSRAAiQQYAQsmcgRLi1wkWZYjSvIj/2ojNZYiuV5WuP3Y4kBFkK6SwIkEOgErJnIZbaoa9GtcjiFOtiKZcGt0ez6MpyOaoSd0z1fpDHQvwD0nwRIwPoErJnIdbaozhrV2aOy2FawHFt+8Rh2TFiAYtfPoNJdDbFowBTrR4cekAAJkIAHBDxdxtaDqnK5iCbzVNP/a4sJZ9pswqHWdxsPReev3Yh2K0cjpIDuiUEhARIgAXsSsGaLPINYlKpbAdVPxGF1vQG4fdMbWFehB84cuJDBFTxFAiRAAtYmYLtEruHIUzAfmmz6F9b2+wCNzs3F+epNsGv2fmtHitaTAAmQQDoEbJnIXb42+vRx7P/XApRKOoVS3ZpgzXguj+tiwyMJkIB9CNg6kWuYqg1oi6tLVuNSSDHUHdIOP/fnQ1D7fH3pCQmQgBKwfSJXJ8u2rIoSu1dhf7FY3PqfezC3zRtIvi57GFFIgARIwAYEAiKRa5xCI0uh2pGfsa7Kvei05AX8XPMpXLuabIMQ0gUSIIFAJxAwiVwDnTe0AGJ3T8bqZoPRYfd7WBnZB4lnrwb6d4D+kwAJWJxAQCVyjVVQnmA0WTEOa3q9hdYnp2BXVEecPfC7xcNI80mABAKZQMAlclewG08dio2Dv0TdhOU4WbMdTmw97TrFIwmQAAlYikDAJnKNUoNx92PXmBmofGUHLsa2QvyyI5YKHo0lARIgASUQ0IlcAdR+oQsOTfoJpZOOI2+bFtj94x79mEICJEACliEQ8IlcI1X9sZY4M3UxCuAyit3RCtu/+cUyAaShJEACJMBEnvIdqHxXfVyetxSO4Dwod18bbPl4Pb8dJEACJGAJAkzkbmGqdFt1JC9ehot5iiDy0fbY+N4Kt7N8SQIkQALmJMBEniou5VtURv41y3AuX1lUe7ID1r0Zl6oE35IACZCAuQgwkacRj9INKqLIpqU4USAKtZ7vjHWjf0qjFD8iARIgAXMQYCJPJw4la5RByS1xOFIwGrVf6o41r8xJpyQ/JgESIAH/EmAiz4B/8egwlP5lEQ4Wqon6I3tgTdEOaN2uHRAZCUyenMGVPEUCJEACuUeAiTwT1sWqlET51wdiL6qi/oXFWOtoCBw6BPTvz2SeCTueJgESyB0CTOQecC46fgQq4gj2IBoNsBGr0Qi4dAkYPtyDq1mEBEiABHxLgIncE76HD6MoEhCOQ9iF6ohxJXP5nEICJEAC/ibARO5JBMLDjVJFkIgIxGMnbjGS+drC7T25mmVIgARIwKcEmMg9wTt6NFCokFFSk3kkDhrJ3Ogzf36aJzWwDAmQAAn4jAATuSdo+/YFJk0CIiLgCApCkYiSiHhnCHaGNkT9N+/FuheYzD3ByDIkQAK+IcBE7ilXTebx8ViyaJFxLPpkP0TsmGck83pv3Iv1L073tCaWIwESIAGvEvAkkReQO64V3SK6XXSEKEUIFK1UBOHbJZnfFIu6Y+7BhpdnkAsJkAAJ5DoBTxK5bmrZTrSuaD3RjqJNRClCoFh4EVSSZL6rUAxqj7obG1+dSS4kQAIkkKsEPEnkDrEoMcWqEDmq6meUFALFI4qi4vb52F2oPmqN6IVNr80iGxIgARLINQJBHt4pj5TbIHqz6ETR50VTi0x1hCrCwsJipkyZkvq8Ld4nJiYiNDQ0TV8uHruMSo+8hGpXt2H2I2+jxP210yxn5g8z8s/MdntqG/3zlJQ5y9k5fm3bttUcG5sb5IvJTXRd11oZ3Sw6OtphV4mLi8vQtTN7zzm2F2jguIJ8js3/nJthWTOezMw/M9qcFZvoX1Zoma+sneMnOTXbu9l40rXinrN/kzeayLWfnJIGgZI3F0fpLQuwv0BNVHuhB7a+zSVw08DEj0iABLxIwJNEHib305a4SkHR20R36RtK2gRKRZdAqU0/42D+W1D12TuwdeyCtAvyUxIgARLwAgFPEnk5uY+2wreKrhPVrMSneQIhIyldvQRKbPwZ8fmroerQ7vhlwsKMivMcCZAACWSbQF4PrtQEXt+DciySikCZGiXhWP8zDse0Q5XB3bAtzyzUGtQuVSm+JQESIIGcEfCkRZ6zOwT41WVrlUKxDQtxJF8VVH6yK7a/r3/cUEiABEjAewSYyL3HMt2aytQKQ9F1C3E0X2VEDeqC7RMXp1uWJ0iABEggqwSYyLNKLJvly9YpjcJrF+FYvihEDpRk/sGSbNbEy0iABEjgRgJM5Dfy8Om7cnVLI3TNIhzPF4nIv3dmMvcpbVZOAoFDgIk8l2Ndrl4ZI5mfCIkwkvm29xfnsgW8HQmQgN0IMJH7IaKazG9aG2e0zCsP6oxf3uUDUD+EgbckAdsQYCL3Uyg1mRdeF2c8AK3yVBdsnSDrnFNIgARIIBsEmMizAc1bl+gD0KLrF+Fo/iqoOrgLtrzNGaDeYst6SCCQCDCR+znaZWqXlnHmi3A4fzSqPdsNm8bM87NFvD0JkIDVCDCRmyBipWuGoeSWRThYoAZqvHgH1r/KFRBMEBaaQAKWIcBEbpJQlapWEmW2LcS+QnVQZ8SdWNt9FBAZCQRLiPQ4ebJJLKUZJEACZiPARG6iiJSoUhwVti/A7vx10eDHV7HqkKxX5nAAhw7Jlh2yZweTuYmiRVNIwDwEmMjNEwvDkmKRxRBRKhHbUBONZM/r5WjmtPDSJWD4cJNZS3NIgATMQICJ3AxRSGVDkeO7URV7sVn2um6BlVgq/xpy+HCqknxLAiRAAtIDSwgmJBAejptwWdrk26RNHotW0i6PQxtAPqeQAAmQQGoCTOSpiZjh/ejRQKFCKIA/pE2+GavQGG2xGIsKdoEjWfrMKSRAAiTgRoCJ3A2GaV727QtMmgRERCBf0HU0qnQKK8r0RLtdHyCu4XNIvs5kbppY0RASMAEBJnITBCFNEzSZx8cDycnIc/ggmh79Dsvq/B3tNr6NJTUGIOnq9TQv44ckQAKBR4CJ3CIxD84bjBab3sOyVsPRds8krKncB1cu/GER62kmCZCALwkwkfuSrpfrDgoOQsslo7C8x9tofnwqtkZ2x4UTF718F1ZHAiRgNQJM5FaLmNjb4vtnsPqxjxFzfgHib74Vv+48a0EvaDIJkIC3CDCRe4tkLtfTZNLD2DJ8KqIvbcKFui0Rv+xILlvA25EACZiFABO5WSKRDTsajLoT+yfOR+mkYwhp0ww7vtuRjVp4CQmQgNUJMJFbPII1n2iNs9OXIgRJKNu7BdaOW25xj2g+CZBAVgkwkWeVmAnLR/Woi+TlK3EhXxjqPHMr4gZOM6GVNIkESMBXBJjIfUU2l+st2zQKJXetwP6iDdB6Ym/M7/KusXBiLpvB25EACfiBABO5H6D76paFo0qh6qGF2FipB26f8xQW1B6Ca1c4cchXvFkvCZiFABO5WSLhJTvyFS2ImANTsarhIHTYPh5rwnvj9xOyBC6FBEjAtgSYyG0Y2qC8edB07btY3ecdNDs9A4crt8GRdSdt6CldIgESUAJM5Db+HjT56klse20GKl/ZDjRpjC1fbLGxt3SNBAKXABO5zWNf56Xu+HWqDE8MSkKVh5pj8ZCZNveY7pFA4BFgIg+AmEf1ikG+zetwNPQWtBrfA/Pbv8mlcAMg7nQxcAgwkQdIrEvUKo/KR5ZgQ1Rv3L7oeSyLegAJv14OEO/pJgnYmwATub3je4N3+YoVQuy+b7Ci8yi0PPIVjkS0wKHlh28owzckQALWI8BEbr2Y5chiXQq3+ezh2DryB1S6sheFWsVi/bilOaqTF5MACfiXgCeJvJKYGCeqKzLJ8Ac8JUqxOIF6L3fD+flrkRhSHPWeaYeFDYbCER6B1u3aAZGRwOTJFveQ5pNA4BDwJJEnCY5nRGuINhH9e8prOVCsTCC8Q3WEHViLjcXbo/2msVh5pCISHYWAQ4eA/v2ZzK0cXNoeUAQ8SeQnhMjGFCoJctwpWiHlPQ8WJxBaoSgahu7CErSUX+nVOIWy2I8o4JLMBh0+3OLe0XwSCAwCQVl0M1LKLxWtJXpB1F2kCQdVhIWFxUyZMsX9nG1eJyYmIjQ01Db+qCPanRLkcMivdT1UwlEUwiX8IiFuHLQOSxYtspWvdoyfe4DonzsNa71u27btBrE41tdWa/bSG92Z2Y2io6MddpW4uDj7uRYRIWncWCzRcRxlHJLQjfcrQ1o6rp6/aCt/bRk/twjRPzcYFnspeXV9Zrk1vfOedK3otSGiusi1PgGbLkqxE4HRo4FC0jcuUk46V2pJe3xBUAc0vbYMR8o1wpG52+zkLX0hAdsR8CSRa/fLx6LaNz7OdgToENC3LzBpEhARAUdQEEIiKuK2Lx/Ein/MR+jVMyjVuSE2PvYvZ5udvEiABExHwJNE3lysfkC0nejmFO0sR4qdCGgyj4939onLUZN78xEdcHXtVmwu2gYNPnoCWyr3xOXDp+3kNX0hAVsQ8CSR6yaQ2iqvI1ovRefIkRIABMJjSyP21Gz82HYcqsfPxcXKtRH/AcMfAKGnixYi4Ekit5A7NNUXBELyB6PbosHY8OF6/BpUGpF/74JfWj4OR0KiL27HOkmABLJIgIk8i8ACuXiz/6uNUgfWYVrloai5/N84UaYuzkxfGshI6DsJmIIAE7kpwmAdI0pXyo87972FHwYvwZUrQIm72mBv16edE4is4wYtJQFbEWAit1U4c8cZGdiCnuNaImnDVnwX9gSqzn4Hp8rURsJMXZKHQgIkkNsEmMhzm7iN7hdd/yb0PPY+PnsoDgmJQSh8Rzsc7iSTe3/7zUZe0hUSMD8BJnLzx8jUFobIVLF+n7VBwvKt+KTUs6gw72P8Vr4GLn42lePOTR05GmcnAkzkdoqmH32p37wQ+h59Ex/2W4P9l8vhpr/djZMNuxpj0/1oFm9NAgFBgIk8IMKcO07mzy9rHH8qa/6sXoM3y41H6IYluHpzDSQMHwNcvZo7RvAuJBCABJjIAzDovnY5pnFeDD70ND4ZuhNzkzui8Osv4veI2nDMm+/rW7N+EghIAkzkARl23zutfedPvlUJ1XdMx7O15uLXUw4EdeqIhNtk8cwDB3xvAO9AAgFEgIk8gILtD1erVwfe2NIRy/+1DSMLvI6gn39CUvQtuPbsi0CC7lNCIQESyCkBJvKcEuT1mRIIlm/Z3wbkx4BDw/BSr9346vo9CHl7DK5ERMPxkSysef16pnWwAAmQQPoEmMjTZ8MzXiZQujQwYWoFRC75AvdFrcLG81EIeuxRXK3ZAFiwwMt3Y3UkEDgEmMgDJ9am8bRVK+Dz3U2wZuwKPFRwCo7tli6WDh2QdGtHWShZV0qmkAAJZIUAE3lWaLGs1wjow9DBQ4Lw5sHeGPPATgzBWFxYtA6oXx/Jfe8HDh702r1YEQnYnQATud0jbHL/ypQB/vNFfty7Zgj6NNyPMXgBf3w9DcnR1YBBg4CTJ03uAc0jAf8TYCL3fwxogRBo1AiYt7oYqnw7Bq3L78OkpIdxfeK/kFy5CvCijHA5d46cSIAE0iHARJ4OGH6c+wR0VcW77waW7KuAxLc+RMPQXfj6ck8kj/knkiOjgBEjgN9/z33DeEcSMDkBJnKTBygQzStQABg6VAayHLwZG4f8FzF5t+KHi7cCr77qTOivv84x6IH4xaDP6RJgIk8XDU/4m0DJksDYscCMfbUw88FpiA3agJ8SmgHDhzsT+j//CSRyuzl/x4n39z8BJnL/x4AWZEIgIgL49FPgi20NMKn7LDTCGiy8IJ3qw4YhOSIS0ITOWaKZUORpOxNgIrdzdG3mW40awPTpwIcbGuHdjnPQBKuwMMEtoWuXy4ULNvOa7pBA5gSYyDNnxBImI9BAJoL++CPwzuommHDbHKOF/nNCE2eXi7bQR43iQ1GTxYzm+JYAE7lv+bJ2HxJo3BiYPRv4YF0jTOw8G7FYh3kJLYCXX3Z2uegoF24758MIsGqzEGAiN0skaEe2CcTKXhY//AB8tCkWX/SaaTwUnZXQxhjlcj08EnjlFeD8+WzXzwtJwOwEmMjNHiHa5zGBevWAb74Bvt7dAD8+/D0a5t2EHxLaAyNHIqlSpJHQ8/KhqMc8WdA6BJjIrRMrWuohgapVZdr/f4CZh+th3QvT0Dx0C76/2MFI6LG970Pyy9JCZ5eLhzRZzAoEmMitECXamC0C5coBY2S70HnH6+D4hKnoVH4LZl29HcGjRuJKuUhcHv4aR7lkiywvMhsBJnKzRYT2eJ1A4cLAU08Bsw7XwXaZHfpQvS2Ye6UtCr7+DySWjsKpwTIO/eJFr9+XFZJAbhFgIs8t0ryP3wnkyQO0bn0Gn2+qg8iN32NEt/VYdq0JykwYhrMlbsaGhyfij8Q//G4nDSCBrBJgIs8qMZa3BQFZ9hyvzIxB7MnZmPz4cuwLjkbMpwNxsmg1fNt9MvbvTbaFn3QiMAgwkQdGnOllOgTCwoC+HzRHw8TFWDdqPq4VLo57frwfv0fH4vkGC/Dtt8DVq+lczI9JwCQEmMhNEgia4V8CwXmC0HB4B1Q5tx7n3puMysXO4Y1NHVDk3k5oW2YHnnySu9D5N0K8e0YEmMgzosNzgUcgOBglBt6HYid3I/nNt3HrTauw7EId3DJxIG6tfwZ16wLjxgGnTgUeGnpsXgKeJPJPxPxfRbeZ1w1aRgJeJpA/P4KffQYh8fuQ54kBGBD0IY4WrIr7f5+IZ5+5jgoVgM6dgcmTOeDFy+RZXTYIeJLIP5N6O2ajbl5CAtYnUKoU8P77CNq6FQWaxeDZQwNx8ZZYvH/fSuzYAdwv+0TrvqN9+8rwxlnAHxz0Yv2YW9ADTxL5UvGLGyZaMLg02YsEdA3dBQuAKVNQIOEMBnzZHAdvfQwrZ583kvi8eUC3bkDZssCjjwI//QRcu+bF+7MqEsiAgOyS6JFESilpb6BWBqX7yzlVhIWFxUyRL7wdJVF2pAkNDbWja4ZP9C/z0Oa5fBkRn3+OSlOn4lrRotg7aBCON2+LDRtLYNGi0lixohQuXcqLIkWuoUWLMzJ2/TTq1z+PkBBH5pXnsATjl0OAfry8bdu2G+T2sb40IVIq3+bpDaKjox12lbi4OLu6ZvhF/7IQ3o0bHY6YGIekZ4ejWzeH4/hx4+LLlx2OGTMcjr59HY7QUOfpokUdjvvvdzimT3c4EhOzcI8sFmX8sgjMRMUlv673NMemLudJ10rqa/ieBEhACeisotWrnRuLardLzZrAV1+hQH4H7rgD+O9/gdOnnZtg9OwJzJkD3HknoN3uev4TGUbA0S/8KnmDABO5NyiyjsAlkDcvMGSIc5B5tWrOp569egFnzxpMChQAunZ17jl68iSk60X6H6UDcvNm4JFHAF3Yq4lsbjR6NCDPU6VpH7go6Xn2CXiSyL+W6leJyrcUR0Xl60chARK4gYAm8eXLnRtB6z50deoAcXE3FAkJAdq2lS3q3gHi44FNmwDdxChZVgN46SUYY9QrVXIm+u+/58KMN8DjmwwJeJLI+0gN0m6AfA1RUfRjUQoJkEBqAroq1/PPO7tb9IF4+/bGPqJpDV8JkmEGuhGG7EqHtWuB48flP5b8z9LWuW6OoV0wJUsCbdo4fxs2bnQm/NS35HsSUAKeJHKSIgESyAoB3R16gwxA+NvfgNdfB9q1A06cyLAG7WJ5+GHgu++AM2ecjfmhQ517SA8bBsTEOMer33uvM+Fri55CAi4CTOQuEjySgDcJaItcm9g69VOb0/pgdOlSj+6QL5+zJa6bYmj3i/4GfPmlcybpkiXOcepRUUCVKsD//Z+zBa/975TAJcBEHrixp+e5QeC++4A1awAZb260zHWhliw+0dRJRjqDVIauG10w22QgsPaz16rlTOJ9pPNTW/Q6Z2n8+KpG10wmfwDkhue8Ry4SkEfuFBIgAZ8S0Iy7bp2zq+WZZ2RGhmTiDz8EtOmdRdG+dR3lqKorMiYlOVvt+lxVdcGCMpg501npzTcDLVs6tUULQN/r9RT7EWAit19M6ZEZCRQpAshMUGOYysiRwP79wLRpzkHlObBXRz82bOjU554DFi5cIY3/1tAuGO3J+eEH59BHvUXp0kDz5k5t1gzQrnxZG4xiAwJM5DYIIl2wCAFZItdI5NWrO1vnjRsDc+cC0dFecyBPHgdiY2Wet6g2/nVo465dwLJlkKUDgJUrAR3aqKJ/EGjXfdOmgJqiGhnJVrsBx2L/MJFbLGA01wYEtFNbn1Z27+5sHuuUT21W+0D0t0P7zlX1waiKPhhdtco5SlKP2sszYYLznLba1ZRGjf5s6etMVIq5CTCRmzs+tM6uBHTAuDaRO3RwzhLSZvJtt+WKt/rwVJcMUFXRVRp/+cX5TFafy+q4dv1tcT2TDQ93tvB1CKSqdsnoFnkU8xBgIjdPLGhJoBGoWtXZ19GpE9Cli3OM4T335DoFnXGqyVn18cedt79wwTlqcv16WclJVIfFT5/+p2kVZWqgdsvopCbXkd0yf/LJ7VdM5LlNnPcjAXcCOm5Qn0xqN4sOVdTmsY419LPos9k2bZzqMuW335xrxOiweE3sOsZ99uw/Z5zqNboygW6Hp1q7tnOIpI1XfXah8fuRidzvIaABAU9Ax5hrX4Ym8wcfdI4p7NfPdFiKFftrcr90ydkto4uAbdniVB3vLsv2/08qV3YmdB2FqapDJ/X5ri4oRvEOASZy73BkLSSQMwI33eTcK65HD+dcfR0grlsNmVwKFfpzxIvLVB0pc+iQczVHXdFx+3ZnstffKnVLRR/C6rh214NYPd5yi6zMVw1QFJSsEWAizxovliYB3xEoWNA58FtXzHrsMWeT1QTdLFl1WJO0DspR1XXXXaL7me7eDWOvU03uqjt3On+/XAley+oKkDpCU1UTu0t1w2tK2gSYyNPmwk9JwD8EtL9Bnyrqw0/tXilc+MZs6B+rvHJXHbeu/eaq7qIJft8+Z1LXMe+qmuA//fTGLhpFU758rPFwVZ8Tu1Rb9joSJ5BnrTKRu3+j+JoEzEBAM9aMGc7hiHff7ew/1yVxbSqa4F1dLO4u6vBHXTNGW/F79jh11aorMlQy1Jix6t6K1y4eTei6kJiq9su7jjp8MhurIbibYvrXTOSmDxENDEgC2hLXTuXWrZ0tct1aSGfpBJBoC7t8eafqhhwqixdvk9E0bYy+9sOHgb17na15bdGramtesV296iyv/2pXjw6XdHX3uI6RkTBmsuo9dCl5KwsTuZWjR9vtTaBECeCnn5yzP7t1c07F1CxEga4xo61u1dtvvxGIPmzVjTp0OZuDB52qr+PjnTj1nLtoXdovHxHxp2or3qV6Tlv8ZhYmcjNHh7aRgI4z1yamLoii/eY6G7R4cXLJgICrBa6tcP2DJrVcueIcVaMjazS5q+pr1YULnT8C+mPgLrpbkyZ0Va1X1fVaH8Lqe38meyZy92jxNQmYkYAO33D1md91FzBvnv07fX0YB30E4RoJk9ZtdE6Wttq160aT+5EjTtX3qvpbeu7cX6/Ucfaa1FW1u8Z1dHUP6W+yPpTVmbTeFiZybxNlfSTgCwLatPzkE+CBB5yrX+nrQB6m4QvGKXVqonV1s+h67mmJToQ6ehQ4dsx5dL3W96o6tFIXJ7t+/a9X6yJkmtRdqsldX+dEmMhzQo/XkkBuEtAx5fp0T9cz19WrBg7MzbvzXm4EtBtFZ6dmtAKxJvHTp52JXUffqGpL3/VajzrMUhO+/hWQE2Eizwk9XksCuU3glVeci5wMHuxc0CS9JmNu28X7/YWAjoTR1rZqRqL98efP52yPkeCMbsBzJEACJiOgT/J0J2YdrtGrl/PvepOZSHOyRkBDqg9TcyJM5Dmhx2tJwB8EdJEtXb9cO2r14af7oGl/2MN7+p0AE7nfQ0ADSCAbBHQqpC4zqLtADBuWjQp4iZ0IMJHbKZr0JbAI6OJa+sBz/HjnwuCB5T29dSPARO4Ggy9JwHIE3nrL+dCzXz/n8AjLOUCDvUGAidwbFFkHCfiLgM5u+fZbZ3+5Dk9Ma+Cyv2zjfXONABN5rqHmjUjARwR0muLEibqiFMK//tpHN2G1ZibARG7m6NA2EvCUwEMPAbLkbaQ+ANU91ygBRYCJPKDCTWdtS0Cn60urPEmXv9X+ct2tgRIwBJjIAybUdNT2BGQRj91Dhji3uh892vbu0sE/CTCR/8mCr0jA8gTOtmjhXFhLE/mGDZb3hw54RoCJ3DNOLEUC1iHwzjtAmTKA9puzi8U6ccuBpUzkOYDHS0nAlAR044l//9u5lurYsaY0kUZ5l4Cnibyj3Ha36D7RF7xrAmsjARLwOoGuXZ3rsOiStwcOeL16VmguAp4kclmMETJIFZ1EZYEH9Ek5yoFCAiRgWgLaxaK7JDzxBKBb0lNsS8CTRN5IvNeWuP6s65imb0TvEKWQAAmYmYDuNaYPPefPB6ZMMbOltC2HBGTwaaYiix5Du1YeTSkpe02hsWjq7Un6y2eqKrVEtxmv7PePbNSEM/Zz638e0b//obDkC8bPkmEzjJYpupCJAL6RXlLtR25VayJ/3+19Wi/Xp/WhTT6zs28aIvpn7S8q42fd+GU7dp50rRwTLpXc2FSU1/oZhQRIgARIwAQEPEnk68TOqqJRovlE7xWdKUohARIgARIwAQEdkZKZyNagkK27MVl0kOh/RaeJZiYbMitg4fN29k3DQv8s/OVk/CwdPLv/37N0cGg8CZAACZAACZAACZAACZAACZAACZAACaQikNm0/fxSXvagMiYTrZFjpKiVJDP/+okzp0U3p6hrnL28Nb18Ihb+KrotHUt1fsG7ojoRbKtoA1ErSWb+tRFnfhd1xe4fVnJObNVRZHGiO0S3iz4lmlqsGkNPfGsjzlo1fgXE9rWiW0Q1diNEU0uu5U59SLpftLKojmRRo2qIussT8ubDlA90pIsmdauIJ/71E2cyG09vVn9biWGanNNL5J3l3FxRTQZNRPWH2EqSmX9txJlZVnIola3l5L3rx7WwvN4jmvr/n1Vj6IlvbcRfq8ZP/0+FiqqEiOr/Lf0/5i5Zzp2eDD90v4HrtSfT9nUa/+cpF3wnx/ai6oQVxBP/rOBHejYulRPn0jspn2vsvhDVBTpWixYT1f9gVpHM/LOKH+nZeUJObEw5mSDHnaIyH/8GsWoMPfHtBkct9kb/TyWm2KyJXDX1QjhZzp3ZTeT6pTmSYowejoqm/iK5l0mS8/qnUElRK4i77WpvWv7p53eJateD/lC5T5qSt5YWT/23spNNxXj9S1L/8qhpYUcixfb6oqn/arJDDNPzTdyFleOnf/Frt96vogtEM4qdR7kzu4lc7h3w8qMQiBStI6rBcP31IS8pJiegrdkI0bqi74nOELWi6J/o00SfFr1gRQcysDkj36wev+vidz1RnSWvf/3XEs2RZDeRH5O7urdA1SD9zF3cy+SVE0VFz7oXMPFrd9vVzLT8U1+upvjwkRxjUl7b4eCJ/1b2U5Oe68/bOfI6RLSUxRxSmzWJTxadnobtVo5hZr7ZIX4ast9E40R1YIW7uMfOp7lTKz8gGiXqetiZ+s/Tv8s594edU+S9VcQT/9z7jHuKY9qXbCWJFGO3pWNwF/nc/WHn2nTKmfnjSDEuPf/KyjnX8xptER12ey8vTS9q+xeiEzKw1Kox9MQ3K8cvTGJWLCVuBeW4TLRrynvXIVdzZ2e5qz4t19Erw1MsGCnH7imvdZjNVNF9opoIKotaSTLzb4w4o8OHtJ9Vf1Wri1pFvhZD9aHSNdGjoo+IDkhRORhJbqIcNba/iMaKWkky82+gOOOKnf4AN7OSc2JrC1F9QKbPZzanqH5f7RBDT3yzcvzqSJw2iWrstKHhGvpqp9wpblFIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgARIgAR8SeD/AZ8IeEHSk87QAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1 簡単な行列計算(25点)\n",
    "\n",
    "関数\n",
    "$$\n",
    "f(x) = \\frac{4}{1+x^2}\n",
    "$$\n",
    "を多項式\n",
    "$$\n",
    "F(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + a_4 x^4\n",
    "$$\n",
    "で補間することを試みる．\n",
    "与関数のx=[0, 0.25, 0.5, 0.75, 1.0]\n",
    "での値から直接逆行列から多項式補間する手法を試す．\n",
    "連立方程式の係数行列$\\mathbf{A}$(ヴァンデルモンド行列と呼ばれる)およびデータベクトル$\\mathbf{y}$は\n",
    "```\n",
    "A = [[1.         0.         0.         0.         0.        ]\n",
    " [1.         0.25       0.0625     0.015625   0.00390625]\n",
    " [1.         0.5        0.25       0.125      0.0625    ]\n",
    " [1.         0.75       0.5625     0.421875   0.31640625]\n",
    " [1.         1.         1.         1.         1.        ]]\n",
    "y = [4.         3.76470588 3.2        2.56       2.        ]\n",
    "```\n",
    "となる．ヴァンデルモンド行列$\\mathbf{A}$の逆行列をデータベクトル$\\mathbf{y}$に掛けることで，\n",
    "係数の値を求めよ．\n",
    "\n",
    "以下は$\\mathbf{A}$，$\\mathbf{y}$を求めるコードと，与関数，補間関数のプロットである．\n",
    "```python\n",
    "import scipy.linalg as linalg   # SciPy Linear Algebra Library\n",
    "import numpy as np\n",
    "\n",
    "nn = 5 # x_i number\n",
    "\n",
    "def func(x):\n",
    "    return 4.0/(1+x**2)\n",
    "xx = []\n",
    "yy = []\n",
    "for x in np.linspace(0,1,nn,endpoint=True):\n",
    "    xx.append(x)\n",
    "    yy.append(func(x))\n",
    "\n",
    "print(xx)\n",
    "print(yy)\n",
    "\n",
    "a_matrix = []\n",
    "y_vector = []\n",
    "for i in range(nn):\n",
    "    for j in range(nn):\n",
    "        a_matrix.append(xx[i]**j)\n",
    "    y_vector.append(yy[i])\n",
    "        \n",
    "A = np.array(a_matrix).reshape(nn,nn)\n",
    "y = np.array(y_vector)\n",
    "print(A)\n",
    "print(y)\n",
    "\n",
    "inv_A = linalg.inv(A)\n",
    "print(np.dot(inv_A,y))\n",
    "```\n",
    "\n",
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 0.25, 0.5, 0.75, 1.0]\n",
      "[4.0, 3.764705882352941, 3.2, 2.56, 2.0]\n",
      "[[1.         0.         0.         0.         0.        ]\n",
      " [1.         0.25       0.0625     0.015625   0.00390625]\n",
      " [1.         0.5        0.25       0.125      0.0625    ]\n",
      " [1.         0.75       0.5625     0.421875   0.31640625]\n",
      " [1.         1.         1.         1.         1.        ]]\n",
      "[4.         3.76470588 3.2        2.56       2.        ]\n",
      "[ 4.          0.15529412 -5.39294118  4.29176471 -1.05411765]\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg as linalg   # SciPy Linear Algebra Library\n",
    "import numpy as np\n",
    "\n",
    "nn = 5 # x_i number\n",
    "\n",
    "def func(x):\n",
    "    return 4.0/(1+x**2)\n",
    "xx = []\n",
    "yy = []\n",
    "for x in np.linspace(0,1,nn,endpoint=True):\n",
    "    xx.append(x)\n",
    "    yy.append(func(x))\n",
    "\n",
    "print(xx)\n",
    "print(yy)\n",
    "\n",
    "a_matrix = []\n",
    "y_vector = []\n",
    "for i in range(nn):\n",
    "    for j in range(nn):\n",
    "        a_matrix.append(xx[i]**j)\n",
    "    y_vector.append(yy[i])\n",
    "        \n",
    "A = np.array(a_matrix).reshape(nn,nn)\n",
    "y = np.array(y_vector)\n",
    "print(A)\n",
    "print(y)\n",
    "\n",
    "inv_A = linalg.inv(A)\n",
    "print(np.dot(inv_A,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAppElEQVR4nO3dd3hUVf7H8fcJBCIGAhI6JAGkSE8BQpMEFQuCihVjYV3l57qWFXUtrIKs6OoqoK4uy9pYF1dBiiBFQRJKkBaadEUDRIoUhYRIzfn9cRNBBDIJydy5k8/ree6Tmcydyfdw9ZOTc88911hrERER7wpxuwARETk3CnIREY9TkIuIeJyCXETE4xTkIiIepyAXEfG48r7sZIzJBLKB48Axa21CaRYlIiK+8ynI8yVba/eUWiUiIlIsGloREfE448uVncaY74AfAQv8y1o7+jT7DAAGAISFhcVHRUWVcKmBIS8vj5CQ4P39p/Z5m9rnXZs2bdpjra1RnPf6GuT1rLXfG2NqArOAB6y18860f7NmzezGjRuLU0/AS0tLIykpye0ySo3a521qn3cZYzKKe/7Rp19t1trv87/+AEwCOhTnh4mISMkrNMiNMecbYyoXPAZ6AmtKuzAREfGNL7NWagGTjDEF+39grZ1ZqlWJiIjPCg1ya+23QFs/1CIifnD06FGysrI4dOiQ26UUWUREBOvXr3e7jHMSFhZG/fr1CQ0NLbHPLMo8chEJAllZWVSuXJmYmBjy/9L2jOzsbCpXrux2GcVmrWXv3r1kZWXRsGHDEvvc4JzHIyJndOjQIapXr+65EA8GxhiqV69e4n8NKchFyiCFuHtK499eQS4i4nEKchEJWFdddRU//fTTWfd55plnmD17drE+Py0tjauvvrpY7w0kOtkpIgHHWou1lunTpxe679ChQ/1QUWBTj1xEXDF8+HBatWpFq1atGDlyJJmZmTRr1ow77riDVq1asW3bNmJiYtizx1l09a9//StxcXF07dqVfv368fLLLwPQv39/Pv74YwBiYmIYPHgwcXFxtG7dmg0bNgCwZMkSOnXqRGxsLJ07dybYlhBRj1ykDPvTn2DlypL9zHbtYOTIs++TkZHBu+++y+LFi7HW0rFjR7p3787XX3/NmDFjSExM/NX+S5cuZcKECSxcuJCwsDDi4uKIj48/7WdHRkayfPly3nzzTV5++WXeeustmjdvzvz58ylfvjyzZ8/mqaeeYsKECSXT4ACgIBcRv1uwYAHXXXcd559/PgB9+/Zl/vz5REdH/ybEAdLT07nmmmsICwujcuXK9O7d+4yf3bdvXwDi4+OZOHEiAPv37+fOO+/k66+/xhjD0aNHS6FV7lGQi5RhhfWc/a0g2M9FxYoVAShXrhzHjh0D4OmnnyY5OZlJkyaRmZkZdCsoaoxcRPyuW7duTJ48mdzcXA4ePMikSZPo1q3bGffv0qULU6dO5dChQ+Tk5PDpp58W6eft37+fevXqAfDee++dS+kBSUEuIn4XFxdH//796dChAx07duTuu++mWrVqZ9y/ffv29OnTh06dOnHllVfSunVrIiIifP55f/7zn3nyySeJjY39pZceVAqm+ZTk1rRpUxusUlNT3S6hVKl93uZL+9atW1f6hZSC7Oxse+DAAXvw4EEbHx9vMzIy3C6p2E53DIBltpiZqzFyEfGEAQMGsGbNGo4cOcKdd95JXFyc2yUFDAW5iHjCBx984PnVD0uLxshFRDxOQS4i4nEKchERj1OQi4h4nIJcRPyuc+fOhe4zcuRIcnNzS72W9957j/vvv/+s+6SlpbFw4cIif/bJi36VJgW5iPidL6FYnCA/fvx4cUs6q+IGub8oyEXE78LDwwEnIJOSkrjhhhto3rw5KSkpWGt57bXX2L59O8nJySQnJwPw+eefc8kllxAXF8eNN95ITk4O4PR6H3/8ceLi4hg/fjxJSUk89NBDtGvXjlatWrFkyRIA9u3bx7XXXkubNm1ITExk9erVv6lr6tSpdOzYkdjYWC699FJ27dpFZmYmo0aNYsSIEbRr14758+eze/durr/+etq3b0/79u1JT08HYO/evfTs2ZOWLVty991341znU/o0j1ykLHNrHduTrFixgrVr11K3bl26dOlCeno6Dz74IMOHDyc1NZXIyEj27NnDc889x5QpU6hduzYvvvgiw4cP55lnngGgevXqLF++HIBRo0aRm5vLypUrmTdvHnfddRdr1qxh8ODBxMbGMnnyZObMmcMdd9zBylPa3rVrVxYtWoQxhrfeeouXXnqJV155hXvvvZfw8HAeffRRAG699VYefvhhunbtytatW7n88stZv349zz77LF27duWZZ55h2rRpvP322yXxL1ooBbmIuKpDhw7Ur18fgHbt2pGZmUnXrl1/tc+iRYtYt24dPXv2JCQkhCNHjtCpU6dfXr/55pt/tX+/fv0AuPjiizlw4AA//fQTCxYs+GUN8h49erB3714OHDjwq/dlZWVx8803s2PHDo4cOULDhg1PW/Ps2bNZt27dL88PHDhATk4O8+bN+2Xp3F69ep11/ZiSpCAXKcsCYB3bgmVn4ddLz57MWstll13G6NGjT3tl56nL3556p3pf71z/wAMPMHDgQPr06UNaWhpDhgw57X55eXksWrSIsLAwnz63tGmMXEQCUuXKlcnOzgYgMTGR9PR0Nm/eDMDBgwfZtGnTGd/70UcfAc4NLCIiIoiIiKBbt26MHTsWcMbmIyMjqVKlyq/ed/Jyt2PGjDltLQA9e/bk9ddf/+V5wRDNxRdfzAcffADAjBkz+PHHH4vV9qJSkItIQBowYABXXHEFycnJ1KhRg/fee4+77rqLNm3a0KlTp1/ux3k6YWFhxMbGcu+99/4yTj1kyBAyMjJo06YNTzzxxK+CusCQIUO48cYbiY+PJzIy8pfv9+7dm0mTJv1ysvO1115j2bJltGnThhYtWjBq1CgABg8ezLx582jZsiUTJ04kKiqqhP9VzqC4yyaebdMytt6l9nlbMC9ja621Bw4cKHSf7t2726VLl/qhmuIr6WVs1SMXEfE4newUkaCSlpbmdgl+px65SBlk/XShivxWafzbK8hFypiwsDD27t2rMHeBtZa9e/eW+LRFDa2IlDH169cnKyuL3bt3u11KkR06dChg5m4XV1hY2C8XQJUUBblIGRMaGnrGKxYDXVpaGrGxsW6XEXA0tCIi4nE+B7kxppwxZoUx5tPC9q28aRPExED+VVQiIlJ6itIjfwhY7/PeW7bAgAEKcxGRUubTGLkxpj7QCxgGDCxs/2OUZzdVMbmWQ7c/ycGn3yevfjQVWzSi7nWJhHVrD5UqnWPpIiICYHyZgmSM+Rh4AagMPGqtvfo0+wwABgA0pXL8a3TiKOUJ5RjVzT6i7BZq4pwlP0p5tlRryQ8JiZibEjncOBp8XJ3MbTk5Ob8sih+M1D5vU/u8Kzk5OcNam1Cc9xYa5MaYq4GrrLX3GWOSOEOQnyzBGLus4El0NMc3Z5KVBRvS97Jz8iJCFqXTOGseiXYhIVi2RzQnt08/Gr1wDyH16hSnHX5TcEeTYKX2eZva513GmGIHuS9j5F2APsaYTOBDoIcx5r8+fXqlSjBsGOXKQXQ0XH5rde4c14vbtz5P7MEFfPbOdt7r8Cabc2pz4fuDOV4/inXt+vHjtMC9N56ISKApNMittU9aa+tba2OAW4A51trbCv3k6GgYPRpSUk778nnnwZW/q03/xX+gw8FUpo3YxOR691Nv1XSqXd2FjfV6sPPjBUVsjohI2VMq88izmzaFzMwzhvipKlaEXn9qwo1ZI9i+5HvGdRpBxPZ11L6xG2vr92T7lGWFf4iISBlVpCC31qYVNj5+ri5qH85NC//EsY3fMrHzy9T4fiW1r+nAsra/J2fzrtL80SIinhSwV3bWb1qJvumPcGTtN8xo8ShtVr9PXpOmLO43grwjv72nn4hIWRWwQV6gfosq9Fr7Eus/+oq1EZ3p+OFANkZ2ZsujrztXj4aE6CpSESnTAj7IC7S9qRkd90xnzoAPqZH9HXVeeYT5Wxpw1IboKlIRKdM8E+QAIeUMPf51M8dr12MJHejGAjbRlC1EQW4uDBrkdokiIn7nqSAvUGvXarqSzkISqcsOLmAfc7kYu2WL26WJiPidJ4OcqCgAOrOIQ1RkM43ozjwWlOvO/u9zXC5ORMS/vBnkw4b9suhWHXbRhtXMD+lO5+Pz2d2wA+sn+r5Io4iI13kzyFNSnKtGo53FtkKio+n2n3tYN3IWVY/vocH17Zlz7zi3qxQR8Qvv3uotJeU3V462BvYkrWBL95vo8a+b+WzJcnosHEZoWDl3ahQR8QNv9sjPIrJtPZrvSGVRu3u5fMWLLK13LXu+PeB2WSIipSboghyg3HkVSFzxT5b0f5MO+2bwY/NENkzb7HZZIiKlIiiDvECHd//A5n/OIvLYLiJ7J7J4hJbHFZHgE9RBDtDs3mQOz11EbmhV2g7swewBOgkqIsEl6IMcoHa3Jlyw8Us2V03g0n/fzIykF8k7Xvgt7kREvKBMBDlAeEwkzbbNZmnjW7hy7hPMbvkQRw/nuV2WiMg5KzNBDlA+PIyEjWNZ1Plhem58nYUx/cjZe9jtskREzkmZCnIAUy6ExPThLL7h73TfOY4NDa9g77f73S5LRKTYylyQF+g4/lGWP/w+bbMXsLNlD3as3u12SSIixVJmgxwgbvhtbHhhMo0OreNgwsVkzt/mdkkiIkVWpoMcoPUTvdgy+nNqHttO+aSubJy6ye2SRESKpMwHOUDze7qxZ3waYfxM1WsuZu2HX7ldkoiIzxTk+RpdH8vPM+dhQ8pR59YkVr29zO2SRER8oiA/SYPLmpOXNp+D5aoQc/clLH893e2SREQKpSA/Rd2ujai4eD77KtSm2YM9WfpSqtsliYiclYL8NGrG1afKinnsCGtIq8evYumwz90uSUTkjBTkZ1C9RS2qr0pl23lNaf2XPiwePN3tkkRETktBfhbVmtag5ldz+K5SS2KHXsviiJ5079EDYmJg7Fi3yxMRARTkharauDp1n7+fr2lC7IE0ltj2sGULDBigMBeRgKAg90HEiGepzzY20ZQ4lrOIDpCbC4MGuV2aiIiC3CdbtxJBNlFsYQPNiS8I861b3a5MRERB7pOoKACqkEM0maznIuJZzpLKl7hcmIiIgtw3w4ZBpUqAE+YxfMd6LnLGzB+f4HJxIlLWKch9kZICo0dDdDTWGKpEVyf61YGsD29P7Eu3sPQJhbmIuEdB7quUFMjMZO6cOZCZScSD/YleN5P14e1p9+ItLHtqotsVikgZVWiQG2PCjDFLjDGrjDFrjTHP+qMwL4hoUIWotTNZf34CbV+4mYynJ7tdkoiUQb70yA8DPay1bYF2wBXGmMRSrcpDqkZVocHamWyoFE/r525i+ZApbpckImVMoUFuHTn5T0PzN1uqVXlMtegI6q/9jI2VYmn17A2s+OunbpckImWIsbbwTDbGlAMygAuBN6y1j59mnwHAAIAaNWrEjxs3roRLDQw5OTmEh4ef9rWD3/9Mg9//hWaH1zDt9y9zwW2t/VzduTtb+4KB2udtwdy+5OTkDGttQrHebK31eQOqAqlAq7Pt17RpUxusUlNTz/r6nq/32bVhcfYQFezKv83wT1ElqLD2eZ3a523B3D5gmS1CHp+8FWnWirX2p/wgv6JYvzXKgOoXVqPmqllsDmtJsyeuZfXLWgJXREqXL7NWahhjquY/Pg+4DNhQynV5WmTTC4hcMZvvKl5Ek8euYfUrs9wuSUSCmC898jpAqjFmNbAUmGWt1dm8QtRsfgEXLJ9NZsVmNHm0D1+N/MLtkkQkSJUvbAdr7Wog1g+1BJ1aLapjl81ma3wPGj/cmzXlPqXVAz3cLktEgoyu7CxltVtFUjXjC7ZVaEyjB69m7T90D1ARKVkKcj+o1aoGEUu/IKtCIxo+0Iu1b6S5XZKIBBEFuZ/UblOTykvm8H2FhsTc34u1b851uyQRCRIKcj+q07Ym4YvnsL1CDDF/vEphLiIlQkHuZ3Xa1SJ88Rx2hEYT88erWPOPNLdLEhGPU5C7oE67Wpy/JJXtFWJo9MBVfPWaToCKSPEpyF1Sp10tKi9NJatCIxo/1IvVI+e4XZKIeJSC3EW129QkYtkcsio2psnDvVj1sq4AFZGiU5C7rFbrmlTNmMPWik1p9lhvVrww0+2SRMRjFOQBoGbLGlRfNYfvwlrQ4qlrWDZEKyCIiO8U5AEisll1aq35gm8qtaHNs31Z0uc5iImBkBDn69ixbpcoIgFKQR5ALmhcjXprZ7GxYlvipg7hyy11wFrYsgUGDFCYi8hpKcgDTNWYqkRH5rCGlnRgCQvo7LyQmwuDBrlbnIgEJAV5AKqyfSNN+JqVtKMrC5lHV+eFrVvdLUxEApKCPBBFRXE+P9OSNSwhgYtZQCpJEBXldmUiEoAU5IFo2DCoVIkwjtCOlXxJR5JJY855vbB5hd8sW0TKFgV5IEpJgdGjITqaCuY4HRrsIr3WdfTY8Cap7f9M3nGFuYicoCAPVCkpkJkJeXmU2/odnbI+Zn6bP9Jj+cvMbXEvxw4fd7tCEQkQCnKPCCkfQtcVrzP/4kEkbxrN4kb9OHTgiNtliUgAUJB7iAkxdJv7HAuufZku28ezOqYPB3YcdLssEXGZgtyDuk56hEX3vE38j7PIvPBSfli/1+2SRMRFCnKPShx9F6sGjadp7goOtO1G5vxtbpckIi5RkHtY3HN92fzGZ9Q89j2hSZ1Z9/E6t0sSERcoyD2u5X3d2TtxHqEco/aNXVkyfIHbJYmInynIg0DDa9uSt2AhByrUoM0jl5J6/wS3SxIRP1KQB4nanRpSfUM6myPi6P7GjXzW6zWsrhsSKRMU5EGkcsNImmz5guUNruXy6Q8xq/VAjh7ShUMiwU5BHmQqRJxH/Lfj+bL9A/RcO4LFUTeyf0eu22WJSClSkAchU74cnZa8xqJ+r9J592S2Nkpi29KdbpclIqVEQR7EEj94kDV/nUyjQ2shsSOr/rPK7ZJEpBQoyINcm7/04Yfx8wg1x2h8ZxfSBk5xuyQRKWEK8jKg4Q3xVFi5lKzwi7h4xLV8dslLWgpXJIgoyMuIC1rVpdG2uWQ0vJHL5zzO/Ia3k/3Dz26XJSIlQEFehlSoWomEbz4k/arn6LbtA7ZFd2XLAt0HVMTrFORljAkxdJk2iNVDP6HBoa+pdHECy4bPc7ssETkHhQa5MaaBMSbVGLPOGLPWGPOQPwqT0tXu6d78+NkSckKr0e6RHnwR9yg2KpruPXpATAyMHet2iSLiI1965MeAR6y1LYBE4I/GmBalW5b4Q1TP5tT4dgnLq13CJSteYeG2+uTYSrBlCwwYoDAX8YhCg9xau8Nauzz/cTawHqhX2oWJf4TXi6B9+Abm0o1EFrGL2mymIeTmwqBBbpcnIj4wtggrKxljYoB5QCtr7YFTXhsADACoUaNG/Lhx40qwzMCRk5NDeHi422WUqO49emCsZTntaEAWlcjlK1rR0Sxl7pw5bpdXooLx+J1M7fOu5OTkDGttQnHe63OQG2PCgbnAMGvtxLPt26xZM7tx48bi1BPw0tLSSEpKcruMkhUT4wynADuoxU7qEMtKvgztRvwPM6lQtZK79ZWgoDx+J1H7vMsYU+wg92nWijEmFJgAjC0sxMWDhg2DSk5Y12EXrfiKWaYnnY7OZ1udDmybscblAkXkbHyZtWKAt4H11trhpV+S+F1KCoweDdHRWGMIja7PZe/fQfoznxF+eA+RV7Vn+T3/RAuciwQmX3rkXYDbgR7GmJX521WlXJf4W0oKZGY6Y+KZmZCSQpdne3J4yWpWRiQR99Z9rGp0HT9v3e12pSJyCl9mrSyw1hprbRtrbbv8bbo/ihP3RSXUJGHXNKYmD6d55gwONmpN5ps6/CKBRFd2SqFCK4bQe87DZIxaxg+mJjF/7MVX3f6Azc5xuzQRQUEuRdD5/1oT+e1SJjR6lJYL/sWOWm3ZM1GX94u4TUEuRVKzQUX6fvN3Pnl4LocOwQXXJ/H11X9yLiASEVcoyKXIjIHrhnfjWMZqPq5xH02mvcquWq3JnpLqdmkiZZKCXIqtaez5XPf9P3jvzlSycwyVr+nB1isHwE8/uV2aSJmiIJdzEhoK/d9LInvBat6JfIx6M9/mp7otOPjeeM07F/ETBbmUiNgulUjJeolR/Rez+ec6nP+7m9jZ/mpnTrqIlCoFuZSYihXhj+8mwKLFvFRnBOEZczl8YQuyB70Ahw+7XZ5I0FKQS4mL71ieh7f8iXceXc+MvCuo/PxT7I9ujZ35mduliQQlBbmUitBQePDvDWi+biKPtZrBD7ss5soryL6sL3z7rdvliQQVBbmUqubN4cVVV7Dgn2sYGvY8ZvbnHGt6EUcfewqys90uTyQoKMil1IWEwO/urci9W57kLzds5IPjNxP68gscim6KfettOH7c7RJFPE1BLn5TsyaMHF+PmLn/4daGX7L8x4aYe+7mcMs4mDXL7fJEPEtBLn538cUwZmMii19J587zxvH9xmzo2ZNjl14BK1e6XZ6I5yjIxRWhofDwQMNL393IC7evZyCvcGDOUoiNJS/lNvjuO7dLFPEMBbm4qlYt+Pd/KnLL4oH0a7+ZF3iCI/+bQF7TZvDAA7Bzp9sligQ8BbkEhA4dYOaiqjT+6AW61/2G0cfu4vgb/ySvUWN46inYt8/tEkUCloJcAoYxcNNNMPebeuT8fRTtwzfwv5+vI++Fv5EX0xCefRb273e7TJGAoyCXgBMWBo8+CrO+u5DlA/9LfPnVfHLwUhgyxAn055/XHHSRkyjIJWBVrw6vvAKTv2nFlDsmkGAy+Dy7Mwwa5AT63/4GObrdnIiCXAJedDS8+y78Z00co/t8SgcW88WBDvDkk+RFxziBrh66lGEKcvGMFi1g4kQYldGB166YTiJf8kX2SYH+/PNw4IDbZYr4nYJcPCcuDqZOhVcXJTLysul0YDGzsxOdIZfoGHjuOZ0UlTJFQS6e1bEjTJsGby7twBtXTSOBpczM7gpPP+0E+rPP6rZzUiYoyMXzEhLgk0/grRUJ/OeGKSSYDD7NToIhQzgeFQODB8OPP7pdpkipUZBL0GjXDj78EP63MY6pd02iffkVfJJ9CQwdyrEGMTB4MOV1UlSCkIJcgk6TJvDvf8OUre1Y+sQEuoSvYtLBnjB0KAk33kre04M15CJBRUEuQatOHXjhBZi5vQ3bR47nyrqr+PTw5YQ8N5RDdWL4edBfNctFgoKCXIJe5crw0EPw6dY2rB0yhDvbrWLGoWTOe/4Zcmo2ZNfDf4ODB90uU6TYFORSZpQrB92772HMijbELJ/Es72XMf9oIrVGPsneCy4k4643OJJzxO0yRYpMQS5lUmwsDJ4ST8LOaYz9wwK+CWlK/Lv3szOiGR/1Gcvmr/PcLlHEZwpyKdNq1ICUN7vQPieNpc99xtHK1bh56m3sb5rA43Gz+OgjOHzY7SpFzk5BLgKElDO0H9STxvuWse/1sTSquo8XV/Skyi1XklxrHQ8+qLvQSeBSkIucLCSEC+6/lao7N5L30stcev6XzD/QhoveuJ9LY/fQti0MHw67drldqMgJhQa5MeYdY8wPxpg1/ihIJCBUrEjIY48QmvkN5e67l3vNKLLOa8Jt+9/gsUeOU68eXHUVjB2rCS/iPl965O8BV5RyHSKBKTIS/vEPzOrVhHWO57Et93PwogT+cetC1q2D225z7juakgKffgpHNOlFXFBokFtr5wG6YaKUbS1awKxZMG4cYdl7uPf9Lnx36T0snPYjKSkwcyb07g21a8Pdd8Pnn8PRo24XLWWFsdYWvpMxMcCn1tpWZ9lnADAAoEaNGvHjxo0rqRoDSk5ODuHh4W6XUWrUvsKV+/lnoseMocH48RyNiODrBx5ge5dkMpZfwJw5NUlPjyQ3tzxVqhyla9c9dO++m9jYHwkNLfz/tXOl4+ddycnJGdbahGK92Vpb6AbEAGt82ddaS9OmTW2wSk1NdbuEUqX2FcHy5dbGx1sL1vbube327dZaa3/+2drJk61NSbE2PNx5OSLC2ttus3biRGtzckquhFPp+HkXsMz6mLGnbpq1IlJcsbGwaJFzY9FZs6BlS/jgA8IqWq65Bv77X9i927kJxnXXwfTp0LevM+x+zTXwzjua/SIlQ0Euci7Kl4eBA51J5s2aOWc9b7gB9u4FICwMrr7auefozp0wZw4MGODs/vvfOwt7JSbCsGGwejX4MNIp8hu+TD/8H/Al0MwYk2WM+X3plyXiMc2awYIFzo2gp06FNm0gNfVXu4SGQnIyvPoqZGbCihXOTYzy8uAvf4G2baFBAyfoJ03SwoziO19mrfSz1tax1oZaa+tba9/2R2EinlOuHDz+uDPcEh4Ol1wCgwaddvqKMc6NMJ5+GpYsge3b4e23nd75hx86QzDVq0NSkvO7YflyJ/BFTkdDKyIlLS4OMjLgd7+D55+HHj1gx46zvqVOHbjrLvj4Y9izx+nMP/qocw/pJ5+E+HhnvvottziBn5npn6aINyjIRUpDeLiTuGPHOt3p2FiYN8+nt1ao4PTEX3jBGX7ZsQPef9+5knTuXGeeesOG0Lgx/N//OT34nTtLtzkS2BTkIqXp1lth8WKIiHB65sOHF/mMZu3azhWkY8Y4QzBr1jjj7K1aOSHer5/To2/RAkaMaMKHHxb6B4AEmfJuFyAS9Fq1gqVLnaGWRx5xknjUKKfrXUTGOLMcW7aEBx+EY8ecXntqqrPNmlWLKVOcfS+8ELp1c7auXZ3nxpRw2yQgKMhF/KFKFRg/3pmmMnQobN4MEyY4k8rPQfny0L69s/35z/DFF+lERHRn7lxnJOeTT5ypjwA1a0KXLs7WubMzlF+xYgm0TVynIBfxl5AQJ8ibN3d65x07wowZ0LRpif2IcuUsCQmQkOB0/vPyYMMGmD8f0tNh4UJnaiM4fxDExkKnTk4pHTtCTIx67V6kIBfxt379nLOVffo43ePp050udSkICXHGzlu0cE6MgnNi9MsvnVmSX37pjPKMHOm8VrOmU0qHDid6+uf4R4P4gYJcxA2JiU4XuWdP5yqhSZPgssv88qNr13aWDLjuOuf50aPw1VfOOdnFi5157dOnnzgnGxXl9PDj450tLs65RZ4EDgW5iFuaNHHGOq68Enr1cuYY3nyz38sIDXXCOS4O/vAH53sHDjizJpctc7aMDJg48cR76td3hmXatTvxVcMy7lGQi7ipTh1ncnifPs5UxaNHnbmGLqtSxZnLnpR04ns//eSsEbN8uRPsK1bAtGknrjitUsVZmaBtW2dr3dqZsBOkq84GFAW5iNsiIpyxjD594I47nDmF/fu7XdVvVK3623DPzXWGZVauhFWrnG3MGMjJObFPo0ZOoBdsLVs653fDwvxbfzBTkIsEgvPPd+4Vd+21zrX6x445l3AGuEqVTsx4KZCXB1u2OKs5rl4Na9c6YT99utMscE7CXnjhiROxLVrARRc5a4+df747bfEyBblIoDjvPGfid9++cM89Tpc1AIZZiiokxJmU07Chs+56gSNHYONGWLfOCfe1a2H9euf3V0HAg7MCZPPmztas2YmtXj3/t8UrFOQigSQszDmr2KuXM7xSufKv09DDKlRwxs1bt/71948cgW++cUJ9wwZnW7/euZDp5CGasDCoWzeB2FjnPHHBduGFzkycsnyiVUEuEmjCwmDyZGc64k03OWMSl1zidlWlpkKFE8MrJ7PWWTNm40bYtMnZvvzyEF99Fc4nn/y6F1+pkhPojRs7W6NGJ75GRRVrNQRPUZCLBKLKlZ0A797d6ZHPmeNcpVOGGAN16zpbcrLzvbS0NSQlJXHsGGzdCl9/7fTmC7YNG5x/tsOHT3xOSIgzXbJguKdgi4lxtrp1naXkvUxBLhKoLrgAPv/cufqzd2/nUsyGDd2uKiCUL+/0ths1gssv//VreXnOKpGbN8N33znb5s3OGu6ff+68dupnNWgA0dEntqioE1uDBk6PP5ApyEUCWZ06ThezUydn3Dw9HapVc7uqgFbQA69f3/mD5lSHDjmzarZsccI9M/PE8y++cIL+1LsxVa/uBHqDBic+u+BxvXrOVzfDXkEuEuiaNz8xZn799TBzZvAP+paisLATM2FO5+hRJ8y3bnXCfds2Z9u61dnS02Hfvt++r2pVJ9Tr1XOGawq+Fmx16jgnZUNDS75NCnIRL+jeHd55B26/3Vn96p13yvY0jVIUGnpiiKVbt9Pvk5sLWVnw/ffO14LHBdvatc7iZMeP//a9kZFOqBdstWs7X8+FglzEK267zTm7N3Sos3rV/fe7XVGZVamSc3Xq2VYgPn4cdu92gn3HDmfbvv3E4x07nGmWO3ee9v7cRaIgF/GSwYOdRU4efthZ0ORMXUZxXblyTm+7du2z75eXBz/+eG7LBeuenSJeEhLirJLYqBHccIPzN714WkiIczL1nD6jZEoREb+JiHDWL8/NdU5+njxpWsokBbmIF7Vo4SwzuGQJPPmk29WIyxTkIl7Vt69zwnPECGdhcCmzFOQiXvb3vzsnPfv3d6ZHSJmkIBfxsrAw+OgjZ7z8tttOP3FZgp6CXMTrmjWDN96AtDSi/vc/t6sRFyjIRYLBnXfCTTcRM2aMc781KVMU5CLBwBh44w2OVa7sjJcfOeJ2ReJHCnKRYBEZycaBA507IQ8b5nY14kcKcpEgsrdrV2dhrWHDICPD7XLETxTkIsHm1VehVi1n3FxDLGWCglwk2FSrBv/6l7OW6iuvuF2N+IFPQW6MucIYs9EY840x5onSLkpEztHVVzvrsAwdCt9+63Y1UsoKDXJjTDngDeBKoAXQzxjT4uzvEhHXvfqqc5eE++5zbkkvQcuXHnkH4Btr7bfW2iPAh8A1pVuWiJyzevWck56ffQbjxrldjZQiX24sUQ/YdtLzLKDjqTsZYwYAA/KfHjbGrDn38gJSJLDH7SJKkdrnbadv3y23OJv3BfPxO8NdRAtXYncIstaOBkYDGGOWWWsTSuqzA0kwtw3UPq9T+7zLGLOsuO/1ZWjle6DBSc/r539PREQCgC9BvhRoYoxpaIypANwCTCndskRExFeFDq1Ya48ZY+4HPgPKAe9Ya9cW8rbRJVFcgArmtoHa53Vqn3cVu23GalqSiIin6cpOERGPU5CLiHhcsYO8sMv2jTEVjTEf5b++2BgTc06V+pkP7etvjNltjFmZv93tRp3FYYx5xxjzw5nm+hvHa/ltX22MifN3jefCh/YlGWP2n3TsnvF3jefCGNPAGJNqjFlnjFlrjHnoNPt48hj62DbPHj9jTJgxZokxZlV++549zT5Fz05rbZE3nJOem4FGQAVgFdDilH3uA0blP74F+Kg4P8uNzcf29Qf+4XatxWzfxUAcsOYMr18FzAAMkAgsdrvmEm5fEvCp23WeQ/vqAHH5jysDm07z36cnj6GPbfPs8cs/HuH5j0OBxUDiKfsUOTuL2yP35bL9a4Ax+Y8/Bi4xxphi/jx/C+plCay184B9Z9nlGuA/1rEIqGqMqeOf6s6dD+3zNGvtDmvt8vzH2cB6nCuwT+bJY+hj2zwr/3jk5D8Nzd9OnXFS5OwsbpCf7rL9U/+xf9nHWnsM2A9UL+bP8zdf2gdwff6frR8bYxqc5nWv8rX9XtYp/8/bGcaYlm4XU1z5f3bH4vTsTub5Y3iWtoGHj58xppwxZiXwAzDLWnvGY+drdupkZ/FNBWKstW2AWZz4DSqBbzkQba1tC7wOTHa3nOIxxoQDE4A/WWsPuF1PSSqkbZ4+ftba49badjhXyXcwxrQ6188sbpD7ctn+L/sYY8oDEcDeYv48fyu0fdbavdbaw/lP3wLi/VSbPwT1sgzW2gMFf95aa6cDocaYSJfLKhJjTChO0I211k48zS6ePYaFtS0Yjh+AtfYnIBW44pSXipydxQ1yXy7bnwLcmf/4BmCOzR+994BC23fKeGMfnLG8YDEFuCN/5kMisN9au8PtokqKMaZ2wZijMaYDzv8HXulkkF/728B6a+3wM+zmyWPoS9u8fPyMMTWMMVXzH58HXAZsOGW3ImdnsVY/tGe4bN8YMxRYZq2dgnMw3jfGfINz4skza2j62L4HjTF9gGM47evvWsFFZIz5H86Z/0hjTBYwGOekC9baUcB0nFkP3wC5wO/cqbR4fGjfDcAfjDHHgJ+BWzzUyQDoAtwOfJU/1grwFBAFnj+GvrTNy8evDjDGODfsCQHGWWs/Pdfs1CX6IiIep5OdIiIepyAXEfE4BbmIiMcpyEVEPE5BLiLicQpyERGPU5CLiHjc/wN//VgOdzre2gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "    \n",
    "for i in range(0,5):\n",
    "    plt.plot(xx[i],yy[i],'o',color='r')\n",
    "\n",
    "# Lagrange補間\n",
    "f = interpolate.lagrange(xx,yy)\n",
    "# print(f)\n",
    "x = np.linspace(0,3, 100)\n",
    "y0=func(x)\n",
    "y1 = f(x)\n",
    "plt.plot(x, y0, color = 'b', label=\"original\")\n",
    "plt.plot(x, y1, color = 'r', label=\"interpolated\")\n",
    "plt.xlim(0,3)\n",
    "plt.ylim(0,5)\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAAXNSR0IArs4c6QAAAERlWElmTU0AKgAAAAgAAYdpAAQAAAABAAAAGgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAABiqADAAQAAAABAAABBgAAAABgCvtYAAA5sUlEQVR4Ae2dB3hUZfbGX3qREgihaOgBhAQIBEKHIChlV3QVV1DBSlFYRd1VFAurICp/K6KCgGBlsSOKIJjQe4cgCIhIpJdgQGrmf85MRidhmEwyd2buvXnP87y5/Su/7zKH7/vuPRegkQAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkICVCBSyUmH9LWtkZKSjVq1a/p6e7byTJ0/isssuy7bPqhusi/lajm1ivjbREtmlXQKtx5o1aw4LjihztpJxpbpWkpoYExPjyK8lJyfn91LTXce6mK5JHGwT87WJlsgu7RJoPeT3c7W3n+PC3nZaeN/XUvaB5cuXt3AVWHQSIAESMBcBuzkKc9FlaUiABEjABgToKGzQiKwCCZAACQSTQNFgJs60SYAESMBMBM6dO4e9e/fi9OnT2Yqlw9Vbt27Nts+KG/7Wo2TJkoiOjkaxYsX8qiYdhV+YeBIJkIAdCKiTKFu2LPSpyEKF/nro8/fff3fut3od/amHzN3jyJEjTodZu3Ztv6pshaEnfVZ1mugd0a1+1YonkQAJkIAXAtqTkMfnszkJL6fZepc6SGWQs1flq9LhchRTpFAHRZtzFK67bG8T7RANzzp2gyw/FQ0Q9craxwUJkAAJ5IuAZ08iXwnY4KK8MgiXo5gqrNUpeFoR2Rgv6iFqJOqbtYyW5a8itQuuRXD+zhm9GinPHELm+czgZMBUSYAESMCCBMI1R7FQWNXKwStRtrUnsStr/3RZXifaK1JnsV7ky7ENlOMq59hbSkqKrubJprx1FjPSbsL3Easx7Km9iEqMyNP1Zjs5IyMD+eFgtnpoeexSF7vUw6ptopO9Oo6f0y5cuOB1f87zzL6dl3ro0JMVfh9qCXTPoafesj3JoyH6yfobIp2jeFf0lsivOYqEhIR8vf6ZeSHTMbrHDEdEoWOOkjjleKFHsuPcH+fylZYZLgr0LU0z1MFdBrvUxS710HaxYl1SU1Pdt1S25YkTJ7JtW3XjUvWQp70uqpI3FvIba9k3s09K4e8U3Sv6UOTLnCE80tPTfZ1zyWOFChdC20eikLr2DHpU24BHZyehVcWfsGGGTpvQSIAESCBwArt378aVV16JW2+9FQ0bNkTv3r1x6tQpzJ8/H82aNUPjxo1x11134cyZM1i1ahVuuEGnaYGvvvoKpUqVwtmzZ50T0XXq1HHu37lzJ7p37w75DzK6deuGH3/80bn/jjvuwODBg9GqVSs88sgjzn35/ROuoSdv5U2TndU9Duhwk+4LuVWLr4LP9lbGZ/9ZhiGvxKDFzREYPi4FT3zTBiXKlQh5eZghCZBAEAgMGyYD2jqiDZSSoScU0WnSAC0+Hnj11VwT2bZtGyZPnox27do5ncLLL7+MCRMmOJ1F/fr10b9/f7z11lsYOnSoFNFVxkWLFiEuLs7pPM6fP+90AJrRwIED8fbbb6NevXr44YcfcN999zmXekwfB166dKlULbC6+Rrz13xCaasks3oifbC3uKiPaKYoL2ZYrCftXfR+qQ22/lQUt9ZdgVGLkxAfJdAnbMpLeXguCZAACVxEoHr16k4noQduu+02p4PQdxrUSajdfvvtWLhwIYoWLYq6des6XwZcuXIlHnroIed+dRodOnRwzt2pI7jpppsQL05qmDi/ffv2OdPQP7o/UCeh6YSrR/Gx5J0kqiTSyeqnRZNFQ0VzROr+9BHaLaK8mA49XZvfoSdvGVWsWwFTd7RHX3kiauDTVdF+cG38660FGP1dAspULePtEu4jARKwAgGP//n/EeIX7nI+nhoREeF8Cc4bto4dO2L27NnOt6i7du0KHVLSSeuxY8ciMzMTeq2715HzhTujPpkQrh6FPvpaTVRMpENM6iTUvhWpS60rGi3KqxnWo8iZcbcRLbB5T3kMabwI4zZ0QFz0ccwdsybnadwmARIggVwJ7NmzB8uWLXOe99FHH6FFixbQuYsdO3Y4973//vvo1KmTc117Dq+KU2vTpg2ioqKcDkWHrnQYqly5ctCeyCeffOI8V2assWHDBue6kX/C5SiMrINnWgFNZnsm5G297OVlMW5jJywcvxkli5xFt8cTcGe9RTj283Fvp3MfCZAACXgl0KBBA4wfP945mX3s2DE8+OCDePfdd51DRTqZXbhwYedEtF6sk9EHDhyA9izUmjRp4pzwdvdKPvzwQ+d8R9OmTZGYmOic9HaeaOCfcA09GViFbElpj+JreVZ6QLa9Bm+0v68J1t9yGs/0TMGLy9rju5gjGP/wctzwYmuDc2JyJEACdiSgcw8ffPBBtqp16dIF69aty7ZPN/RJJ30Cym0TJ050rzqX2qP47rvvnOueQ09Tp07Ndl4gG3brUQTCIk/XlowoieeWJmHVRztQtcQx3Di2NXpHL8P+jRqZhEYCJEAC9iFgN0cR1KEnb83erO+VWHm4Lp67JgWz0pqhUXwxTBuwGI5Mh7fTuY8ESKCAE9DItZs3e75rbH4gdnMUQZvM9tWUxUoXw2NzkrD+m9/QqMyvuGNSe/SovAa/LNnr6zIeIwESIAFLELCbowgr9Ct71sHCo3EY13sBFh+5ErHtI/DGTQsYZDCsrcLMSYAEAiVgN0cR8qGnnA1QuGhhDP2kE7YsPo72kT/iX592QseKm7FttjvWYc4ruE0CJEAC5iZgN0cRlqEnb01cs100Zh9MwNR7FiM1ozqa9rwcY7ql4Nypc95O5z4SIAESMC0BuzkKU4HWMCC3v9MeWzeex7VXrMPjcyXIYKUdWPfxj6YqJwtDAiRAAr4I0FH4omPQsSpxUfhkbxsJMrgc+85URMtbYvB42xScPn7aoByYDAmQgBUJ6JvUGobDSNOAgUab3RxF2OcofDWQvpCXuqME+tdbhjHLJMhgld+w+M2Nvi7hMRIgAZsR0FAd+ma2RojVMBzPPvssWrZs6Xzj+umnNeydy9577z3nPn3jul+/fs6dGufp008/dZ+CMmVc8eb0A0Qa6uPmm29Go0b6gVBjraixyYU9tZC8mR1ILSvUjsCU7R3Q9/k1GPhkZXQY0gRD3l6AMd81h4YIoZEACYSGgEeUcQmyVyqUUcbx008/Ydq0aZAPDTl/+DUyrPYuevXq5YwOGxkZiVGjRjlDhFeqVAlHjx7NFcratWuxfPlyZ3iPXE/O4wl261HksfrhO/3q4QnY9GsF3N90Ad7cJEEGa6RDv9lNIwESsD+BmjVronXr1pg7d65T+sGi5s2bOz86pE5EvyuhIcLVSahVrFgxVyga50lf5guG2a1HEQxGQUtTw5S/tr4TbpZvXNx9f2l0f6IF+k9ZjFfmxkLDm9NIgASCR8Ajyrh8L/sPlC0buh69O/y39iIee+wxDBo0KFtFx40bl23bvaExotxzGrrUr925zZ2me9vIJXsURtLMZ1ptBzXGugNXYES7FHy0qxUa1juPTx92hSDOZ5K8jARIwAIE9NOlU6ZMcX6ASIublpaGgwcP4qqrrnKGDj9y5IizFu6hJ+0xrFmzxrlv5syZkG9hO9eD/cdujuJaATbRyA8XBbsB3OlrkEH9it6q6bsQXfIwbnq5DW68Qp6SWn/AfQqXJEACNiNwzTXX4JZbbnF+a0LDi+v3szUCbGxsLEaMGOH8JoVOZuuX7dQGDBiABQsWQPfp9yyC2YuwGeqLqyMfGZceXf4sOTk5fxcaeNW5P845nu+e7CiBPxwRhY45pty50JF5ITPPOZihLnku9CUusEtd7FIPbSYr1iU1NdXrHSaTyl73W21nXurhjYX8mnqdKLVbj+Jir2HBPUVLFsWjs5Ow8bt9aFz2F9z1bgd0i1qL3YsZZNCCzckik4DlCdBRmLgJ63erjZQjjfFmn4VYdrQ+4jpE4PUbF+DC2QsmLjWLRgIkYDcCdBQmb1ENMnjvxx2xZekJdIzaigc+74QOkanYOmunyUvO4pGAOQnIcJI5CxbCUuWVAR1FCBsnkKxqtLkC3+xvgfcGLca2k1cg/tpojL6aQQYDYcprCx6BkiVLQp8kyusPpZ1Iad2VgbLw1/gehb+kTHCeBhns93Z7dPvXIdzfYw2emJeEGZHbMOWdTCTc1tAEJWQRSMDcBKKjo7F3714cOnQoW0FPnz6dpx/ObBebaMPfeqiTUBb+mhUcRR2pzAhReVFvfytm5/Mqx0Zh+p4o9H18Be59oSYS+0Xh3+NSMHJ2K5SqWMrOVWfdSCAgAsWKFUPt2rUvSkNjJenb0Va3YNUj2ENPUwT8QVHOD8R2l33bRDtEw0W+TL/4c7evEwrqseuea4XUXaVwV4MleHFlEppW3Y8Fr60vqDhYbxIggSARCLajmCrlVqfgaUVkY7yoh0jDHPbNWjaW5awcqizbNB8EImqWxzs/dsS8F9figqMwkobF497YhTix94SPq3iIBEiABPwnUMj/U/N9Zi25Uh1AXFYKbWQ5UtQta/uxrOWYrOWlFp/KAV9DTwPluApVqlRJmD59uq7m2TIyMv4M3Zvni8N8wZljZ/HZ8DOYvL0nqhXej4f+mYzmg/wfhwxz8X1mb+V28ayYXeqhdWJdPFvWHOuBtknnzp01PkiLcNSmlmS62SNj/bGf5LGtgdbf8NjOuRopO94W6fOgbqeS8xz3tjOER0xMjEzs58+s+LZpzpoun7TJEVtiuzwD6HDcWmux49CPh3OeYrltO7SLQrdLPVgXc/4TCvT+kh9Sy76ZfUQKP1hUV5Rbr8M038yWsobNWt0dh7WHa+KBpl9ixu6WaNjQgen3L4Ujk8+Ph61RmDEJWJhAsOcovKFJk53VPQ7o2IjuM8KcPQorBgU0ovKeaRQvUxzXvxqBNZ/uRu3SB9B3XFtcf8VKpK3e53ka10mABEggVwLhcBSrpFT1RPqMWnFRH9FMkRHGHkUOio1vrC/hP67E//09Bd/vb4xGLUvjnf4L2bvIwYmbJEAClyYQbEfxsWS9TNRApBHt9DHX86KhojmiraIZoi0iI4w9Ci8UixQvgoe/liCD8w6hecQuDHy/I7pErsfOH37xcjZ3kQAJkEB2AsF2FH0lu2qiYiIdYposUvtWVF+k8w6jRUYZexQ+SMZ0qYn5h5piwq0LseZ4HTTuEoWXr0thkEEfzHiIBEgACLajCDVj9ihyIa5BBgd+IEEGV55Cl8qb8fDMJLStuBWbv/gplyt5mARIoKASsJujYI/Czzs5umU1zNzXEh8NXYpdp6qi+Q018d/OKTib8dc3eP1MiqeRAAnYnIDdHAV7FHm4YTXIoD4NlSozRDfVXIWRKUlIqPQLVk1LzUMqPJUESMDuBOzmKNijyMcdG9WwEj7c3Q5fP7kSx86VQes7GuDfLVJw6vCpfKTGS0iABOxGwG6Owm7tE9L6/P2ZRGz5uTQGNFyCl9YkoUm1Q0h5lUEGQ9oIzIwETEjAbo6CQ08B3mTla5TH26kdkfyKy0F0fjAegxouRPqe9ABT5uUkQAJWJWA3R8GhJ4PuRI1Cu3GffOdChqAm/dgOjWqfcg5NGZQ8kyEBErAQAbs5CguhN39RS1cqjbGrkrB86jZEFv8dvUYl4pZaS3Bo62HzF54lJAESMIwAHYVhKO2bUMvbG2H1oVrOx2c//UWCDMYWwkdDljAMiH2bnDUjgWwE7OYoOEeRrXmN29Agg0/9kIR1X+5BTOl9uPXNduhVbRX2rtpnXCZMiQRIwJQE7OYoOEcR5Nss9roYLDna0Bn6Y/7BODRKvMwZEiTzfGaQc2byJEAC4SJgN0cRLo4FKl8NMvjgl0nYnHwYLSvswOCPJMhg1AbsmM8ggwXqRmBlCwwBOooC09TGV7ROUg3MO9xMwpYvwloNMti1sjOc+fnTGiCYRgIkYBcCdBR2ackw1UPDgNwzrQNSV53CNVU34j/fSJDBStuw6bPtYSoRsyUBEjCaAB2F0UQLaHpXtKiGL9MSnZ9c3X2qMpr3ro2nO6XgzIkzBZQIq00C9iFgN0fBp57CeG9q7+Lm19pi67Yi6FN7BZ5ZmITmlX/F8kmbw1gqZk0CJBAoAbs5Cj71FOgdYcD1kfUq4v1d7fHNyFU4cb402g5ohIcSUnDy4EkDUmcSJEACoSZgN0cRan7MzweBnk+3xJbdZTA4djFeWZuExlccwfyxa31cwUMkQAJmJEBHYcZWsVGZykWXw5ubO2LB6xtQtNAFdH2kOQZcuRDHf2GQQRs1M6ticwJ0FDZvYLNUr+O/mmLD/qp4JDEFU7ZJkME6f+Crx1eYpXgsBwmQgA8CdBQ+4PCQsQRKVSyFF1YkYcW0bYgqno7rx7TCzTWW4sDmQ8ZmxNRIgAQMJWAVR3G91Pod0f9E1xhKgImFnECL/hJk8EgdPNslBV/+moBGTYrgg3sZZDDkDcEMScBPAqFwFFOkLAdFOZ+R7C77tol2iIaLfNmXcnCAaLDoZl8n8pg1CBQrXQxPzEvC+q/3osFlaej3djv8repq7FmWZo0KsJQkUIAIhMJRTBWe6hQ8rYhsjBf1EDUS9c1aNpblrByqLNtue0JW9DqaTQg0/HtdLDrSCK/dsAALDjVCbFuZ/O6zAAwyaJMGZjVsQaBQiGpRS/JRBxCXlV8bWY4UdcvafixrOSZrmXOh5Xxe9L1oXs6DWdsDZalClSpVEqZPn561O2+LjIwMlClTJm8XmfRsq9Xl2Kbf8foTUUg50RptLluDB5/6FVGJEU66VqvLpW4Ju9RD68e6XKqVw7c/0Dbp3LnzGil9i3DVoJZk7Dn01Fu2J3kUpp+sv+GxnXP1ftmhFXhbpMNPPi0hIcGRX0tOTs7vpaa7zop1ybyQ6Xj37kWOiELHHCXwh+P57smOc3+cc1ixLt5uCLvUQ+vGunhr4fDuC7RN5Id1tbcf11AMPXnLN6/7XpcLEkTqJNRZXMoYwuNSZCyyX8OA3DGpPbauP4u/Xb4ew79LQquKPyHth6MWqQGLSQL2IxAuR6EzltU9cEbLOmcxPYAU9NWqTSrjs7TW+PTfy5B2JhK3P9sLI9ql4PTx0wUdDetPAiEnEC5HsUpqWk9UW1Rc1Ec0UxSoMdZToARNdv2NY9sgdXsx3Hh5Mp5bmoRmVdKwdMImk5WSxSEBexMIhaP4WBAuEzUQ7RXdLTovGiqaI9oqmiHaIgrUOPQUKEETXl+xbgXc+2ExfDdqNf64UALtB8fi/qYLkLE/w4SlZZFIwH4EQuEo9NHXaqJiIh1imixS+1ZUX1RXNFpkhLFHYQRFk6bRbUQLbN4bgaFNFuGNjR0QF30cc8foMw40EiCBYBIIhaMIZvlzps0eRU4iNtsuU7UMXt/QCYve3IySRc6i2+MJuLPeIhz7+bjNasrqkIB5CNjNUbBHYZ57K6glaXdvE6w/cDkeb5uC93e0QaOYM/j8keVBzZOJk0BBJWA3R1FQ27FA1rtkREmMXpKE1dN3olqJo7hxbGv0jl6G/Rs1YgyNBEjAKAJ2cxQcejLqzrBQOvE3N8CKwzEY0y0Fs9KaoVF8MUy9ZzEcmQ4L1YJFJQHzErCbo+DQk3nvtaCWTIMM6st5G779DbFl9+DOye3RvfIa7F6sD9rRSIAEAiFgN0cRCAteawMCDXrUwYIjjfHGTQuw9EgDxHWIwLjeDDJog6ZlFcJIwG6OgkNPYbyZzJJ14aKFMWRGJ2xenI72kT/i/s86oWPFzfjx211mKSLLQQKWImA3R8GhJ0vdfsEtbM120Zh9MAHTBixGakZ1NP3bFXjumhScO3UuuBkzdRKwGQG7OQqbNQ+rEygBDTLYf6IEGdx4HtdFr8WI75OQWGkn1n6oAQFoJEAC/hCgo/CHEs+xPIEqcVGY8Wsb57sW+89UQOJt9fBYmxT8cYxBBi3fuKxA0AnYzVFwjiLot4y1M/jHC62RuqMEbq+3DM8vT0J8lX1Y/OZGa1eKpSeBIBOwm6PgHEWQbxg7JF+hdgQmb++A719Yi7OOougwpInEj1qA33/73Q7VYx1IwHACdnMUhgNigvYl0PWR5tj0awUMi0/Bm5s6ILbGCcx+1usHvuwLgTUjAT8I0FH4AYmn2JeABhl8ZV0SlkzYgjJFTqPnUy3Qv+4SHPmJX9Szb6uzZnklQEeRV2I835YE2gxsjHWHovFk+2R8vCsRDRtkYsaDyxgGxJatzUrllYDdHAUns/N6B/D8PwmUKFcCzyzqjDUzdqFGqYO4+dU2+McVK/Hb2v1/nsMVEiiIBOzmKDiZXRDvYoPr3OSmBlh+pD5e7JmCOfuboFFCKUy+cxF7FwZzZnLWIZAXR9FWqnWLqL+HrFNTlpQE8kCgaMmi+M83Sdg4Zz+alv8Z90ztgKuj1uHnhb/mIRWeSgL2IOCvo3hfqvt/ovailllqIUsaCdiaQL1raiP5cBO81XchVh6NQVyninjthgW4cPaCrevNypGAJ4Ginhs+1tUpNBIxwL8PSDxkTwIaZHDwRx3xtwd+w6Br0zDsi074X+QmTPqwFBr1irFnpVkrEvAg4G+PYrNcU9XjOq6SQIEjUL3V5fhmfwt8cO8SbD95OZpdVx2jujLIYIG7EQpghf11FJWETapojmimh2SVRgIFh4AGGbz1zXZI3ZSJf9RYgyfnJ6FF5C6seV//edBIwJ4E/B16GhnG6jeUvB8QqbOaL3pLRCOBsBKoHBuF6b9Eoe+Ilbj3+Rpo1b8SHn4jBSNnt0KpiqXCWjZmTgJGE/C3R7FAMvam3MozRU7QL93r0JWndZeNbaIdouGeB7ysazzowaJ/itp5Oc5dJBA2AteNTkTqrlK4s8EyvLgyCU2r7seC19aHrTzMmASCQcBfR3FHPjOfKtepU/C0IrIxXtRDpBPkfbOWjWU5K4cqy7ZaL9E3om91g0YCZiIQUbM83vmxA+aNXYcLjsJIGhaPwbELkf7rCTMVk2UhgXwTKOTHlU/JOfVFt/lxrrdTaslOdQBxWQfbyHKkqFvW9mNZyzFZS18LdRZ/u8QJA2W/ClWqVEmYPn36JU7zvTsjIwNlypTxfZJFjrIuoW+os8fO4LPhZzFpe09ULXwAw29bgsZ3Rv1ZELbJnyhMtWKXdgm0Hp07d14jDdMir40zUS74WORvz8Nb+rVkp+fQU2/ZnuRxYj9Zf8NjO+dqkux4XTRBNETky5whPGJiYhz5teTk5PxearrrWJfwNcmKyZsccSW2yePkDkffmkscB7cedhaGbRK+NvGVs13aJdB6yI/ram8/sLk5AB0WelaU6e3iEO1LkXzuFw0S6ZCVL2MID190eCxkBBLvisOaQzXx36RkfPpLCzSUQdYPhzLIYMgagBkZSiA3R6H/Q58hqmtgrmmSVnWP9KJlXfcZYc4eRXp6uhFpMQ0SCIhA8bIl8FRyZ6z7fDfqlU7DbePb4LkbS2Dvqn0BpcuLSSDUBHJzFClSoD6iDwws2CpJq56otqi4SNPXdzOMMPYojKDINAwlEPuP+lh8pBFe6fUDlh5vitjE0pjYj0EGDYXMxIJKIDdHoZnr/MKN+SyFzm8sEzUQ7RXdLTovGiqaI9JHX7XHskVkhLFHYQRFpmE4gSIlimLYV1fh/ZfnoUXETgz6oAO6RK7HruRfDM+LCdqLQHKyDOnImM6pU+Grlz+OQkv3Wz6LqHMc1UTFRDrENFmkpo+51hfpkNZokVHGHoVRJJlOUAhUbFYe8w7HY+JtC7H6eF00vqoSXrs+mUEGg0LbHok+9BCwa5e8eLYtfPXJzVE8klW0cbLUJ49yKnwl954zexTeuXCviQgUKlIYA97viNSVJ5FUOVV6Gp3RIXILfpy53USlZFFI4C8CuTkKHRpS00em1niRHjOTsUdhptZgWXwSiG5ZDbP2tcD79y3DtpPRiL+uBsZ0nY/zp876vI4HSSDUBHKL9aQ/vGrTXAvT/9UexbV86sn07cQCZhHQIIP6NNTVQ45gaPcNeHx+F3wamYop71xA09sakxMJmIJAbj0KdyF1PkFfvpsr+sFDsmoqY4/CVM3BwvhLoEqjSHyypxU+Hb4aaWcroUW/KzGi1TycPhrGGUx/C8/zbE/AX0fxiZBYJ3pC9B8PySqNBEjAKAI3jmmB1J0lcVv9VXhuZVc0qXoAKa/oPz0aCYSPgL+OQh9p1fDeK0WecxXhK7n3nHXoaSKHnrzD4V5rEKhYqxze3dYW817aICERCqPzQ81wz5WLcHT3CWtUgKW0HYHcHEVFqbFKh3Q0zpI+6urep0uzGYeezNYiLE++CXR5qCk2/haFRxOTMXVbGzSsewb/e2SNRI/Kd5K8kATyRSA3R6G9B33i6XbRv0VLsrZ1n4pGAiQQRAKlK5XG8ys6Y/UH21CjxAH0GZuAa6PXYs+6I0HMlUmTQHYCuTkKDbNRRyQhzZwRXjfIUr/KMk4UK6KRAAmEgED8rbFYfqQ+Xu4+F8m/NUCj5iXwWr/VuHCe3YsQ4C/wWeTmKNyA9PFY/STp6yJ1Euo4zPjILOcopGFo9iRQpFRxPDj7GmyZk4YO5Tdh2Act0DZqO7bMY5BBe7a4eWrlr6OIkyLfI5KoI04NkKXuM5txjsJsLcLyGE6g1jX18e3hRHzUbzZ2Ha+IhKsr4KXeS3HhXDi/BmB4NZmgiQj46yjWSplbe5S7laxzjsIDCFdJIJQEChUtgr7v9cCWFSfRI2o1/v1ZW3SO2oRd838OZTGYlxcC42TM5cknvRyw8C5/HUWC1HGpaHeWNCJsS9Em0UYRjQRIIAwEKifWwuf722HaXQuwIb0WmnSNwsQb58BxTp9op4WDwP3ymbVRo8KRc/DyzC2Ehzvn7u4VLkmABMxFQMOA9J/cCUmDD+Cunvsx6PNu+LLSEkz6rAIu76rTiTQSCIyAvz2KXyQbXwqsFLyaBEggYAI1WlbB3P1NMO6utUg50QxxV1fF9BtmAGfOBJw2EyjYBPx1FFahxKeerNJSLGdQCBQuUghDJzfH+uVn0KDSUfT94p+4ufIPODKHU4pBAV5AErWbo+BTTwXkxmU1fROo36oCFu2Lwejbt+GLE10Q1/0KzPjbNDgyTvq+kEdJwAsBuzkKL1XkLhIomASKygzk41MbYNXis6ha6Txu/vZ2tK+0FcvH6WfraSTgPwE6Cv9Z8UwSsCSBpu3KYPX+6pj0yHbsOl8Dbe5viT51VmL3xhOWrA8LHXoCdBShZ84cSSDkBIoUAe5+oT5+2lcWT7adj5k/x+HKpsXx6A3bkZ4e8uIwQ4sRoKOwWIOxuCQQCIEyUaXwzJIu2P7NDvSpMBdjv4hBTOUTGD/mBM6dCyRlXmtnAlZxFJdJI+hjG3+3c2OwbiQQKgLRPZtg6oEeWH3fu4g7txZDHy+HJjXTMetrB8OYh6oRLJRPsB3FFGFxULQ5BxN9gW+baIdoeI5j3jYflZ3yQDiNBEjAMALFiqH5+Lvxw+Yq+KrBI8jctx/X9iqEPr1O4vRpw3JhQjYgEGxHMVUY5XyrW0ZLMV7UQ6SvjfbNWjaW5awcqizbV4tSRepwaCRAAgYTKNSoIXptGYPNr8zD6GIjMWPWZejWOA3HjxaMIIMbJQhRRob+LNEuRaDQpQ4YuL+WpKUOIC4rzTayHCnqlrX9WNZyTNYy52K07NChJ3Uqf4j+IfJ2Bw+U/SpUqVIlYfr06bqaZ8vIyECZMmXyfJ0ZL2BdzNcqZm+Tkvv3Y9PjmzDk56dRr+RuPP/8RpRtWskrSLPXxWuhvezs3DkJ9esfw4QJ+rmdwE3TU0tOTnEuA/0zYEACduwoi4kTV6NevQyfyQXaJp07d14jGbTwmUmQDtaSdDd7pN1b1id5bPeT9Tc8ti+1eocc8GuOIiEhwZFfS05Ozu+lpruOdTFdkzgs0SaZmY7v//2dowxOOGrgF0fqQ+84HOfOXQTTEnW5qNQX74BzVubi/fndY3R68fH68VuHY+3a3EsUaJvIb+xqbz/AwR568pZnfvdNlQtn5XIxQ3jkAoiHSSBXAoUKoevYblg49wzOlCiL9i//A8ti7wE2GPM/7lzz5wmmIxAOR5EmFKp7kIiWdd1HIwESMBGBZldXwtLUCoisVhxdtr+Jr5s/DTzxBDjTbaJGClFRwuEoNH5APVFtUXFRH9FMkRHGWE9GUGQaJJBFoE4dYMmGsohrVhzXOz7HpNH7gWbN5Os0+nkaWkEhEGxH8bGA1I8cNRDtFd0tOi8aKpoj2irSx163iIwwDj0ZQZFpkIAHgago4IeFRXFNt8IYINOLz+67B4527RHz+uv6uJDHmVy1K4FgO4q+Aq6aqJhIh5gmi9S+FdUX1RWNFhll7FEYRZLpkIAHAX0QcKb0+/v3B55Kfxj3xS1E1S9kR1wcMHeux5lctSMBiS9pK9MexbXpDF5jq0ZlZcxBQN7Pw9SpwOWXQx6bbY9ZlQ4j5uA2RHb7FZH1UxDZsxUqVS+FyEhkU3WZkSxVyhx1YCnyR8BujkJ7FF+XL19+QP5w8CoSIAFfBOSBKIyRN57qy3jA5Mnncd7REqnba+PIdoeoGC54ubi2zEbu3AnotYGaPii6RQaqtSNDCx2BYA89ha4mzIkESCBkBO68Exg1ajMWLSmM1EOVcWDdPpyLT8QxRGBHtyFY8e0RfCsDzA88APz8s7xItdmYos2eDTRu7HIWxqTIVPwhYDdHwclsf1qd55CA0QTi41Fo1UpEPP8Y6qZMRuKt9dDjwFQ8OEy6AGI//GBMhup81E6ccC35NzQE7OYoOJkdmvuGuZDAxQT0k3qPSvxOfTEvNhaQbkfNQd1Rt+Y5wxzF/PkXZ8s9wSdgN0cRfGLMgQRIwDeBBvI0/IIFEvpTYn/K+xZXpX2AFHnL+/xZbyHafCfleTQtDfjxR889XA8VAbs5Cg49herOYT4k4ItAYflpue8+52TCVbEHcOJ0CaxrIc+YbNVXp/Jn7E3kj5sRV9nNUXDoyYi7gmmQgFEEatRA5+9kOEps/o6agMxl4LnnkJ/P6c2bZ1ShmE5eCdjNUeS1/jyfBEggyASqVC3kfJz1h0RxGNdfD4wYAbRsCaxd63fO+lisOoqKFf2+hCcaSICOwkCYTIoESMA7gauuAhavLIEz7/0P+OIL4MABIDFRvm8pH7j84w/vF3ns1bmJffuALl08dnI1ZATs5ig4RxGyW4cZkYD/BNRRqD9YsUKu0V5Faipwxx3ACy+4hqMWLfKZmHvYiY7CJ6agHbSbo+AcRdBuFSZMAvkn0KkToPPbf05IV6ggny+bBHz/PXD2LNCxIzBkCPD7714zUUehkWz1LW9a6AnYzVGEniBzJAESyJVARATQvLmXF++6dnW9tj1sGPDWW673L/T1aw87fx5ISQH0VFp4CNBRhIc7cyWBAkdAh42WLwdOnsxR9csuA155xfWNi7JlgZ49XWFqjxxxnrh6tetNbA475eAWwk06ihDCZlYkUJAJ6DyF9g4WL74EhdatXU9CPfkk8PHHQKNGwCefYP48VxgQvd5t+hQULXQE7OYoOJkdunuHOZFAngi0aycfpinmZfjJM5USJYBnngHWrJEPJlcH/vlPzHttM5rFnUWlSp4ncj2UBOzmKDiZHcq7h3mRQB4I6AiTdhr+nND2dW2TJs5xqlOjXsbSw/Xlm91va1xzgF0JX9SCdsxujiJooJgwCZBA4AR0nkHfszt2zI+0JMjg4pYP4ixKoGv9PcA997jeu/DjUp5iLAE6CmN5MjUSIAEfBHSeQTsFGjPQH9PHYnW4qv3SF11PRW3b5rrsf/Li3oUL/iTBcwwgQEdhAEQmQQIk4B+BVq1cn0X19/sU6ijatgUuKys/VYMHA+9MdGX0+mviPdq7XtzzL2ueFQABOooA4PFSEiCBvBEoXhzo0MG/eYrDh4H163O8PxFV2ZXhU08DP/0ENGsGPPus66W9vBWFZ+eBAB1FHmDxVBIggcAJ6DyFRvDYv993WsnJrmEqr+9PdOvmSuSGG4CnngJatABWrfKdII/mm4AVHEWS1G6RSB57gK7TSIAELEzA/T6EOgJfpk9H6ft3GmjWbYUKuddkWVl6F/q+xVdfAfpynj5S9cgjwKlTHidx1QgCwXYUU6SQB0U5P63eXfbprNQO0XCRL5OpL2SISor2+jqRx0iABMxPQEeLypfP5X0KqYbOTyQlAfqFVZ/Wq5fzA0m4+25g7FigaVP/Z8t9JsyDbgLBdhRTJSN1Cp5WRDbkG4noIZJXL9E3a9lYlrNySP7L4OxN6LmPiv4ropEACViYQBH5BVAH4GtCe/duYOfOHPMTvuqswaQmykS3dkMy5ZOrmsG997pif/i6jsf8IpCbr/YrER8nLZRjtXIcT5Rt7Unsyto/XZbXicaI/p61z9tCn7wu4e1A1r6BslRh7969EkQsRVfzbBkZGfm+Ns+ZBfkC1iXIgPORPNvEBa169StkxKgepk9fjqpVT19E8ptvqsq+K1Gu3Er59/jXUNKGDRVkf1N5F2OtBJ09cdF1GqK2sHyru/aUKYgWx3Hms8+w/cEHcbRNm4vP/XNPknMtv78Zfybz54qx6WVkJEjKZbFagl6lp+vgyqXNyvdXLamW59BTb9me5FHVfrL+hsd2zlWZrcIEkTw4jSSRL7tWDk6MiYlx5NeSk5Pze6nprmNdTNckDraJq002b9a3KRyOyZO9t1GfPg5HtWoOR2Zm9uNz57quW7w4+36vW8uXOxyxsa4Lbr3V4Th0yOtprilzr4fytdPo9OLjXVVYuzb34gR6f8nvp4RgvNgKX7zLdHs+lxINEt0sShH5Mobw8EWHx0jAJAQ03p/ORXsbftKRIx1B0qedsk1e57Xs+tKGvgb+9NPAjBlAw4aQLozrdzyvaRXw88PhKNKEeXUP7tGyrvuMMGePIj093Yi0mAYJkECQCKgDuOoql6PQ/y972mYZfzh0yKDPnhYvDowc6QoyqF896itTovqFvTSjfnI8S27f9XA4Cn3YuZ5IWg3SiugjmikywtijMIIi0yCBEBBQR6HfwXZH5XBnqU87qXl9f8J1KO9/G8uzMsuWAf/3f66v6mmX5p132Lvwk2SwHcXHUg5pHTQQ6aOt8vwazouGiuaItoqkT4gtIiOMPQojKDINEggBAbcjyBlNVrcbyC+GRhk31PRxq4cfBjZudH1ub+BAg72RoaU1VWLBdhTSz0M1UTGRDjFNFql9K6ovqisaLTLK2KMwiiTTIYEgE9CRoJo1s89T6OezNWCg24kEpQgxMa5JkAkTXENS7kwYZNBN4qJlsB3FRRkGeQd7FEEGzORJwCgC7nkKfUNbJ7DVVqxwfSo16N/HlsdooT2KLR6DGfoI7aZNroLwbzYCdnMU7FFka15ukIC5Ceg8hX6bYsMGVzl12El/w5OSQlTuaB3oyLKffwYSElyT32fOuPdyKQTs5ijYo+BtTQIWIqCOQs09T6ET2fpbXUHfqwu1bZUp05tukvgP/3UVQrs3NCcBuzkK9ih4Y5OAhQhcfrm8f32la57i999dQ0++hp0Ceq8iNy76Ue4PPwS+lp+R48cBHYp66CHXWFhu19r8uN0chc2bi9UjAfsR0F7FwoWuXsV5eSYyqBPZ/uD7u0QS0jjogwYBr7wC6Pe7vb0Z6E9aNjnHbo6CQ082uTFZjYJDQB3FyZMS7G2MhIguCbRrZ4K6lyvn+vRqSopr0kS914ABrp6GCYoX6iLYzVFw6CnUdxDzI4EACSQluUJ1rFzp+rqpOovcLOfb3Lmdn+/jnTq5Ztr/8x9AAg0iNlZeDzbq/eB8lyrkF9rNUYQcIDMkARIIjEBkJBAf70oj7MNO3qpSujTw4ouuCRQt7HUS7LqPBJQ4eNDb2bbcR0dhy2ZlpUjAWgTcTz/5msgORo3y1DPRz61KqG/nN7q/+EK+ptPINfmdp0SCUYvgp2k3R8E5iuDfM8yBBAwnMFSC+jz7rCuyhuGJG5mgBhl84glg3TqJWFcPuO02+YqOTH7/+quRuZguLbs5Cs5RmO4WY4FIIHcCtWq5fn/1ZTtLmPYmFi8GXn1VPn6Q4pq7eOutv14xt0Ql/C+kVZrF/xrxTBIgARIIBQENMvjAA/JZNomLrt++uO8+oHPnUOQc8jzoKEKOnBmSAAnYioBGN5w7V0KeSsxTdywSraC+FGITo6OwSUOyGiRAAmEkoK+M33WX60U9dzFat87uONz7Lbiko7Bgo7HIJEACJiWgMUncphPc+qTUk08CFg8yaDdHwaee3DcplyRAAuEloGFAbrkFGDUKaNbM9YW98JYo37nbzVHwqad83wq8kARIwFAC+nLetGnA7NmuGCUam2TYMEsGGbSbozC0nZkYCZCAuQgENXpssKravbvrySh9Kuq114C4OMD9YfBg5WlwunQUBgNlciRAAsEnYLmXocuWBd54wxUmV1/au/pq4O67LRNkkI4i+Pc0cyABEjApgZA7nA4dXE9CDR/uGpZq2BDQcCAmNzoKkzcQi0cCJGAzAhoeV2Oqa7jcqlWBG24A/vlP4MAB01bUCo5CyzhaNE50u2lJsmAkQAIkkBcCzZu7nMVo+Xn76itAexfvvQeEvJuTe6GD7SgkgDsOiuQd92wmszvYJtohkj6YT7tOjkaLzon2+jyTB0mABEjASgSKFQMef9w1HKWO4nb5v3DPnsCePaaqRbAdxVSprToFT5MAKRgv6iFqJOqbtWwsy1k5VFm2G4iWiuTjtbhXRCMBEiABexHQD4cvWgS8/rprqR9IGi8/k5mZpqhn0SCXYqGkXytHHomyrT2JXVn7p8tSew0yaAeJ13uRaS/ibNbeCxcd/WvHQFlVYe/evRLQMUVX82wZGRn5vjbPmQX5AtYlyIDzkTzbJB/QPC5Zvz5CtuIlyvc6XLiQ7nEkf6uu3+Ek58X5/c24OOcA0mvcGCUnTUL9l15CRYm9fnzCBJxKT5EsKsqnMFYjPT3j4uw89lj5/qol9fAceuot25M86tZP1t/w2M65Wlp2TBbpHMWQnAe9bSckJDjya8nJyfm91HTXsS6maxIH2ySwNpk3TwfwHY6FCwNLx331hQuu9DRNo8w1yRBgapmZDse77zocERGO+ELrnHVeu+JsrokGen/J7+lqb7+pwR568pZnXvedkgvkgWP8S6RDVr6MITx80eExEiABaxDQNwvvuAPYuhUoV85V5v79XR9MCkMNwuEo0qSe1T3qqhPVuo9GAiRAAiTgSUAfn61dx7Xn0CGgZUtgxAjg9GnPs4K+Hg5HsUpqVU9UW1Rc1Ec0U2SEMdaTERSZBgmQgPkIfPYZ0E9G6p97TqZp4oElS0JWxmA7io+lJstE+uSSTkrrENJ50VDRHJH0qzBDtEVkhHHoyQiKTIMESMB8BMqXB959V3455adTexT6lve/ZET+99+DXtZgOwp99LWaqJhIh5h0UlrtW1F9UV3RaJFRxh6FUSSZDgmYkIA7KKAJ30kLHa1rrnEFGZSnopyP0GqQQXUeQbRgO4ogFt1r0uxReMXCnSRAAt4IWNbhlCnz1zsXpeXBUI1QK5PfRU+c8FbNgPfZzVGwRxHwLcEESIAELENAv3Eh75Q4J7g/+ACJ+qTU0qWGF99ujsJwQEyQBEiABExNQIMM6lf05IW8jLoymh8TY3hx7eYoOPRk+C3CBEmABCxBQJ6E2jh2LFBZIx8Za3ZzFBx6Mvb+YGokQAIkALs5CjYpCZAACZCAwQSKGpxeuJPToadr09MDDxYW7oowfxIgARIwCwG79Sg49GSWO4vlIAESsA0BuzkK2zQMK0ICJEACZiFAR2GWlmA5SIAESMCkBDhHYdKGYbFIgARIwCwE7Naj4ByFWe4sloMEgkjAsqE3AmASzjrL1zFsaRK4Hb/ks2aV5LrD+bzWbJexLmZrEYBtYr420RLZpV0CrUdNYRFlziYyV6lWm6s4AZWGdQkIX1AuZpsEBWvAidqlXYJSD7sNPQV8tzABEiABEiCB7AToKLLz4BYJkAAJkEAOAkVybHPTRWCNjUCwLuZrTLaJ+dpES2SXdrFLPcx5l7BUJEACJEACJEACJEACJEACJEACJEACJEACJBA8AvLhWWwT7RAND142QU95t+SwSbReFJTH5STdYNkUSfigaLNHBhVl/XvRT1nLCh7HzLzqrS4jpcBpIm0bVU+R2a26FDBZlCraInpApGbFdrlUXUZKfazWLvJpO6wUbRBpu/xXpFZbtEKkv2P/ExUX0QwioBP7O0V1RApW4TcSWdF2S6H1xRsrWkcpdHORp6N4UbbdjluXL4isYN7qMlIK/m8rFN6jjNVkXdtEraxou0j/bVixXS5Vl5FSH6u1i74wXUakVkykzqG1aIaoj0jtbdG9zrUA/vDx2L/gJcqqeuBdorOi6aLrRLTQElgo2R3NkaW2w7Ssfbq8Psdxs256q4tZy+qrXPvk4NqsE36X5VbRFSIrtsul6pJVPUstHFLajKwSq6NQ6b6rRJ+K1Az590JH4YKpf/XG//WvTezN2uexyzKrerPMFeljcgMtU+pLF7SKHNJ/4Gr7RbptZRsqhd8o0qEpqwyjuXnXkpVmIv3fq9XbxbMuUh1YsV10JESHMA+KdHhWR0WOi86L1Az5HaOjcMG029/2UiEdKughGiLSIRC7mDpBlVXtLSl4XVG8SJ3fSyKrmA5zfCYaJjqRo9BWa5ecdbFqu1yQdtB7KVqkoyJXigw3Ooq/kOpElk50uU3B6z4rmrvc+r+ML0R6A1nZDkjhdWxZTZdaL6ua1kX/cWeK3hFZpW10WEOdxIeiz0VqVm2XS9XFiu3iaglXLyJZNtqIIkTuT0gY8jtGR+HGDKyS1Xqi2iKdzNbJoJkiq9llUuCyWYXW9WtEnhPDWYcstdB2uD2rxLr8ylKlz15Yt8PTvf8QWaFtCkk5J4u2il4Wuc2K7XKpulixXTTKqzoFtVKiq0XaRuoweovUrP7vxVULk/3tKeXRJzp2ikaYrGz+FqeOnLghS/rInNXq8bGUWYdkzol0fPVuUaRovugn0TxRRZEVzFtd3peC66PLG0X6Q+v5AyWbprT2UiqHSMu8Pkv6b8WK7XKpulixXZpIG6wTabvofzieEqnpb8BK0Q7RJ6ISIhoJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJkAAJGEBgmKRR2oB0mAQJkAAJkIBNCeyWelWyad1YLRIgARIggTwS0LfgvxHpy4768tPTIo1CrC/X6ZuyavqW/DKRRmLVF6E03pDabtGLIj1XX5SKEdFIgARIgARsRuBGqc87HnUqL+u7Re4ehS4XitShqD0qcr9Ju1vW3W/Q95f1WSIaCZAACZCAzQjUl/rsFr0g6iBS2y1yO4q/y/phkTsERqqsa/wktd0iDbugpgHrjjjX+IcETELAHWHQJMVhMUjAsgQ0RlhzUU/RKJHGpvI0DUb3vaiv506PdY2l5DbPdfc+LkmABEiABCxO4HIpv37DWE17D1+KdM5BoxGraaTPPSL3/IMOQWkvRG23aLiuiN0m+tq5xj8kYBIC7FGYpCFYDMsTaCw1GCvS70ycE90raiP6TvSbqLPoDtHHohIitSdE2hNRqyDSKKBnRJfqdcghGgmQAAmQQEEksFsq7Z7LKIj1Z51NTqCwycvH4pEACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZAACZBAHgj8P+cxTHjIxtqqAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2 誤差・精度(25点)\n",
    "\n",
    "最初に紛れ込んだ丸め誤差が次第に拡大されて，最後には真の値と全く違う値となる\n",
    "「不安定な」例として「Numeriacl Recipes in C」で紹介されている漸化式の誤差を検討する．\n",
    "\n",
    "次のような，いわゆる黄金比\n",
    "$$\n",
    "\\phi = \\frac{\\sqrt{5} - 1}{2} = 0.61803399\n",
    "$$\n",
    "の累乗計算を考える．素直に累乗(power)で求めた場合と，\n",
    "\\begin{align}\n",
    "\\phi^{n+1} &= \\phi^{n-1} - \\phi^n \\\\\n",
    "\\phi^0 &= 1.0 \\\\\n",
    "\\phi^1 &= 0.61803398\n",
    "\\end{align}\n",
    "という漸化式(recurrence formula)で求めた場合とで，数値を%20.15fで出力して，\n",
    "それらの誤差をn=30程度までで議論せよ．\n",
    "\n",
    "以下は，それぞれのリスト(phi_power, phi_recur)を片対数でプロットした結果である．\n",
    "\n",
    "「Numeriacl Recipes in C, C言語による数値計算のレシピ」, W.H.Press他(技術評論社, 1993), p.44.\n",
    "\n",
    "``` python\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "phi1 = 0.61803399\n",
    "phi_recur = [1]\n",
    "phi_recur.append(phi1)\n",
    "phi_power = [1]\n",
    "phi_power.append(phi1)\n",
    "step = [0, 1]\n",
    "\n",
    "for i in range(2,31):\n",
    "    # ...\n",
    "    #ここを考える\n",
    "    # ...\n",
    "\n",
    "plt.plot(step, phi_power, color = 'r', label=\"power\")\n",
    "plt.plot(step, phi_recur, color = 'b', label=\"recur\")\n",
    "plt.legend()\n",
    "plt.xlabel('step')\n",
    "plt.ylabel('phi^n')\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 2:    0.3819660128,    0.3819660100\n",
      " 3:    0.2360679789,    0.2360679800\n",
      " 4:    0.1458980349,    0.1458980300\n",
      " 5:    0.0901699447,    0.0901699500\n",
      " 6:    0.0557280907,    0.0557280800\n",
      " 7:    0.0344418542,    0.0344418700\n",
      " 8:    0.0212862366,    0.0212862100\n",
      " 9:    0.0131556177,    0.0131556600\n",
      "10:    0.0081306189,    0.0081305500\n",
      "11:    0.0050249989,    0.0050251100\n",
      "12:    0.0031056201,    0.0031054400\n",
      "13:    0.0019193788,    0.0019196700\n",
      "14:    0.0011862413,    0.0011857700\n",
      "15:    0.0007331375,    0.0007339000\n",
      "16:    0.0004531039,    0.0004518700\n",
      "17:    0.0002800336,    0.0002820300\n",
      "18:    0.0001730703,    0.0001698400\n",
      "19:    0.0001069633,    0.0001121900\n",
      "20:    0.0000661070,    0.0000576500\n",
      "21:    0.0000408564,    0.0000545400\n",
      "22:    0.0000252506,    0.0000031100\n",
      "23:    0.0000156057,    0.0000514300\n",
      "24:    0.0000096449,   -0.0000483200\n",
      "25:    0.0000059609,    0.0000997500\n",
      "26:    0.0000036840,   -0.0001480700\n",
      "27:    0.0000022768,    0.0002478200\n",
      "28:    0.0000014072,   -0.0003958900\n",
      "29:    0.0000008697,    0.0006437100\n",
      "30:    0.0000005375,   -0.0010396000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAydklEQVR4nO3dd3gVVfrA8e+bBAhNQgcNkkgACQkEEpq0ACrgz7YKCqsgFhCVVXR3FburIpZdG2JBQOwsthULikBCb6GHUEREBemCEJGa9/fHXDTGkITk3sy9k/fzPPMkM3fmzHsYuC/nnJkzoqoYY4wxJxPmdgDGGGOCmyUKY4wxBbJEYYwxpkCWKIwxxhTIEoUxxpgCRbgdQCDUqlVLY2JiinXsL7/8QuXKlf0bkEusLsHHK/UAq0swKmk9li5dultVa+fd7qlEISIXARfFxcWRkZFRrDLS09NJTU31a1xusboEH6/UA6wuwaik9RCR7/Lb7qmuJ1X9RFWHVKtWze1QjDHGMzyVKIwxxvifJQpjjDEF8tQYhTHGFOTo0aNs2bKFQ4cO/WF7tWrVWLt2rUtR+U9R6xEZGUl0dDTlypUrUrmWKIwxZcaWLVuoWrUqMTExiMhv2w8cOEDVqlVdjMw/ilIPVWXPnj1s2bKF2NjYIpUb9F1PIlJZRF4XkVdF5Cq34zHGhK5Dhw5Rs2bNPySJskZEqFmz5p9aVQVxJVGIyAQR2SkimXm29xKR9SKyUURG+DZfBryvqoOBi0s9WGOMp5TlJHHCqf4ZuNWimAj0yr1BRMKBMUBvIB7oLyLxQDTwg2+344EM6suRGaQ/vIucYzmBPI0xxoQUV8YoVHW2iMTk2dwW2KiqmwBEZBJwCbAFJ1msoIDEJiJDgCEAdevWJT09/ZTjmvDSESZv7ctXURkMf2ALtdtGnXIZwSQ7O7tYfw7ByCt18Uo9IDTrUq1aNQ4cOPCn7cePH893e6g5lXocOnSo6NdPVV1ZgBggM9d6H2BcrvUBwAtAZeA14CXgqqKUnZycrMWRczxHR/aerFGyVyM5qE/0TtOjvx4tVlnBIC0tze0Q/MYrdfFKPVRDsy5ZWVn5bt+/f38pRxIYJ6vH0aN//h7L788CyNB8vlODfjBbVX9R1WtV9SZVfbugfUXkIhEZ+/PPPxfrXBImnHNnbbKWHaZ3/ZXcNTWVdjW+ZuXk9cUqzxhj8tq8eTNnn302V111Fc2aNaNPnz4cPHiQGTNm0KpVKxITE7nuuus4fPgwS5Ys4bLLLgPg448/pmLFihw5coRDhw5x1llnAfDNN9/Qq1cvkpOT6dmzJ+vWrQNg0KBBDB06lHbt2nHnnXeWKOZguj12K9Ag13q0b1upq59Ulw+21OGDfy7glmfiSLkyihGj07nvsw5UOK2CGyEZY/xt+HBYsQKAisePQ3h4yctMSoJnny10t/Xr1zN+/Hg6duzIddddx9NPP80rr7zCjBkzaNKkCQMHDuSll15i2LBhrPDFOGfOHBISEliyZAnHjh2jXbt2AAwZMoSXX36Zxo0bM3PmTG6++WZmzpwJOLcDz58/n/AS1i2YWhRLgMYiEisi5YF+wJRTKUD9ONeThAl9/tOBtV9HcFWjRTw6N5Wk2luY/8rqEpdtjCnbGjRoQMeOHQG4+uqrmTFjBrGxsTRp0gSAa665htmzZxMREUGjRo1Yu3Ytixcv5o477mD27NnMmTOHzp07k52dzfz58+nbty9JSUkMHz6cbdu2/Xaevn37ljhJgEstChF5F0gFaonIFuBBVR0vIsOAL4FwYIKqrjnFcn+bPdZfajSqzsSNneg/MoMhD9aj09BY/vbSLEZ+kUyVelX8dh5jTCnL9T//X0v5gbu8t6dGRUWxZ8+efPft0qULU6dOpVy5cpx77rkMGjSI48eP89RTT5GTk0NUVNRvrY68D9z5a+p0V1oUqtpfVeurajlVjVbV8b7tn6tqE1VtpKoji1FuwGaP7XlvCpnfV+OWxDmMXtmZhOh9TBu11O/nMcZ43/fff8+CBQsAeOedd0hJSWHz5s1s3LgRgDfffJOuXbsC0LlzZ5599lk6dOhA7dq12bNnD+vXrychIYHTTjuN2NhY3nvvPcC5OWnlypV+jzeYup5KrKSD2YWpenpVRq/qyuwxmUSGH6HnPclc23gOe7/dF5DzGWO8qWnTpowZM4ZmzZqxd+9ebr/9dl577TX69u1LYmIiYWFhDB06FIB27dqxY8cOunTpAkCLFi1ITEz8rVXy9ttvM378eFq2bEnbtm35+OOP/R5vMA1ml5iqfgJ8kpKSMjiQ5+l0cwtW/PUQD1+QzpMLOvFF3B7G/H0hlz3ZPpCnNcZ4REREBG+99dYftvXo0YPly5f/ad+KFSty+PDh39bHjh37h89jY2P54osvgD92PU2cONFv8XqqRVGaIqMieWx+Kkve2Ui9Cnu5/Kn29IlewPZVO90OzRhj/MpTiSLQXU/5adX/bBbvbsRj56fz6dZWxCeV4/XBc9EcLbUYjDGhIyYmhszMzMJ3DCKeShSBHMwuSLlK5bj7y1RWfPYj8VV+YNC4TvSus5Tv5m0p1TiMMSYQPJUo3Hb2BWcx+6cERveZxdw9Z9O8UxQv9J1lkwwaY0KapxKFG11PeYVFhDHsva6smbuPTjXX8bf3u9KlRibrp25yLSZjjCkJTyUKt7qe8tOwYzRTdyYz8Ya5ZGU3oOUFpzOqZzpHDx51OzRjjDklnkoUwUbChGte7cTaVce46Izl3DMtlXa1NrL83XVuh2aMMUVmiaIU1E2ozXtbOvDBPxey7XAN2vw1jnvOSefQvqK/itAY4z2qSk6Of8cwjx075tfywGOJIhjGKApy2ZPtydpYgYGNFzBqQSpJdX9k7our3A7LGFOKNm/eTNOmTRk4cCAJCQk88sgjtGnThhYtWvDggw/+tt8bb7xBixYtaNmyJQMGDACcqcPff//93/apUsWZby49PZ3OnTtz5ZVXEh8f7/eY7cnsUlY9NooJGzrT//GlDLm/Dp1vacEtL89i1BetqXp66U1KZkxZl2uWcY4fr1ias4zz9ddf8/rrr7N//37ef/99Fi9ejKpy8cUXM3v2bGrWrMmjjz7K/PnzqVWrFj/99FOhZS5btoyFCxeSmJhY4nrk5akWRSg5b0Qyq3+ozq0tZ/Hi6s4knPkzX47McDssY0wpaNiwIe3bt2fatGlMmzaNVq1a0bp1a9atW8fXX3/NzJkz6du3L7Vq1QKgRo0ahZbZtm1bYmJiAhKvp1oUoaZKvSo8t6IrV76ymutvrUSv+1IYOGEuz0xrTo1G1d0OzxhPy/0//wMHfi3VacZPTP+tqtx9993ceOONf/h89OjR+R4XERHx25hGTk4OR44c+VOZgWAtiiBwzo2JLN9xBvd2TOedTe1o1vgY7/99gdthGWMCrGfPnkyYMIHs7GwAtm7dys6dO+nevTvvvffeb++oONH1FBMTw9KlzusNpkyZwtGjpXO7vacSRbAPZhckMiqSR+emsmTSJqIjd9P36Q5cfsZCtq3Y4XZoxpgAOf/88/nrX/9Khw4dSExMpE+fPhw4cIDmzZtz77330rVrV1q2bMkdd9wBwODBg5k1axYtW7ZkwYIFAW1F/IGqem5JTk7W4kpLSyv2sf5y9Nej+nivNK3Arxole3XCtbM153jOKZcTDHXxF6/UxSv1UA3NumRlZeW7ff/+/aUcSWCcSj3y+7MAMjSf71RPtSi8IiIygrumprLqi20kVv2O617rTM/ay9g81yYZNMaUPksUQaxJz1jS9yTyYr/ZLPipCQmdo3j+8lkcP3Lc7dCMMWWIJYogFxYRxk3vdmHN/P10qb2W2z7sSueaWaz99Bu3QzMmJDk9LGXbqf4ZWKIIEWd2OIPPtqfwxo1zWf/LGSRdFM3I82ySQWNORWRkJHv27CnTyUJV2bNnD5GRkUU+xp6jCCESJgx4uRM9/7aLW3sv5b7pqUyuuZ4Jr+aQfHUzt8MzJuhFR0ezZcsWdu3a9Yfthw4dOqUvzmBV1HpERkYSHR1d5HKDPlGIyFnAvUA1Ve3jdjzBoE7z2kz6vjb971nETU80pO2A2vxjdDoPTW1HxRoV3Q7PmKBVrlw5YmNj/7Q9PT2dVq1auRCRfwWqHgHtehKRCSKyU0Qy82zvJSLrRWSjiIwoqAxV3aSq1wcyzlB1yWPtyNpUkeuazuPJxam0rLedWc+tcDssY4zHBHqMYiLQK/cGEQkHxgC9gXigv4jEi0iiiHyaZ6kT4PhCXlTDary6rgvTn1zGcQ0jdXgSNzWfzf4t+90OzRjjERLoQR0RiQE+VdUE33oH4CFV7elbvxtAVUcVUs77BXU9icgQYAhA3bp1kydNmlSseLOzs3+bujfUHN57hA9GHGb8hguoH7adO65Io/WNRe+HDGahfF1y80o9wOoSjEpaj27dui1V1ZQ/fZDfU3j+XIAYIDPXeh9gXK71AcALBRxfE3gZ+Aa4u5BzXQSMjYuLK/LTiXmF4tOmeS0ct1qbV9igoHpVzFzdtW632yGVmBeui6p36qFqdQlGJa0HofpktqruUdWhqtpIC2l1aBC9M9tN7a5PYNnuhtzW8n9M3tyGZs2USbfOR3PK7i2BxpjicyNRbAUa5FqP9m0rsVCeFNDfylcpz6XPRrH0/c3EVtpB/9HncOkZi9masc3t0IwxIcaNRLEEaCwisSJSHugHTPFHwdai+LPEy5uw4Kez+feF6Xy1PZH4NpV4deBsa10YY4os0LfHvgssAJqKyBYRuV5VjwHDgC+BtcBkVV3jp/NZiyIf4eXD+fsnqayavovWUZsY8mYXetRcwTczv3M7NGNMCAhoolDV/qpaX1XLqWq0qo73bf9cVZv4xh1G+vF81qIoQFyPhszY1ZJXrprN0n1nkdijNk9fkm6TDBpjChT0g9mnwloUhQuLCGPIW11Ys/ggPepk8vcpqZxTYy2ZH33tdmjGmCDlqURhLYqii25Tnynb2vDOsPlsOliP1pc15F/d0jmSfaTwg40xZYqnEoW1KE6NhAn9R59D1hro23AJD6WnklzrO5a8nuV2aMaYIOKpRGEtiuKp3awWb2/uyCf3L2bv0Sq0H9SUf6Skc3D3QbdDM8YEAU8lClMyFz7cljXfVmJws3n8Z2kqLervIv3ZFW6HZYxxmacShXU9lVy1M6vxclYX0p5ZAUC325O4sdlsfv7e/kyNKas8lSis68l/UocnsWpbbf6Rks64dR2Jjz3IJ/cvdjssY4wLPJUojH9VqlWJp5aksnDiemqWP8DFj7blrzHz2LV2t9uhGWNKkSUKU6g218STsSuGf3VL5/3v2tCsufDOLfNsGhBjyghPJQobowic8lXK88DMVJb/73viKm3jqhc7cnH9JWxZYpMMGuN1nkoUNkYReM0viWPeT814+pJ0ZuxMIL5tZV65ajY5x3LcDs0YEyCeShSmdISXD+f2/6WSmbabNtU3MvSdLvSovZKNM2ySQWO8yBKFKbazUs9k+u5WvDpwDsv2nUXiuXX494XpHDt0zO3QjDF+ZInClIiECTe83pmsJQc5v94q/vlZKufUWs/qDza4HZoxxk8sURi/OCOlPv/b2pZJt85n88E6tO4Ty4Nd0zm8/7DboRljSshTicLuenKXhAlXPncOa9eH0y92EQ/PTqV1nR9YOC7T7dCMMSXgqURhdz0Fh5qNa/Dmpk589tAS9h+rxDmD47kjOZ1fdv7idmjGmGLwVKIwweWCB9uwZnMVhjafyzPLUkk8Yw8znlrmdljGmFNkicIE1GnRp/FiZhdmPb+SCDnOuXe2ZvDZs9n3nXUPGhMqLFGYUtHlby1Zub0ed7ZNZ8L6jsSf9Ssf37PI7bCMMUVgicKUmoo1KvLEolQWvb6e2uV/5tJR7bjyzPnsyNzldmjGmAKERKIQkUtF5FUR+a+InO92PKZkUgbGk7HnLB7pkc7/fkgmvkU4b91kkwwaE6wCnihEZIKI7BSRzDzbe4nIehHZKCIjCipDVf+nqoOBocCVgYzXlI5ylcpx3/RUVnyyhaaVtzLg5Y78X70Mvl+w1e3QjDF5lEaLYiLQK/cGEQkHxgC9gXigv4jEi0iiiHyaZ6mT69D7fMcZj2h2YSPm7InnuctmMWtXPM3POY0X+82ySQaNCSKiGvjmvojEAJ+qaoJvvQPwkKr29K3fDaCqo05yvACPA1+p6vST7DMEGAJQt27d5EmTJhUr1uzsbKpUqVKsY4NNqNVl7+oDPH9fbdL3t6dD5aXc/sAP1G4bBYReXU7GK/UAq0swKmk9unXrtlRVU/70gaoGfAFigMxc632AcbnWBwAvFHD8rcBS4GVgaGHnS05O1uJKS0sr9rHBJhTrknM8R1+7fo5GyV6twK/6eK80Pfrr0ZCsS368Ug9Vq0swKmk9gAzN5zs1JAazVfV5VU1W1aGq+vLJ9rMpPEKfhAmDxnVi7Yoj/N/pKxjxRSrtanzN1pk/uR2aMWWWW4liK9Ag13q0b5sxANRrUYcPtrbn/X8sYOvhmlzzyMXc2zGdQ/sOuR2aMWWOW4liCdBYRGJFpDzQD5hS0kLV5nrynMuf6kDWhnJcfnoaj81PpVXdrcx/ZbXbYRlTppTG7bHvAguApiKyRUSuV9VjwDDgS2AtMFlV1/jhXNb15EE1GlXnprfL8cWjGfx6vAKdhjbn1pazyN6e7XZoxpQJAU8UqtpfVeurajlVjVbV8b7tn6tqE1VtpKoj/XQua1F4WM97U8jcEsWwFnN4YVVnEqL3MW3UUrfDMsbzQmIwu6isReF9VepV4fmVXZnzYiaR4UfoeU8y1zaew95v97kdmjGe5alEYS2KsqPjTS1YseN07jknnTc3diA+7jAf3rnQ7bCM8SRPJQpTtkRGRTJyXioZk76hfoWfuPyp9vSJXsD2VTvdDs0YT/FUorCup7Ip6cqmLNodx6ie6Xy6tRXxSeWYeMNcm2TQGD/xVKKwrqeyq1ylcoz4IpWVn/9I86rfc+34TvSqs5TNc7e4HZoxIc9TicKYpr3PYtaeRF7oO4v5e5qS0DmK0X1skkFjSsJTicK6ngxAWEQYt0zuSubcn+lUcx23ftCVLjUyWff5JrdDMyYkeSpRWNeTya1hx2im7kzm9cFzycpuQMv/O4PHzk/n6MGjbodmTEjxVKIwJi8JEwaO7cTaVce4JHoZ936VStta37Ds7bVuh2ZMyLBEYcqEugm1mfxDBz68cyHbD1en7dWNubtDOr/utUkGjSmMpxKFjVGYwvzlifZkbazANY0X8PjCVJLqbmPui6vcDsuYoOapRGFjFKYoqsdGMX5DZ756YhlHNILOt7RgWItZHPjxgNuhGROUPJUojDkV597ZmtU/VGd4Ujovru5M8zP3M/WRDLfDMiboWKIwZVqVelV4Znkq815ZQ5XwQ1zwQAoDG81jz9f2Rj1jTrBEYQzQYUgiy3dFc3+nNN7d1JZmTXOYfPsCmwbEGDyWKGww25REhdMq8PCcbiydvIkzK+7kymc78JczFvPjsu1uh2aMqzyVKGww2/hDi75NWbinCU9ekM6X21sQn1yR8dfOsdaFKbOKnChE5BwR+auIDDyxBDIwY9wUERnBPz9LZdWX22lZ7VtumNiZ82ov59vZP7gdmjGlrkiJQkTeBP4NdALa+JaUAMZlTFBofH4sabtb8FL/2Sz+KY6ErjV47rJZHD9y3O3QjCk1EUXcLwWIV1Vre5syJywijKHvdOH/bvuRGy/ayvCPuvLfmqsZ93ZF4i+Oczs8YwKuqF1PmUC9QAZiTLBr0O50Ptuewls3zWPDL6fT6pIGPHquTTJovK+oiaIWkCUiX4rIlBNLIAMzJhhJmHDVix3JWp3DX85cyv0zUkmpuYmlb2a5HZoxAVPUrqeHAhlEQUSkGXAbTrKaoaovuRWLMSfUaV6bSd/Vpv+9i7np8TNpN7AWf38hnYemtqNijYpuh2eMXxWpRaGqs/JbCjtORCaIyE4RycyzvZeIrBeRjSIyopBzr1XVocAVQMeixGtMablkZFuyNlXk2qYLeHJxKi3rbWfWcyvcDssYvyrqXU+Diln+RKBXnrLCgTFAbyAe6C8i8SKSKCKf5lnq+I65GPgM+LyYcRgTMFENq/Hqus5Mf2o5xzWM1OFJDG0+m59/2O92aMb4hRR2I5OIPAA0UdWri3UCkRjgU1VN8K13AB5S1Z6+9bsBVHVUEcr6TFX/7ySfDQGGANStWzd50qRJxQmX7OxsqlSpUqxjg43VpfQd2XuYD0YcYdyGC6gXtoMRV88j8drav30eKvUoCqtL8ClpPbp167ZUVf/86IOqnnQBxgLvAmEF7VdIGTFAZq71PsC4XOsDgBcKOD4VeB54BbilkHNdBIyNi4vT4kpLSyv2scHG6uKeReNXa0KF9Qqq/RvO051rd6tq6NWjIFaX4FPSegAZms93a2FdT/2BR1Q1p9gpqoRUNV1Vb1XVG1V1TCH72hQeJii0vS6Bpbsa8q/UNN7/LoVm8fD2MJtk0ISmwhLFRcBkEWnkx3NuBRrkWo/2bSsxmxTQBJPyVSvwQFo3ln+4mcaVtnL1mA48dnkFtizZ5nZoxpySAhOFqqYD/YC3/HjOJUBjEYkVkfK+8v3yTIa1KEwwav6XJszdE88zF89k/r6WNG9bibEDbJJBEzoKvetJVTOBy4tTuIi8CywAmorIFhG5XlWPAcOAL4G1wGRVXVOc8vM5n7UoTFAKrxDB8I+78+bT00mJ+oYb3+pMj5or2JT2nduhmSCXlgaNGsHBg+7FUNTnKH4sTuGq2l9V66tqOVWNVtXxvu2fq2oTVW2kqiOLU/ZJzmctChPUarSqxvTdSYy9ejYZ+xqR2L0Wz12aZpMMmpO64w7YtAnWr3cvhgIThYjc6fs5WkSez7uUTohFZy0KEwokPIzBb3Yha/EvpNbJYvjH3ehccw3rpmxwOzRj8lVYi2Kt72cGsDSfJahYi8KEkug29fl0Wwpv3ryA9b9Ek3TJmYw6dwbHDh5xOzRj/qDAuZ5U9RPfz9dLJ5ySEZGLgIvi4mzqZxMaJEy4ekwHzrtlD8N6reSeGT14v2YWE149TsurE90Ozxig6FN4NBGRsSIyTURmnlgCHdypshaFCVV142vy3vfteH9EBluP1CJlwNnc2246h35ycQTTGJ+iTjP+HrAcuA/4Z67FGONHl49KIeubSK5usoTHFp9Li3o7SH9mudthmTKuqInimKq+pKqLVXXpiSWgkRWDDWYbL6gRcxqvrT+H6f9ZSQ5hdLujFTecPYefNtskg8Ydhd31VENEagCfiMgtIlL/xDbf9qBiXU/GS3rc0ZJVP9bmrrZpTFzfgWaNDvPfO5diLyQ2pa2wFsVSnDuergH+AczzrZ9YjDEBVKlWJR5f1I2Mt9ZzZoUd9HsqmYuil/H98j1uh2bKkMKm8IhV1bNw3hvxArASWAGMBpoHPDpjDABJVzVn4Z4mPN1rGmk/NiW+dQWeG5DB8WPWvDCBV9QxiteBZjjTfY/GSRxBd8usjVEYLwuvWJ7bp57Pmi+30rnaaoa/lcI5tTewZrpNMmgCq6iJIkFVb1DVNN8yGEgIZGDFYWMUpiyIOb8Jn+9uyzsDprJpXw2Sz6vOf/rM5/hR194GYDyuqIlimYi0P7EiIu2wMQpjXCMR4fR/ozdrFv1C79oZ/OODc+hWezWbZnzrdmhl3ujRcP/9bkfhX0VNFMnAfBHZLCKbcWaEbSMiq0VkVcCiM8YUqE7bGD7c3pHXr5vFyp9jaHFubcZe/iV69JjboZVZt94Kjz7qdhT+VeAUHrn0CmgUxphikzBh4PiupA7dwXUXbOfGD3vyv1rzGPdBdU4/N97t8IwHFHWa8e8KWgIdpDGmcGe2qcu07S0Yfd0y0ve3IuG8eky6bDIcPux2aCbEFbXrKSTYXU+mrAsLF4aNb82KhYdpWusn+n90BVfWmcmeL21I0RSfpxKF3fVkjKNJu+rM2RbHyGvW89H+HiT0OoPJ//c6mv2L26GZEOSpRGGM+V1EBNwzsSlL5h6hXq1jXPn5NXSqtZaFo5e4HZoJMZYojPG4lh2rkLG9AePu3MCmY2fS4dY29DtrMZtX2SSDpmgsURhTBoSHw/VPNOHrbVW5/5wZTPk2gbNblueuyzZgQ3qmMJYojClDqtSuyMPzerDhs430qz6Npz6KI67OfsaM2s/Ro25HZ4JVSCQKEaksIhkicqHbsRjjBdEXtGDijt5k3PwaCUeXMeye02jR8Gc+/URtGnPzJwFNFCIyQUR2ikhmnu29RGS9iGwUkRFFKOouYHJgojSmjCpXjtZjrmdmZl0+bnonOdu2c9HFQr+Lf+HQIbeDM8Ek0C2KieR5qltEwoExQG+cWWj7i0i8iCSKyKd5ljoich6QBewMcKzGlEkS34yL14wi85npjCz3EJM/rUzPxK3s+6lsTDK4ahVkZ4e7HUZQEw1wO1NEYoBPVTXBt94BeEhVe/rW7wZQ1VEnOX4kUBknqfwK/EVV//Q3WESGAEMA6tatmzxp0qRixZudnU2VKlWKdWywsboEn2CvR+T27ay+ZzW3fPsgjSM38/jjq6jasla++wZ7XYqqW7dUmjTZyyuvrPRbeQBpael+KW/w4GQ2bqzK2LEZNG6cXeC+Jb0m3bp1W6qqKX/6QFUDugAxQGau9T7AuFzrA4AXilDOIODCopwzOTlZiystLa3YxwYbq0vwCYl65OToV//4QquwX8/kO82641XVo0f/tFtI1KUIwFmCtbykJKe8ZcsK37ek1wTI0Hy+U0NiMBtAVSeq6qcF7WNTeBjjByKc+1RPZk87zOEKVen09F9Y0PwGWOmf/3Gb0ONGotgKNMi1Hu3bZowJIq3Oq8X8rOrUrF+eHhte5JPWD8J992Ej3WWPG4liCdBYRGJFpDzQD5jij4LV5noyxq/OOgvmraxKQqvyXKofMm7kdmjVCubPdzs0U4oCfXvsuzgvOWoqIltE5HpVPQYMA74E1gKTVXWNn85nXU/G+Fnt2jBzdgTn9wxjMON4ZNsNaMdOxD3/PGQXPLhqvCGgiUJV+6tqfVUtp6rRqjret/1zVW2iqo1UdaQfz2ctCmMCoEoVmDIFBg6EB37+OzcnzKbeR1MgIQGmTXM7PBNgRX3DXUgQkYuAi+Li4twOxRjPKVcOJk6E00+Hxx/vxKe1dhO3cz01e/5AzSbp1LygHbUaVKRmTf6wNGgAFSu6Hb0pCU8lClX9BPgkJSVlsNuxGONFIjBqFDRpAuPHH+OYtiFrQyx7Nih7NpTjeD7HxMbCN984x5aUKqxZ4zRkTOkJmdtjjTHB49pr4dFHM5kzL4ysXXXYsXwbR5PaspcoNva8hUWf7+Hzz+G22+DbbyEzs/Ayi2LqVEhMdJKFKT2eShQ2mG2MS5KSkCWLiXr8bhqlj6ftVY3pvWMitw93Zn6YOdM/p/n8c+fnfnuVRqnyVKKwwWxjXBQRAXfd5TyY17w5XHstDW/sRaOGR/2WKGbM8E855tR4KlEYY4JA06YwaxaMGQPz59N961ukTzvMsSMlm2Rw61ZYt85PMZpT4qlEYV1PxgSJsDC4+WZYs4buzXew/1AFlqcMhrVri12ktSbc46lEYV1PxgSZM8+k2xd3ATBjY0NISoLHHqM4r9ObPt3PsZki81SiMMYEn7r1hIQEmNn2Lrj0Urj3XmjTBpYtK3IZqk6iqFEjcHGak7NEYYwJuO7dYe7iChx+47/w0UewYwe0bQsjRsCvvxZ6/Lp1sG0b9OhRCsGaP/FUorAxCmOCU/fuTj5YtAinVZGVBYMGwRNPON1Rc+YUePyJbidLFO7wVKKwMQpjglPXrs749m8D0tWrw7hx8NVXcOQIdOkCt9wCBw7ke/z06c5MtrGxpRez+Z2nEoUxJjhFRUHr1vk8eHfuuc5j28OHw0svOc9fTJ36h12OHYP0dGdX4w5LFMaYUtGjByxcCL/8kueDypXhmWecd1xUrQoXXOBMU7tnDwAZGc6T2Nbt5B5LFMaYUtG9u9M6mDv3JDu0b+/cCXX//fDuuxAfD++9x4zp+tvxJ6gGPl7zO08lChvMNiZ4dezoTFVe4HQeFSrAww/D0qXO/ORXXMH05zJplXCEWrVKLVSTh6cShQ1mGxO8Kld2Gg1FesK6RQtYuJCDjz7N/N1N6LHhZRg/3poSLvFUojDGBLcePZzepb17i7BzRARz29zOESpwbpPv4YYbnOcuTKmzRGGMKTXduzuNglmzirb/9OlOd1Wn+U86d0WtX+988N//wvH8XpNkAsEShTGm1LRr57wWtajTjk+fDuecA5WrhsHQofDqWOeD55+DTp2cB/dMwFmiMMaUmvLloXPnoo1T7N4NK1bkeX6idh3n5wMPwtdfQ6tW8MgjzkN7JmAsURhjSlWPHk5DYPv2gvdLS3O6qfJ9fqJnT6eQyy6DBx6AlBRYsiQg8ZoQSBQikioic0TkZRFJdTseY0zJnHgeIi2t4P1mzHCev2vT5vdtIrl2qFPHed7i44+dh/Pat4c774SDB/0ec1kX0EQhIhNEZKeIZObZ3ktE1ovIRhEp7DYGBbKBSGBLoGI1xpSOVq2gWrXCxymmT4fUVOcNqwW6+GJYswauvx6eegpatiz6aLkpkkC3KCYCvXJvEJFwYAzQG4gH+otIvIgkisineZY6wBxV7Q3cBfwrwPEaYwIsPNxJAAUlis2b4ZtvTmF+p6goGDvWaYbk5DgnuOkmZ+4PU2KF5eoSUdXZIhKTZ3NbYKOqbgIQkUnAJao6CriwgOL2AhVO9qGIDAGGANStW5f09PRixZydnV3sY4ON1SX4eKUeULK6NGhwBh9/3JhJkxZSr96hP33+2Wf1gLM57bTFpKf/3pW0cmV1oCXLli3jyJF8kkBYGGFjxhA7YQLRY8dy+IMP2HD77fzUoUMB0aQC+PG6+Le87OxkoCoZGRn8/HN2IfsG6O+XqgZ0AWKAzFzrfYBxudYHAC8UcPxlwCvAf4HUQs51ETA2Li5OiystLa3YxwYbq0vw8Uo9VEtWl8xMVVAdPz7/z/v1U61fXzUn54/bp01zjps7twgnWbhQtXlz54CrrlLdtSvf3Zwh81OLvyD+Li8pySlv2bLC9y3p3y8gQ/P5bg36wWxV/VBVb1TVK1U1vZB9bQoPY0JAfLwzFp1f91NOjtOD1KNHnsHrU9WunfMY+IMPwuTJ0KwZTJpk04AUgxuJYivQINd6tG9bidmkgMaEBhHn7qeZM//8vZ2ZCbt2+Wla8fLl4aGHnEkGY2Ohf3/nDXtb/fKVU2a4kSiWAI1FJFZEygP9gCn+KNhaFMaEju7dnfdgn5iV44SAvPY0MREWLIB//9t5q158PLz6qrUuiijQt8e+CywAmorIFhG5XlWPAcOAL4G1wGRVXeOn81mLwpgQcSIR5H1Ke8YMaNrUmWXcr8LD4e9/h1WrnNftDRlib0MqooAmClXtr6r1VbWcqkar6njf9s9VtYmqNlLVkX48n7UojAkRsbHQsOEfxymOHHEegQjo93dcnJONXnnF6ZI6wSYZPKmgH8w+FdaiMCZ0nBinSEtzBrABFi1yXpUa8Pdjh4U5LYo1uTozOnSA1asDfOLQ5KlEYS0KY0JL9+7OuylWrnTWZ8xwvsNTU0spgOjo33//9ltITnYGvw8fLqUAQoOnEoW1KIwJLSfmfToxTjF9uvNdXb26C8GsXQt9+8K//uUEsWiRC0EEJ08lCmtRGBNaTj8dzj7bGac4cMD5bi6o26lEz1UUplYtePtt+OQT2LfP6Yq64w6nL6yM81SiMMaEnu7dYfZsp1Vx7FgQ3Ih04YXOFOY33gjPPOO8v7uob1ryKE8lCut6Mib0dO/u/Kd91CiIjISOHd2OCDjtNOfVq+npzqBJjx4weLDT0iiDPJUorOvJmNCTmup0KS1e7LzdNDKy8GNK7Tm5rl2dkfZ//hMmTIDmzWGKX54PDimeShTGmNBTsyYkJTm/u97tlJ9KleDJJ50BlJo14ZJLoF8/2LnT7chKjSUKY4zrTtz9FPDnJ/I4pZZJSgpkZDjv6P7oI2cakLffLhPTgHgqUdgYhTGhadgw5/u3dWu3IylE+fJw332wfDk0bgxXX+0Mfv/wg9uRBZSnEoWNURgTmmJinO/fsFD5RoqPh7lz4dlnnQHv5s2dwe8Tj5h7TKhcFmOMCS7h4XDbbc686O3awc03Q7dubkcVEJYojDGmJGJjYdo0GD/+97lIwHkoxCMsURhjTEmJwHXXOQ/qndC+/R8TRwizRGGMMf5y+um///7DD86dUvffH/KTDHoqUdhdT8aYoJGVBX/9Kzz6KLRq5bxhL0R5KlHYXU/GmKBRsya8/jpMnerMUdKxIwwfHpKTDHoqURhjvC2gs8cGSq9ezp1RN98Mzz0HCQm/vxg8RFiiMMaEnJB7GLpqVXjhBWea3PLl4bzz4PrrQ2aSQUsUxpgyq9QTTufOzp1QI0Y43VLNmjnTgQQ5SxTGGFOaIiOdOdUXL4Z69eCyy+CKK2DHDrcjO6mgTxQiEiYiI0VktIhc43Y8xhjjF61bO8li5Ej4+GOndfHGG0HZrxbQRCEiE0Rkp4hk5tneS0TWi8hGERlRSDGXANHAUWBLoGI1xphSV64c3HOP0x3VrBlccw1ccAF8/73bkf1BoFsUE4FeuTeISDgwBugNxAP9RSReRBJF5NM8Sx2gKTBfVe8AbgpwvMYYU/rOPhvmzIHnn3d+Nm8OY8YEzSSDEYEsXFVni0hMns1tgY2quglARCYBl6jqKODCvGWIyBbgiG/1+MnOJSJDgCEAdevWJT09vVgxZ2dnF/vYYGN1CT5eqQe4U5cVK6KAJJYvX87x4yV/sNb5Hk4F8GNdSlBeYiKR48bR5D//ocawYex75RUO/pwO1CAjI4Off84u8PCAXRNVDegCxACZudb7AONyrQ8AXijg+ErAeGA0cEtRzpmcnKzFlZaWVuxjg43VJfh4pR6q7tRl+nRVUJ092z/lHT/ulAf+KU/VT+Xl5Ki+9ppqVJQmyXIF1WWLjhR6WEmvCZCh+XynBv1gtqoeVNXrVfVvqjqmoH1tCg9jjCeIwKBBsHYtnHaas23gQOeFSS5wI1FsBRrkWo/2bTPGGJNbvXoQe5bz+65d0KYN3HsvHDpUqmG4kSiWAI1FJFZEygP9gCn+KFhtridjjFd98AEMGACPPQZJSTBvXqmdOtC3x74LLACaisgWEbleVY8Bw4AvgbXAZFVd46fzWdeTMcabqlWD116DL790WhSdO8Pf/gYHDgT81AFNFKraX1Xrq2o5VY1W1fG+7Z+rahNVbaSqI/14PmtRGONhJyYFDMJn0krP+ec7kwwOG+bcQpuQ4CSPAAr6wexTYS0KY8ypCNmEU6XK789cVKrkzFA7aBAR+/cH5HSeShTWojDGlCkdOzp3Qt17L7z1Fm0HDYL58/1+Gk8lCmOMKXMiI5236GVkkN2oEcTF+f0UnkoU1vVkjCmzkpJY9dRTUKeO34v2VKKwridjjPE/TyUKY4wx/uepRGFdT8YY43+eShTW9WSMMf7nqURhjDHG/yxRGGOMKZCnEoWNURhjjP95KlHYGIUxZUPITr1RAm7WWdSDf+Iisgv4rpiH1wJ2+zEcN1ldgo9X6gFWl2BU0no0VNXaeTd6MlGUhIhkqGqK23H4g9Ul+HilHmB1CUaBqoenup6MMcb4nyUKY4wxBbJE8Wdj3Q7Aj6wuwccr9QCrSzAKSD1sjMIYY0yBrEVhjDGmQJYojDHGFMgSRS4i0ktE1ovIRhEZ4XY8xSUim0VktYisEJEMt+M5FSIyQUR2ikhmrm01ROQrEfna97O6mzEW1Unq8pCIbPVdmxUicoGbMRaFiDQQkTQRyRKRNSJym297yF2XAuoSitclUkQWi8hKX13+5dseKyKLfN9j/xWR8iU+l41ROEQkHNgAnAdsAZYA/VU1y9XAikFENgMpqhpyDxCJSBcgG3hDVRN8254EflLVx30JvLqq3uVmnEVxkro8BGSr6r/djO1UiEh9oL6qLhORqsBS4FJgECF2XQqoyxWE3nURoLKqZotIOWAucBtwB/Chqk4SkZeBlar6UknOZS2K37UFNqrqJlU9AkwCLnE5pjJHVWcDP+XZfAnwuu/313H+YQe9k9Ql5KjqNlVd5vv9ALAWOIMQvC4F1CXkqCPbt1rOtyjQHXjft90v18USxe/OAH7Itb6FEP0LhPOXZZqILBWRIW4H4wd1VXWb7/ftQF03g/GDYSKyytc1FfTdNbmJSAzQClhEiF+XPHWBELwuIhIuIiuAncBXwDfAPlU95tvFL99jlii8qZOqtgZ6A7f4ukA8QZ2+0lDuL30JaAQkAduA/7gazSkQkSrAB8BwVd2f+7NQuy751CUkr4uqHlfVJCAap1fk7ECcxxLF77YCDXKtR/u2hRxV3er7uRP4COcvUCjb4etbPtHHvNPleIpNVXf4/nHnAK8SItfG1wf+AfC2qn7o2xyS1yW/uoTqdTlBVfcBaUAHIEpEInwf+eV7zBLF75YAjX13DJQH+gFTXI7plIlIZd8gHSJSGTgfyCz4qKA3BbjG9/s1wMcuxlIiJ75Yff5CCFwb36DpeGCtqj6d66OQuy4nq0uIXpfaIhLl+70izo04a3ESRh/fbn65LnbXUy6+W+KeBcKBCao60t2ITp2InIXTigCIAN4JpXqIyLtAKs50yTuAB4H/AZOBM3Gmj79CVYN+kPgkdUnF6d5QYDNwY65+/qAkIp2AOcBqIMe3+R6cvv2Qui4F1KU/oXddWuAMVofj/Kd/sqo+7PsOmATUAJYDV6vq4RKdyxKFMcaYgljXkzHGmAJZojDGGFMgSxTGGGMKZInCGGNMgSxRGGOMKZAlCmMCRESGi0glt+MwpqTs9lhjAiSUZ/E1JjdrURjjB74n4j/zvRsgU0QeBE4H0kQkzbfP+SKyQESWich7vvmGTrw/5Elx3iGyWETi3KyLMXlZojDGP3oBP6pqS9+7J54FfgS6qWo3EakF3Aec65uwMQPnvQEn/KyqicALvmONCRqWKIzxj9XAeSLyhIh0VtWf83zeHogH5vmmhb4GaJjr83dz/ewQ6GCNORURhe9ijCmMqm4QkdbABcCjIjIjzy4CfKWq/U9WxEl+N8Z11qIwxg9E5HTgoKq+BTwFtAYOAFV9uywEOp4Yf/CNaTTJVcSVuX4uKJ2ojSkaa1EY4x+JwFMikgMcBW7C6UL6QkR+9I1TDALeFZEKvmPuw3lPO0B1EVkFHMaZydSYoGG3xxrjMruN1gQ763oyxhhTIGtRGGOMKZC1KIwxxhTIEoUxxpgCWaIwxhhTIEsUxhhjCmSJwhhjTIH+H2tlBwSwuRF0AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "phi1 = 0.61803399\n",
    "phi_recur = [1]\n",
    "phi_recur.append(phi1)\n",
    "phi_power = [1]\n",
    "phi_power.append(phi1)\n",
    "step = [0, 1]\n",
    "\n",
    "for i in range(2,31):\n",
    "    recur = phi_recur[-2]-phi_recur[-1]\n",
    "    power = phi1**(i)\n",
    "    phi_recur.append(recur)\n",
    "    phi_power.append(power)\n",
    "    step.append(i)\n",
    "    print(\"%2d: %15.10f, %15.10f\" % (i, power, recur))\n",
    "\n",
    "plt.plot(step, phi_power, color = 'r', label=\"power\")\n",
    "plt.plot(step, phi_recur, color = 'b', label=\"recur\")\n",
    "plt.legend()\n",
    "plt.xlabel('step')\n",
    "plt.ylabel('phi^n')\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 数値積分，解，収束性(25点)\n",
    "\n",
    "数値積分に際して，関数値の計算にコストがかかる場合，限られた点での値から高精度で計算するNewton-Cotesの公式を使う．\n",
    "4区間5点の閉区間(closed)公式(Boole's formula)は次のとおりである．\n",
    "\\begin{align}\n",
    "h = & \\frac{b-a}{4} \\\\\n",
    "\\int_a^b f(x) dx = & \\frac{2}{45}h\\left\\{7f(a)+32f(a+h)+12f(a+2h)\\right. \\\\\n",
    "& \\left. +32f(a+3h)+7f(a+4h)\\right\\}\n",
    "\\end{align}\n",
    "\n",
    "これに従って求めた\n",
    "\\begin{align}\n",
    "\\int_0^1 \\frac{4}{1+x^2} dx & \\,(= \\pi)\n",
    "\\end{align}\n",
    "の値は以下の通りである．\n",
    "\n",
    "数値積分の中点則を用いて求めた場合，同じ程度の精度を達成するには何点が必要となるか．\n",
    "``` python\n",
    "import numpy as np\n",
    "def func(x):\n",
    "    return 4.0/(1+x**2)\n",
    "\n",
    "n = 4    \n",
    "h = 1.0/n\n",
    "\n",
    "boole = 2/45*h*(7*func(0)+32*func(h)+12*func(2*h)+32*func(3*h)+7*func(4*h))\n",
    "print(\"%26s : %15.10f\" % (\"from Boole's formula\",boole) )\n",
    "print(\"%26s : %15.10f\" % (\"actual Error\", boole-np.pi))\n",
    "print(\"%26s : %15.10f\" % (\"Estimated Error f^(6)(0.5)\", h**7*8/945*1311.8)) # 1311.8=$f^6(0.5)$\n",
    "```\n",
    "\n",
    "```\n",
    "      from Boole's formula :    3.1421176471\n",
    "              actual Error :    0.0005249935\n",
    "Estimated Error f^(6)(0.5) :    0.0006778067\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      from Boole's formula :    3.1421176471\n",
      "              actual Error :    0.0005249935\n",
      "Estimated Error f^(6)(0.5) :    0.0006778067\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "def func(x):\n",
    "    return 4.0/(1+x**2)\n",
    "\n",
    "n = 4    \n",
    "h = 1.0/n\n",
    "\n",
    "boole = 2/45*h*(7*func(0)+32*func(h)+12*func(2*h)+32*func(3*h)+7*func(4*h))\n",
    "boole = 2/45*h*(7*func(0)+32*func(h)+12*func(2*h)+32*func(3*h)+7*func(4*h))\n",
    "print(\"%26s : %15.10f\" % (\"from Boole's formula\",boole) )\n",
    "print(\"%26s : %15.10f\" % (\"actual Error\", boole-np.pi))\n",
    "print(\"%26s : %15.10f\" % (\"Estimated Error f^(6)(0.5)\", h**7*8/945*1311.8)) # 1311.8=$f^6(0.5)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1:   3.2000000000    0.0584073464\n",
      "  2:   3.1623529412    0.0207602876\n",
      "  3:   3.1508492099    0.0092565563\n",
      "  4:   3.1468005184    0.0052078648\n",
      "  5:   3.1449258640    0.0033332104\n",
      "  6:   3.1439074272    0.0023147736\n",
      "  7:   3.1432933175    0.0017006639\n",
      "  8:   3.1428947296    0.0013020760\n",
      "  9:   3.1426214566    0.0010288030\n",
      " 10:   3.1424259850    0.0008333314\n",
      " 11:   3.1422813577    0.0006887041\n",
      " 12:   3.1421713566    0.0005787031\n",
      " 13:   3.1420857498    0.0004930962\n",
      " 14:   3.1420178234    0.0004251698\n",
      " 15:   3.1419630238    0.0003703702\n",
      " 16:   3.1419181743    0.0003255207\n",
      " 17:   3.1418810041    0.0002883506\n",
      " 18:   3.1418498552    0.0002572016\n",
      " 19:   3.1418234938    0.0002308402\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "def mid(N):\n",
    "    x0, xn = 0.0, 1.0\n",
    "\n",
    "    h = (xn-x0)/N\n",
    "    S = 0.0\n",
    "    for i in range(0, N):\n",
    "        xi = x0 + (i+0.5)*h\n",
    "        dS = h * func(xi)\n",
    "        S = S + dS\n",
    "    return S\n",
    "\n",
    "x, y = [], []\n",
    "for i in range(1,20):\n",
    "    x.append(i)\n",
    "    error = abs(mid(i)-np.pi)\n",
    "    y.append(error)\n",
    "    print(\"%3d:%15.10f %15.10f\" % (i, mid(i), error))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ちなみに問１で求めた補間多項式を積分した値は，3.142183333となりBoole公式で求めた値と大体一致する．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 微分方程式 (25点)\n",
    "\n",
    "\n",
    "Verlet法による小惑星軌道のシミュレーションを次のような条件で行った．\n",
    "```python\n",
    "def force(pos):\n",
    "    x=pos[0]\n",
    "    y=pos[1]\n",
    "    L=(x*x+y*y)**(3/2)\n",
    "    return [-x/L,-y/L]\n",
    "\n",
    "def Verlet(r0,rh):\n",
    "    f=force(r0)\n",
    "    x=2*r0[0]-rh[0]+h**2/m*f[0]\n",
    "    y=2*r0[1]-rh[1]+h**2/m*f[1]\n",
    "    return [x,y]\n",
    "\n",
    "h=0.1\n",
    "dx, dy=0.0, h\n",
    "m=0.2\n",
    "xx=[3.0, 3.0-dx]\n",
    "yy =[0.0,-dy]\n",
    "```\n",
    "\n",
    "1. この小惑星軌道をplotせよ．\n",
    "1. この小惑星の公転周期が規格化した単位でいかほどになるか答えよ．\n",
    "1. 異なる初期条件を使って軌跡を表示し，その飛翔体の振る舞いを解説せよ．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<function matplotlib.pyplot.show(close=None, block=None)>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ+klEQVR4nO3dfYxc913v8c/HjhPY61uK7VWbJvFugFApLfQho9zmgq4KFOT4SjUtFCVaIIWiVaRGBME/qSwVUSkSFRK3LURQq42ayquWCij1pdvrm0BRQPQh4ypJnbgRbmQ3DoFsbGhjcktI8r1/zAwZr8+ZOWfnPM0575d0NDNnfnvO2bOz3/md7+/hOCIEAGi/bXUfAACgGgR8AOgIAj4AdAQBHwA6goAPAB1xSd0HMMmePXtieXm57sMAgLlx7NixZyJiMem9Rgf85eVl9fv9ug8DAOaG7dNp75HSAYCOIOADQEfMHPBtX2X7i7Yftf2I7dsTytj2R2yftP2w7TfPul8AQD5F5PBfkPRbEfE12/9V0jHb90bEo2NlbpR0zXD5b5L+aPgIAKjIzDX8iHgqIr42fP6spBOSrthU7ICkT8bAlyW90vbls+4bAJBdoTl828uS3iTpK5veukLSE2Ovz+jiL4XRNlZt9233NzY2ijw8IL+1NWl5Wdq2bfC4tlZMWaAGhQV82zsl/Zmk34iI72x1OxFxKCJ6EdFbXEzsSgrMJmtgXluTVlel06eliMHj6mpy+bxl+WJAHSJi5kXSDklHJf1myvsflXTz2OvHJF0+bbvXXXddAFMdPhyxtBRhDx4PH55cdmEhYhCWB8vCQvLPLC1dWG60LC1tvWye/ef93YCIkNSPtFid9kbWRZIlfVLShyaU+Z+SvjAs+xZJX82ybQJ+R5UVwCPyBXE7uay99bJ59p/1d+NLAWPKDvg/LikkPSzpweGyX9Ktkm6Nl78U7pL0TUlfl9TLsm0CfgeVGcAj8gXxMmr4Re8/7/lC65Ua8MtcCPgdVGYAz7v9PME0a9mirzDy/j5cCbQeAR/1yxpsygzgo+MoK4eepWzRbQhZzxdXAp1BwEe9ymoozbvt8Z+ps6abdf9Zfres5ytPozJXAXONgI96lZVGGf+Ztgapab9b1vOV5UqAq4BWIOCjeHmCbN40TZsDeBmynK8sX7p5r67QSAR8FKvsnjQoXpa/WZ72AL6QG2tSwGd6ZOR38KD03HMXrnvuucH6JHfeKS0sXLhuYWGwHtVYWZEOHZKWliR78Hjo0GD9yN69yT87vj7PiGI0DgEf+X3rW/nWZwk2KN/KinTqlPTSS4PHzec/yxdzli97po5oLAI+BvL8k2apCW42Ldigflm+mKd92XMF0GgepHyaqdfrBfe0rcDon3S85rawkF4Lz1se7bG8PAjimy0tDb7Ip72P0tk+FhG9pPeo4SN/Tp4UTXdNS/tkSfeR8qkNNXwM/vGSPgf2IAUDjFtbG1QGvvWtQRrvzjtf/rKfVsPn6rB01PC7qOycPLprUnvMtCuAvFeTKBQBv43yNpzRbRJFmZbuy9vDC4Ui4LcROXnUadIVwLSrSfL7pSKH30bk5NFUk3L4Evn9ApDD7xpy8miqSVeT5PdLV0jAt3237adtH095/622v237weHy/iL2ixTk5NFkaSmfLIO6SPfMpKga/ick7ZtS5m8j4o3D5QMF7bc78nzYycljHk26MmUEbyEKCfgRcb+kc0VsCwm28mFnKgPMm0lXpqR7ClFlDv8G2w/Z/oLt16UVsr1qu2+7v7GxUeHhNRgfdnTBpCtTunMWoqqA/zVJSxHxBkl/IOkv0gpGxKGI6EVEb3FxsaLDazg+7OiKtCtTunMWopKAHxHfiYjzw+frknbY3lPFvluBXjfouknpHvL7mVUS8G2/2raHz68f7vdsFftuBXrdoOvozlmIorplfkrSlyS91vYZ2++xfavtW4dFfl7ScdsPSfqIpJuiySO+qpL1MpReN8DWu3PiPzHSti7MGggUY9IMnaMePkkze7YUI22biMtQoBhpKc/9+8ntb0LArwuXoUAx0lKe6+tUqjYhpVMXbgUHlKujkwiS0mkiet4A5Zo2VUMH++0T8OtCzxugXOT2L0LArxPz3QDlIbd/EQJ+GTp6uQg0TlKlqsMdJgj4RWOYN9Bsabn9XbtaX1Ej4BeN/vVAsyXl9nfskJ59tvUVNQJ+0Tp8uQjMhaTc/iteIT3//IXlWlhRI+AXjZktgebbnNs/l3L/ppZV1Aj4RaN/PTB/OpLXJ+AXjf71wPzpSF6fqRUAQBoE8vGZNc+fl84m3Laj4dOfMLVCkehjD7RTB/L6BPw86GMPdEcL8/pF3fHqbttP2z6e8r5tf8T2SdsP235zEfutHH3sge5oYV6/qBr+JyTtm/D+jZKuGS6rkv6ooP1Wiz72QHe0sL9+IQE/Iu6XlJLwkiQdkPTJGPiypFfavryIfVeKPvZAt7Qsr19VDv8KSU+MvT4zXHcR26u2+7b7GxsblRxcZvSxB7ptzit9jWu0jYhDEdGLiN7i4mLdh3Mh+tgD3ZaW1z9/fi4acS+paD9PSrpq7PWVw3XzZ2WFAA901eh/f9Rff9euQSPuqL/+qBF3vGyDVFXDPyLpl4e9dd4i6dsR8VRF+waA4ozn9XfunKtG3EJq+LY/JemtkvbYPiPptyXtkKSI+GNJ65L2Szop6TlJv1LEfgGgVnPWc6+oXjo3R8TlEbEjIq6MiI9HxB8Pg72GvXPeGxE/GBE/EhHNnC+BUbQA8pizRtzGNdrWhlG0APJKasS1B/GjgZVGAv4Io2gB5DXec08aBPvRhJQNrDQyW+bItm0v/6HG2YPGGQCYZHl5EOQ3q3h2TWbLzGLOcnEAGmYOGnAJ+COMogUwi7TK4bZtjekIQsAfYRQtgFkkVRol6cUXG9MRhBw+ABRl/K5Z27YNgv1mJef0yeEDQBXGR+GmdfaoMadPwAeAMjSwI0g3Az4jagGUrYGDsqqaLbM5RiNqR4OsGj67HYA5NT6z5unTyYOyxstVoHuNtg0ZHAGgQyqMOzTajpuDwREAWqYhcad7Ab+BDSkAWq4hcad7AZ8RtQCqljYo6/z5ShtvuxfwGVELoGqjuLN794Xrz56tdPRtIY22tvdJ+rCk7ZI+FhG/u+n9d0v6Pb18H9s/jIiPTdsuI20BtEoFjbeTGm1n7pZpe7ukuyT9tKQzkh6wfSQiHt1U9E8i4rZZ9wcAc6vmxtsiUjrXSzoZEY9HxPOSPi3pQAHbBYB2qbnxtoiAf4WkJ8Zenxmu2+znbD9s+09tX1XAfgFgvtTceFtVo+3/lrQcET8q6V5J96QVtL1qu2+7v7GxMfuemUYBQFPU3HhbRMB/UtJ4jf1Kvdw4K0mKiLMR8e/Dlx+TdF3axiLiUET0IqK3uLg425FxY3IATbOyIu3cefH6Cu6hXUTAf0DSNbavtn2ppJskHRkvYPvysZdvl3SigP1Ox43JATRRTY23M/fSiYgXbN8m6agG3TLvjohHbH9AUj8ijkj6ddtvl/SCpHOS3j3rfjNpyHBmALjA3r3J3TNLbrwtJIcfEesR8cMR8YMRcedw3fuHwV4R8b6IeF1EvCEifiIivlHEfqdqyHBmALhATY237R5pyzQKAJqopsbbdgd8plEA0FQ1NN52bz58AGiKbdtevinKODv9nrhTMB8+ADRRxe2MBHwAqEvafW/37y9ldwR8AKjLyop0yy2DID8SId1zTykNtwR8AKjT+vrFefySGm7bGfCZPwfAvKhwgGj7Aj7z5wCYJxU23LYv4DN/DoB5UmHDbfsCPvPnAJgnFTbcti/gM38OgHlTUcNt+wI+8+cAmDcVZSbaF/CZPwfAvKkoM9G+gC8NgvupU4O5KE6dItgDaLaKGm7bGfABYJ5U1HBbSMC3vc/2Y7ZP2r4j4f3LbP/J8P2v2F4uYr8A0BoVNNzOHPBtb5d0l6QbJV0r6Wbb124q9h5J/xIRPyTpf0n64Kz7BYBWqaDhtoga/vWSTkbE4xHxvKRPSzqwqcwBSfcMn/+ppJ+yx69dCsS0CgDm0a5d+dZvQREB/wpJT4y9PjNcl1gmIl6Q9G1Jm+7tVQCmVQCAVI1rtLW9artvu7+xsZHvh5lWAcC8Ons23/otKCLgPynpqrHXVw7XJZaxfYmk75OU+FtExKGI6EVEb3FxMd+RMK0CgHm1fXu+9VtQRMB/QNI1tq+2famkmyQd2VTmiKRbhs9/XtJfRxk302VaBQDz6sUX863fgpkD/jAnf5uko5JOSPpMRDxi+wO23z4s9nFJu22flPSbki7qulkIplUAMK92pzRrpq3fgkuK2EhErEta37Tu/WPPvyvpXUXsa6LRiNqDBwdpnL17B8GekbYAIJeRWSlKr9eLfr9f92EAQPkm9VTPEadtH4uIXtJ7jeulAwCdNCeNtgCAWc1Doy0AoAAVNNoS8AGgCb773dJ30c6Az3w6AObJ2pr0b/+W/N65c4XtppBumY0ymk9nNMXCaD4die6ZAJpp0vQvBQ4cbV8Nn/l0AMybSdO/FDhwtH0Bn/l0AMybtCmQd+4sNDPRvoDPfDoA5k1ag+1llxW6m/YFfObTATBPKmqwldoY8FdWpEOHpKWlwVDlpaXBaxpsATRRRQ22Uht76UiD4E6ABzAPKmqwldpYwweAeVJRg61EwAeAelXUYCsR8AGgPhU22EozBnzbu2zfa/sfho/fn1LuRdsPDpfNtz8sF9MsAGiq229Pf6+EruSz1vDvkPRXEXGNpL9S+q0L/19EvHG4vD2lTPFG0yycPj24gcBomgWCPoC6ra1JZ8+mv19CV/KZ7nhl+zFJb42Ip2xfLulvIuK1CeXOR8TOvNuf+Y5Xy8uDIL/Z0pJ06tTWtwsAs0qLT9JgSuRnntnSZsu849WrIuKp4fN/kvSqlHLfY7tv+8u2f3bGfWbHNAsAmiot2EvShz9cyi6n9sO3fZ+kVye8dcFogYgI22mXC0sR8aTtH5D017a/HhHfTNnfqqRVSdo7aw5r797kk8o0CwDqtLY2GBialGHZvbu0cURTa/gR8baIeH3C8jlJ/zxM5Wj4+HTKNp4cPj4u6W8kvWnC/g5FRC8ieouLi1v4lcYwzQKAJrr99uRgb5dWu5dmT+kckXTL8Pktkj63uYDt77d92fD5Hkk/JunRGfebDdMsAGiaSY21EaXGp1kbbXdL+oykvZJOS/qFiDhnuyfp1oj4Ndv/XdJHJb2kwRfMhyLi41m2P3OjLQA0zZ496QG/gA4lkxptZ5pLJyLOSvqphPV9Sb82fP73kn5klv0AQCvU0BVzHCNtAaAqk2bGLLGxdoSADwBVqaEr5rhuBnymWwBQtVFXzCQV1O6lts6HP8louoXRjc5H0y1I9N4BUI61NemWW2rpinnBrmbppVO2UnrpMN0CgCptrmQmKTAOlzm1wvxhugUAVTp4cHKwX1qq7FC6F/DTplVgugUAZZjUUFvxyP/uBXymWwBQlUkNtdu3Vz7yv3sBn+kWAFRhWkPtPfdUHne612gLAGWruKF2HI22AFClBjXUjiPgA0DRGtRQO46AvxmjcAFs1draYDbMNDU01I7r3kjbSRiFC2CrpuXta2qoveAQaLQdwyhcAFs16abkIxXEWxpts2IULoCtWFubHuxraqgdR8AfxyhcAHmNUjmTNGRw50wB3/a7bD9i+6XhbQ3Tyu2z/Zjtk7bvmGWfpWIULoA8RoOrJnXB3L27MYM7Z63hH5f0Tkn3pxWwvV3SXZJulHStpJttXzvjfsvBKFwAWY1q9i++mF7m8GHpmWcaE0NmvaftCUly2lwRA9dLOhkRjw/LflrSAUmPzrLv0qysNOaPA6DBsgyualgsqSKHf4WkJ8ZenxmuS2R71Xbfdn9jY6P0gwOA3KY10jY0FTw14Nu+z/bxhOVAGQcUEYciohcRvcXFxTJ2MRsGZgHdNRpY9Yu/mF6m5sFVk0xN6UTE22bcx5OSrhp7feVw3fxhYBbQXVkmRFtYaGywl6pJ6Twg6RrbV9u+VNJNko5UsN/iJeXsnntusB5Au03L2UuNDvbS7N0y32H7jKQbJH3e9tHh+tfYXpekiHhB0m2Sjko6IekzEfHIbIddEwZmAd2UdWBVg4O9NHsvnc9K+mzC+n+UtH/s9bqk9Vn21Qh79yb/0RmYBbTXHA2smoaRtnkwMAvonmmpnAYNrJqGgJ8HA7OA7pmUsm3YwKppCPh5rawMZs586aXB45z8oQFktLnr9a5dyeXmIGe/GfPhA8BIUtfrSy+VduyQ/uM/Xi43p6lcavhlYYAWMH+S8vXPPy+94hWtSOVSwy8DA7SA+ZSWrz93bpCrn3PU8MvAAC1gPmTN17ek6zU1/DIwQAtovpbn65NQwy8Dd84Cmq/l+fok1PDLcOedF0+y1KJaAtAKLc/XJ6GGXwYGaAHN07F8fRICflmyDtCi+yZQvlG+/vRpKWLw+Oyzg3z9uJZfiRPw65T0IVxdJegDRetgvj6JI6LuY0jV6/Wi3+/XfRjlWV5Onn1zaWlwVQCgGNu2DSpVm9mDq/AWsX0sInpJ71HDrxPdN4HiJaVJ6TkniYBfLz6EQLHS0qT79zO1uWa/49W7bD9i+yXbiZcQw3KnbH/d9oO2W5yjyYn59YFipY1yX1+n55xm74d/XNI7JX00Q9mfiIh2dm7dqtGH7eDBQRpn795BsO/YhxAozKQ06cpK5/+3ZqrhR8SJiHisqIPpJLpvAltDrj63qnL4Ien/2j5me+LNIW2v2u7b7m9sbFR0eA1H903gQuTqt2RqwLd9n+3jCcuBHPv58Yh4s6QbJb3X9v9IKxgRhyKiFxG9xcXFHLtoMWbfBC5Ern5Lpgb8iHhbRLw+Yflc1p1ExJPDx6clfVbS9Vs/5A6i+ya6Ki2VOS1Xz21IE5U+eZrt/yJpW0Q8O3z+M5I+UPZ+W2Xv3uQBWuQl0WaTbiTE/8SWzNot8x22z0i6QdLnbR8drn+N7fVhsVdJ+jvbD0n6qqTPR8T/mWW/nbOV7ps08mLeTUpl0qV5ayKisct1110XGDp8OGJpKcIePB4+PLnswkLEoDlrsCwsTP4ZoGnsCz/Do8UevJ/nf6JDJPUjJaYyl04bMUcP5snaWvJYFD7HWzJpLh1ugNJGNPJiXkzK03MjocIxl04bMfgE82JSnp4bCRWOgN9GeRu0aOBFmSZ9vqZdjdLFslAE/DbKUzNiFC/KNO3zxdVopQj4bZW1ZsQoXhQhrRY/7fNF98pK0WjbdTTwYlaTGl6zpGwkZoytCDX8rtvKJTU5f4ybVIvP8vkiT18ZAn7XbaWBl5x/t0z7gp9Uiydl0yxpI7KasDDStiJ5RiwuLSWPflxaquZYUa0so7anfSYYEVspTRhpW3tQn7QQ8Bto2nD3zfhnb75Jf6MsX/BM5dEokwI+KR3kkyfnT/qn+ab9jbI06jNAam4Q8JFPnpzsVrp80iBcrGnnc9rfKOsXPA2v8yGt6t+EhZROQ2VN02wl/UNqoDhZzmeWGSn5m8wVkcNHLfI28OYt3+X2gSy/e5bzmTVH39XzPIcI+KhH3tphniuCvNuep6A17Viz/u5Zzic1+NYpLeBL+j1J35D0sAb3qn1lSrl9kh6TdFLSHVm3T8BvgbK6fOYpW+aXQ9Fli+gGmbfcPH0ZYqoyA/7PSLpk+PyDkj6YUGa7pG9K+gFJl0p6SNK1WbZPwO+YPIE5z9VAWV8OZZTNcqxZf3dq751USUpH0jskrSWsv0HS0bHX75P0vizbJOB3UNbaZp4gXtaXQxllsxxr3i8wau+dMingF9kt81clfSFh/RWSnhh7fWa4LpHtVdt92/2NjY0CDw9zIWv3vjzdQ/OMHcgzmVwZZbMca57fne6SGDM14Nu+z/bxhOXAWJmDkl6QNHOn6Yg4FBG9iOgtLi7Oujm0VZ7BPmV9OZRRNsuxMtAJW5VW9c+6SHq3pC9JWkh5n5QO6pc1tVF3Dj/PsQIJVGKj7T5Jj0panFDmEkmPS7paLzfavi7L9gn4qEWdvXSAGU0K+B68vzW2T0q6TNLZ4aovR8Sttl8j6WMRsX9Ybr+kD2nQY+fuiMg0N2qv14t+v7/l4wOArrF9LCJ6Se/NdMeriPihlPX/KGn/2Ot1Seuz7AsAMBsmTwOAjiDgA0BHEPABoCMI+ADQETP10imb7Q1Jp+s+jqE9kp6p+yAainOTjnOTjnOTbpZzsxQRiaNWGx3wm8R2P62rU9dxbtJxbtJxbtKVdW5I6QBARxDwAaAjCPjZHar7ABqMc5OOc5OOc5OulHNDDh8AOoIaPgB0BAEfADqCgJ+D7XfZfsT2S7bpTibJ9j7bj9k+afuOuo+nKWzfbftp28frPpamsX2V7S/afnT4/3R73cfUFLa/x/ZXbT80PDe/U+T2Cfj5HJf0Tkn3130gTWB7u6S7JN0o6VpJN9u+tt6jaoxPaHC/CFzsBUm/FRHXSnqLpPfyuflP/y7pJyPiDZLeKGmf7bcUtXECfg4RcSIiHqv7OBrkekknI+LxiHhe0qclHZjyM50QEfdLOlf3cTRRRDwVEV8bPn9W0glNuM91lwzvYXJ++HLHcCmsZw0BH7PIdYN6YDPby5LeJOkrNR9KY9jebvtBSU9LujciCjs3M90ApY1s3yfp1QlvHYyIz1V9PEBb2d4p6c8k/UZEfKfu42mKiHhR0httv1LSZ22/PiIKaQsi4G8SEW+r+xjmyJOSrhp7feVwHTCR7R0aBPu1iPjzuo+niSLiX21/UYO2oEICPikdzOIBSdfYvtr2pZJuknSk5mNCw9m2pI9LOhERv1/38TSJ7cVhzV62v1fST0v6RlHbJ+DnYPsdts9IukHS520frfuY6hQRL0i6TdJRDRrePhMRj9R7VM1g+1OSviTptbbP2H5P3cfUID8m6Zck/aTtB4fL/mk/1BGXS/qi7Yc1qFDdGxF/WdTGmVoBADqCGj4AdAQBHwA6goAPAB1BwAeAjiDgA0BHEPABoCMI+ADQEf8fwGGw2l1j4MQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "def force(pos):\n",
    "    x=pos[0]\n",
    "    y=pos[1]\n",
    "    L=(x*x+y*y)**(3/2)\n",
    "    return [-x/L,-y/L]\n",
    "\n",
    "def Verlet(r0,rh):\n",
    "    f=force(r0)\n",
    "    x=2*r0[0]-rh[0]+h**2/m*f[0]\n",
    "    y=2*r0[1]-rh[1]+h**2/m*f[1]\n",
    "    return [x,y]\n",
    "\n",
    "h=0.1\n",
    "dx, dy=0.0, h\n",
    "m=0.2\n",
    "xx=[3.0, 3.0-dx]\n",
    "yy =[0.0,-dy]\n",
    "\n",
    "max_i = 100\n",
    "for i in range(1,max_i): \n",
    "    x1, y1 = Verlet([xx[-1],yy[-1]],[xx[-2],yy[-2]])\n",
    "    xx.append(x1)\n",
    "    yy.append(y1)\n",
    "\n",
    "plt.plot(xx,yy,'o',color='r')\n",
    "plt.show\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-16 2.958663685489252\n",
      "-15 2.977664855859823\n",
      "-14 2.9910895191270233\n",
      "-13 2.99894878154761\n",
      "-12 3.00125172107917\n",
      "-11 2.9980053601418253\n",
      "-10 2.9892146658513203\n",
      "-9 2.9748825770008547\n",
      "-8 2.955010058153888\n"
     ]
    }
   ],
   "source": [
    "d_period = -12\n",
    "width = 4\n",
    "for i in range(d_period-width, d_period+width+1):\n",
    "    print(i, xx[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.8"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(max_i+d_period)*h"
   ]
  },
  {
   "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
}
