{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "
\n", "代数方程式(fsolve)\n", "
\n", "
\n", "
\n", "file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/fsolve\n", "
\n", "https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python\n", "
\n", "cc by Shigeto R. Nishitani 2017-19\n", "
\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 概要\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "代数方程式の解$f(x)=0$を数値的に求めることを考える．標準的な\n", "> 二分法(bisection method)とニュートン法(Newton's method)\n", "\n", "の考え方と例を説明し，\n", ">収束性(convergency)と安定性(stability)\n", "\n", "について議論する．さらに収束判定条件について言及する．\n", "\n", "\n", "二分法のアイデアは単純．中間値の定理より連続な関数では，関数の符号が変わる二つの変数の間には根が必ず存在する．したがって，この方法は収束性は決して高くはないが，\n", "確実．一方，Newton法は関数の微分を用いて収束性を速めた方法である．しかし，不幸にして収束しない場合や微分に時間がかかる場合があり，初期値や使用対象には注意\n", "を要する．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# pythonの標準関数による解法\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pythonでは代数方程式の解は，solveで求まる．\n", "\n", "$$\n", "x^2-4x+1 = 0\n", "$$\n", "の解を考える．未知の問題では時として異常な振る舞いをする関数を相手にすることがあるので，先ずは関数の概形を見ることを常に心がけるべき．" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdXZ/vHvIkApk9aBKIWCoqIUBY0VkVaJP2utonWuirZKNVgcwAnFEaEgFXGoWBURpBXEAecqjqHqq6igaFGcyosIgr4WBQIWhDy/P55Y1AIJ5+yTdc4+9+e6zgXBnHPuneDNyt5rrxXMDBERSY8GsQOIiEiyVOwiIimjYhcRSRkVu4hIyqjYRURSRsUuIpIyKnYRkZRRsYuIpIyKXUQkZRrGeNOtttrK2rdvn9FzV6xYQbNmzZINFImOJf+k5ThAx5KvsjmWmTNnfmZmW9f2eVGKvX379syYMSOj506bNo2ePXsmGygSHUv+SctxgI4lX2VzLCGED+vyeToVIyKSMip2EZGUUbGLiKSMil1EJGVU7CIiKZNIsYcQNg8h3BdCeCeEMCeE0D2J1xURkU2X1HTHG4CpZnZ0CKEx0DSh1xURkU2U9Yg9hLAZsC9wO4CZrTazL7J93fV55hmYNOlHuXhpEZGcWrkSBgyAjz9ukvP3CtnueRpC6AqMAd4GugAzgf5mtuI7n1cBVACUlpaWTZ48eZPf65Zbtufee9tw110v06rVqqxyx1ZeXg5AZWVl5CTJqKqqonnz5rFjZC0txwE6lnzzyCPbcu21HRkx4kW6dVud0WuUl5fPNLM9a/1EM8vqAewJrAG61Xx8AzB0Y88pKyuzTMybZ9agQbVddFFGT88rgPmXPx0qKytjR0hEWo7DTMeST6qrzX78Y7Pddzd79tnKjF8HmGF16OUkLp4uABaY2cs1H98H7JHA6/6Xdu2gR4/PGDMGvvwyF+8gIpK8Z5+Ft96Cs8+GEHL/flkXu5ktBj4KIXSs+aP/h5+WyYkjj1zAkiUwcWKu3kFEJFk33ABbbw3HHVc/75fUPPazgIkhhDeBrsDwhF73v3TpspQuXfwLleXlARGRnPvnP+HRR+H006FJ7q+bAgkVu5nNMrM9zWw3MzvczD5P4nXXJwTo3x9mz4Zp03L1LiIiyRg9GkpKvNjrS0HeeXr88bDVVj5qFxHJV8uXw7hxcOyx0Lp1/b1vQRZ7kybQty88/LD/mCMiko/Gj4dly/yiaX0qyGIH6NcPGjaEP/0pdhIRkf+2dq2fVejeHbp1q9/3Lthib90afv1r/zFn6dLYaUREvu3RR2HuXDjnnPp/74ItdvAvWFUVjB0bO4mIyLddd53fe3PEEfX/3gVd7HvsAfvu66dj1qyJnUZExL3+Ovz973DWWX7KuL4VdLGDj9rnz4cHHoidRETEXXcdNG8Op54a5/0LvtgPPRS2396/kCIisS1aBJMnQ58+sNlmcTIUfLGXlPgNSy+9BNOnx04jIsVu9Gg/NVzfUxy/qeCLHfxfxs03h1GjYicRkWK2YgXcfLNfMO3QIV6OVBR78+Z+w9L99/v0IhGRGMaPh88/h/POi5sjFcUOfvW5QQO4/vrYSUSkGK1d69f69t4b9tknbpbUFPsPfwgnnOA3LH2esyXIRETW76GH/IzB+efHTpKiYgc491w/x3XrrbGTiEixueYan6F3+OGxk6Ss2Lt0gQMO8BuWVme2paCIyCZ76SV/DBjgM/ViS1Wxg/8YtGgRTJoUO4mIFIuRI+EHP4BTTomdxKWu2A88EHbbzb/Q1dWx04hI2r33Hjz4oK8427x57DQudcUeAlxwAbz9Njz+eOw0IpJ2o0ZB48Y+My9fpK7YwZfzbdsWrr46dhIRSbPFi2HCBDj5ZCgtjZ1mnVQWe6NGPkPmuee0zICI5M6NN/pEjdg3JH1XKosdfFW1H/zAz7WLiCRt+XL485/hyCNhxx1jp/m21BZ78+Z+MeOBB/zihohIksaOhS++8Gt6+Sa1xQ5+MeN739OoXUSStXo1XHst7Ldf/e9nWhepLvbSUl/5ccIEWLgwdhoRSYuJE2HBAhg0KHaS9Ut1sYPfsFRdrY04RCQZa9fCH/8Iu+/u983ko9QX+3bbwXHHwS23wJIlsdOISKF76CF491246CK/byYfpb7YAS680BcHu+mm2ElEpJCZwYgRsMMOcNRRsdNsWFEU+667Qq9ecMMNXvAiIpl49ll49VUYODA/FvvakKIodvAfm/71L7jttthJRKRQDR8O224Lv/lN7CQbVzTF3qOHT00aORJWrYqdRkQKzfTpPmI//3yfRp3PEiv2EEJJCOH1EMKjSb1m0i65BD7+2Kc/iohsimHDYMstfX/lfJfkiL0/MCfB10vcAQfAT37iU5XWrImdRkQKxaxZ8OijvpFGs2ax09QukWIPIbQBDgHGJvF6uRKCj9rnzoXJk2OnEZFCMXw4tGwJZ54ZO0ndJDVivx4YCOT91haHHuqzZIYP10YcIlK7d96B++6DM86AzTePnaZuGmb7AiGEXsCnZjYzhNBzI59XAVQAlJaWMm3atIzer6qqKuPnfu3ww1sxdGgnhg6dzX77fZbVa2Ur22PJF0l8X/JBWo4DdCxJueqqnWnceGv22ms606Z9lfXr1cuxmFlWD+AqYAEwD1gMrATu3NhzysrKLFOVlZUZP/dra9aY7bSTWZcuZtXVWb9cRgDzL386JPF9yQdpOQ4zHUsSPvjArKTEbMCA5F4zm2MBZlgdejnrUzFmNsjM2phZe+A44FkzOzHb182lkhK49FJ44w145JHYaUQkX111FTRs6DckFZKimcf+XccfDx06wJAhfpuwiMg3zZvnU6MrKvympEKSaLGb2TQz65Xka+ZKw4Zw8cUwc6Y2vRaR/zZiBDRoUHijdSjiETvASSdBu3YatYvIt330EYwbB7/7HbRpEzvNpivqYm/UyEftL78MTz4ZO42I5IsRI/zXiy6KmyNTRV3sACefDG3bwuDBGrWLiI/Wx471bvjRj2KnyUzRF3vjxn436vTp8MQTsdOISGxXXeWDvEsuiZ0kc0Vf7ACnnOLn2q+4QqN2kWI2f76P1vv08U4oVCp21o3aX3lFM2REitnw4f7rxRfHzZEtFXuNk0+G9u01ahcpVh9+6DNhTj21cM+tf03FXqNRI7jsMpgxA/72t9hpRKS+DRvmK8AW+mgdVOzfctJJfjfq5Zdr5UeRYvLPf8L48X6XaSHOW/8uFfs3NGrkp2Jefx3uvz92GhGpL0OGrLuvJQ1U7N9xwgmwyy4+al+7NnYaEcm1OXPgzjt9E41CWxNmQ1Ts31FS4v96z5kDkybFTiMiuXbFFdC0aWGuCbMhKvb1OPJI6NrV70b9Kvt19UUkT82aBffeC+ecA1ttFTtNclTs69GgAQwd6nujjh8fO42I5Mrll/t2d+eeGztJslTsG3DIIbD33n5a5ssvY6cRkaS9+KJvtHPBBYWzl2ldqdg3IARfM2LhQrjppthpRCRJZjBoEJSWQv/+sdMkT8W+ET17woEHesEvXRo7jYgk5Ykn4Lnn/KbEZs1ip0meir0Ww4fDkiVwzTWxk4hIEqqrfbS+3XZw2mmx0+SGir0WZWVwzDFw3XXwySex04hItu6912fDDBniCwCmkYq9DoYOhX//G/7wh9hJRCQbX30Fl14Ku+7qG9qnlYq9Djp29PWZb73V15QQkcJ0223wwQd+irWkJHaa3FGx19HgwdCwYWHvqiJSzJYvhyuvhH339enMaaZir6PWrf0mhrvvhldfjZ1GRDbVqFHw6adw9dU+nTnNVOybYOBAv+34wgu1GYdIIVm82Ge2HX00dOsWO03uqdg3QcuWPu+1shKmTo2dRkTq6sorYdWqdVvfpZ2KfROdfjpsv72P3rWsr0j+e+cdv2haUQE77hg7Tf1QsW+ixo39TtTZs+GOO2KnEZHaXHihL8t7xRWxk9QfFXsGjjnGFwi79FKoqoqdRkQ2ZNo0ePhhv9O0VavYaeqPij0DIfgV9sWLYeTI2GlEZH2qq+G886BtWxgwIHaa+qViz9A++/jIfeRIXwFSRPLLxInw2mt+wfT734+dpn6p2LMwYgSsWeMzZUQkf6xc6RtTl5X5PsbFJutiDyG0DSFUhhDeDiG8FUJI4erG67f99nD22X4R9bXXYqcRka+NGgULFvivDYpw+JrEIa8BzjOzTsDewBkhhE4JvG5BuPRS2HJLP4enm5ZE4lu40H+aPuoo2G+/2GniyLrYzWyRmb1W8/vlwBzgh9m+bqHYfHNf9fH552HKlNhpRGTQID9FWswTG4IlOMwMIbQHngM6m9my7/y3CqACoLS0tGzy5MkZvUdVVRXNmzfPLmjC1q6Fioo9WbmyhAkTXqVx4+pan1NeXg5AZWVlruPVi3z8vmQiLccBxXksc+a0oF+/Mk444UNOO+1/6yHZpsvm+1JeXj7TzPas9RPNLJEH0ByYCRxZ2+eWlZVZpiorKzN+bi49+6wZmA0fXrfPB8y//OmQr9+XTZWW4zArvmOprjbr3t2stNRs2bLcZ8pUNt8XYIbVoY8TuawQQmgETAEmmtn9SbxmoSkvh8MPh2HDNP1RJIZJk+Cll3x6Y4sWsdPElcSsmADcDswxs2uzj1S4Ro3yc3sXXhg7iUhxWb4cLrgA9twTTj45dpr4khix9wBOAvYPIcyqeRycwOsWnO23979cEyfCCy/ETiNSPIYNg0WLYPTo4pze+F1JzIp5wcyCme1mZl1rHo8lEa4QDRrktzCfdZZWfxSpD++/D9de6yP1YlhrvS70b1vCmjb1Bf1nzfKlQkUktwYM8CUDRoyInSR/qNhz4JhjoGdP3x/1s89ipxFJr0cegcce8yV5S0tjp8kfKvYcCAFuvBGWLvX1KkQkeStX+pIenTr5qU9ZR8WeI507Q//+MHYsvPxy7DQi6XPVVTBvHtx0EzRqFDtNflGx59DgwbDtttCvny6kiiTp/ffh6quhd28/7SnfpmLPoRYt/Gr9a6/BLbfETiOSDmZw5pnQpIlPVJD/pmLPsWOPhQMO8Aupn3wSO41I4ZsyBZ58EoYOhW22iZ0mP6nYcywEv2niyy/h3HNjpxEpbEuX+gXT3Xf3U5yyfir2etCxo9+4NGmSjzREJDOXXup7Dd96KzRsGDtN/lKx15OLLoIdd/RRxpdfxk4jUnhefdVnwJxxBvzkJ7HT5DcVez1p0sQvoP7zn76uhYjU3dq1gYoKn2X2hz/ETpP/9MNMPdp/fzjpJJ+mBZ2AtyMnEikMU6b8kFmz4N57YbPNYqfJfxqx17NRo6BlS4Cx6MsvUru5c2HcuO3o1cv3MZXaqVnq2dZbw/XXA3QHdFlfZGPM4PTToaTEuPlmn2UmtVOxR9C7N8BU4Crmz48cRiSP/eUv8NRTcNppc2nTJnaawqFij8BHHX0BH40kuJ+4SGp88gmccw706AGHHfZx7DgFRcUezXzgEh5/HO68M3YWkfxz1lmwYoXva6BdkTaNvlxRjaZ7d18FctGi2FlE8seUKT4D5oorYJddYqcpPCr2qKoZN87Xle7XT6dkRMA3p+nXD/bYw/cQlk2nYo9s5519MaMHH4S7746dRiS+s8+Gzz+H8eO1znqmVOx54NxzYa+9fCnSTz+NnUYknocegrvu8jVhdtstdprCpWLPAyUlPjqpqoK+fXVKRorT//0fVFRA166+aJ5kTsWeJzp18jUwHnwQ/vrX2GlE6pcZ/P738MUXPnddp2Cyo2LPI+ecAz/7mU/z+uij2GlE6s9dd/lMmCFDYNddY6cpfCr2PPL1KZm1a6FPH6iujp1IJPcWLvSleLt3h/PPj50mHVTseaZDB9/H8emnfe1pkTSrrvZBzOrVMGGCD24keyr2PNS3Lxx8MAwcCG+9FTuNSO7ceKPvKjZqlG9EI8lQseehEGDcOGjRwhcMW7UqdiKR5M2eDRdeCIce6oMZSY6KPU+Vlnq5v/EGXHZZ7DQiyVq1ygctm20GY8dqOd6kqdjzWK9evvrjNdfAs8/GTiOSnIsvhjff9MkCrVrFTpM+iRR7COGgEMK7IYQPQggXJfGa4q65Bjp2hBNP9Bs4RArd44/Dtdf6TJiDD46dJp2yLvYQQglwE/BLfCPP40MInbJ93TSbOHHif37fvn37b338Xc2aweTJsGQJnHKK7kqVwrZoEfz2t75cwDXXxE6TXkmM2PcCPjCzuWa2GpgM/CqB102liRMnUlFR8Z+PP/zwQyoqKjZa7l26+P8Ef/sb3HBDfaQUSV51tW/mvmKFD1aaNImdKMXMLKsHcDQw9hsfnwSM3thzysrKLFNAET8eNFhlsHseZNFDj019XGRgBn3yIEu8R2VlZTb9N6O2TjYzGlJPQggVQAVAaWkp06ZNq6+3TpE+wCzgHqAMWBY3jkid/RQYiv9APy5ylriqqqpy3391af+NPYDuwBPf+HgQMGhjz8lmxJ7Nv3b5oF27duv9V7xdu3Z1ev4LL5iVlJgddZRZdXVus26KQv++fC0tx2GWP8fyySdmrVub7bij2dKlmb1GvhxLEupjxJ7EOfZXgR1DCNuFEBoDxwEPJ/C6qTRs2DCaNm36rT9r2rQpw4YNq9Pze/SAESN8waQbb8xFQpHkrF3rM7r+9S/f6q5ly9iJikPWxW5ma4AzgSeAOcA9ZqYb4Tegd+/ejBkz5j8ft2vXjjFjxtC7d+86v8Z55/ndeuefD6+8kouUIskYPhyeesoHIV26xE5TPBKZx25mj5nZTmbWwczqNvQsYt8s8Xnz5m1SqYPfpXfHHdC6NRx9tOa3S3564gnfjPrEE+HUU2OnKS6687RAbbEF3H+/b6V3/PGwZk3sRCLrzJsHJ5wAnTvDLbdoyYD6pmIvYHvsATffDM8843tEiuSDL7+EI4/08+v33+832Un9qrfpjpIbp5wC06fDH//oG2IfeWTsRFLMzKBfP3j9dXjkEdhhh9iJipNG7Cnwpz9Bt27wm9/AP/4RO40Us9Gj/frPZZf5InYSh4o9Bb73Pf+Rt2VL+NWvfGqZSH175hnft/eww2Dw4NhpipuKPSVat/ZyX7gQjj1WF1Olfs2d63/vOnaEv/4VGqhZotKXP0X23hvGjPG12889N3YaKRbLlvlPimbw8MO6CSkf6OJpyvz2t76BwbXXws47+4UskVxZs8an286ZA1On+mbsEp+KPYWuvhreew/OPtv/R/vFL2InkrQ6/3x47DGfdnvAAbHTyNd0KiaFSkpg0iT48Y/9vOdbWuBBcuDmm31/gAEDfAtHyR8q9pRq0QIefRSaNoVDDoHFi2MnkjR57DE46yyf0qidkPKPij3F2rb1cv/sMy/35ctjJ5I0mDEDjjnGF/WaNMl/QpT8omJPubIyuOceeOMNPy3z1VexE0khmzvXBwmtWvlWjS1axE4k66NiLwIHH+wLMU2dCn37og2xJSOffQa//KUPDh5/HLbZJnYi2RDNiikSp54KH30EQ4b4aGvEiNiJpJAsX+4DhPnz4cknfSqt5C8VexEZPNiX+f3jH2GrrXyqmkhtVq2CI46A117zu5t/9rPYiaQ2KvYiEoIv0rRkCVxwAWy5pa8OKbIhX29t98wzMGGCrwMj+U/FXmRKSnwtjy++8NMzzZr5RVWR76qu9r8j990Ho0b56qFSGHTxtAg1buw/Uu+zD/Tu7et7iHyTGZx5pi/BO3iw1h4qNCr2ItWsmU9X2313n5P8xBOxE0m+MPPrLzffDAMHwuWXx04km0rFXsRatvQpkLvsAocf7rvJS3Ezg0GDfBG5s87y2VPar7TwqNiL3BZbeKHvtJNfGHvyydiJJBYzuPBCnzV1+ulw/fUq9UKlYhe23tpnPXTs6OWu0zLFx8xnSo0cCWecAX/+szbLKGT61gng89qfecZPyxx2GDz0UOxEUl+qq32J51Gj/PTLjTdqpF7oVOzyH1tu6eXetSscdZQv8CTptmYN9Onj9zece64vw6tSL3wqdvmWLbaAp5/2uwtPPBFuvTV2IsmV1avhuOP8xqMrr/Tld1Xq6aBil//SooWvt/3LX/pFtGHDtHBY2ixb5qs0TpniM2Auv1ylniYqdlmv738fHnjAb2C69FK/WWXt2tipJAmLF0PPnlBZCePHwznnxE4kSdOSArJBjRvDX/4CrVv7bInFi+HOO730pTC9/77vgfvJJ/DII/5TmaSPRuyyUQ0a+ObY11/vI/jyci8FKTx//zvsvbcvwVtZqVJPMxW71En//n4+9s03oVs3mD07diLZFBMmwM9/DqWl8PLLsNdesRNJLmVV7CGEkSGEd0IIb4YQHgghbJ5UMMk/RxwBzz/vsyn22cfXmpH8tnatLxFw8smw777w4ouw/faxU0muZTtifwrobGa7Ae8Bg7KPJPmsrMxHfDvsAIceCn/4g9/gIvln+fKG9Orl67307evb2W2uoVdRyKrYzexJM1tT8+F0oE32kSTftW0LL7wAJ5wAl13mq0OuXKmt6vPJ7Nlw+ullPPOM34twyy3QqFHsVFJfgiU0QTmE8Ahwt5nduYH/XgFUAJSWlpZNnjw5o/epqqqiefPmGefMF+Xl5QBUVlZGTpI5M7j33jbcemsHtt12BVdeOYcOHVbEjpWVNPz9mjq1lOuv34mmTb9iyJC36dx5WexIWUvD9+Vr2RxLeXn5TDPbs9ZPNLONPoCngdnrefzqG59zCfAANf9Q1PYoKyuzTFVWVmb83HwCmH/5C9+0aWZbbvlva9LE7PbbzaqrYyfKXCH//VqxwqxPHzMw69nTbMqU/4kdKTGF/H35rmyOBZhhdejYWk/FmNkBZtZ5PY+HAEIIJwO9gN41byxFZr/9YMyYGfz0p/C73/lSBF98ETtVcXnzTZ/pMn68nx57+mnYYovVsWNJJNnOijkIGAgcZmYrk4kkhWiLLb5i6lQYOhTuvhu6dIHnnoudKv2qq31JgJ/8BP71L984ZcgQ39tWile2s2JGAy2Ap0IIs0IItySQSQpUSYkvP/A//+MX6nr29I0b/v3v2MnSad48OPBAOO88v9noH//wj0WynRWzg5m1NbOuNY/TkwomhatbN5g1y0/LXH21LwP80kuxU6VHdbVvhNG5s089HTPG7wreaqvYySRf6M5TyYnmzeG223w3pi+/hB49/O7VZYU/QSOqOXNg//19l6N99vFpjaedppUZ5dtU7JJTBx7opwh+/3vfmWeXXeCee7QM8KZauRIuvtivXbzxxrp/NNu1i51M8pGKXXKuZUu46SaYPh222QZ+/et1hS8bZ+b/EHbqBFdd5TeFvfsunHqqRumyYSp2qTd77QWvvOIj95kz/dx7375aLXJDXnnFd7L69a9hs81g2jS44w5o1Sp2Msl3KnapVyUlvmnHBx/4xsnjxkGHDnDJJfD557HT5Ye33vI9Z7t186/TbbfBa6/5/QIidaFilyi22MLXeH/rLejVC4YP91UHhw4t3oJ/+2046STYdVd46im44gp47z0/7aJ56bIpVOwS1U47weTJPj3yZz/zvTd/9CMYOBA+/jh2uvrxyiu+JPKPf+xr3p9/Pvzv/8LgwX59QmRTqdglL3TpAg8/7AXfqxeMGgXt2/ueq9Onp28WzerVMGkSdO/up1ymTfOlAObP97n/W24ZO6EUMhW75JUuXeCuu/wURL9+8OijXn577ukza5YsiZ0wO+++Cxdd5NMUe/eGzz7zU1Lz5/tSALrJSJKgYpe81KGDF96CBV7oa9f6Rddtt4Vjj4X77vO53YVg4UKfCbTPPrDzznDNNb62y2OPedH37w8tWsROKWnSMHYAkY1p0cJH7v36+Wma8eN9RH/vvdC0qa+Rcsgh8ItfQOvWsdM6M19tcepUP7304ov+5507w8iRvvrlNtvEzSjppmKXgtG1K9xwg59/f/55H7U/+KBfcATYbTefEvjTn/qjvoq+utpv9X/+ed9ZqrJy3YXfrl19ps/RR/toXaQ+qNil4DRsCOXl/hg92u9gffxxePJJuP12P+0BPirefXcv144dfZ/WHXbwG3wyuWvzq6/go498bvkHH/j0xNdf91v8V6xY95777gsHHZRfP0VIcVGxS0ELwUfqu+3mSwR/9ZWX7Ysv+q+zZvmc8DVr1j2nUSMoLfUS3nxzX7CseXP/80WLOjJhAqxaBVVVsHy5r3O+eLFf6Pzm7JwWLfxib58+sMce/lNChw661V/iU7FLqjRq5EsX7LXXuj9bvRo+/BDef99H2osWeVEvXgxLl/qSBsuX+wXaVat+wPe+B02arCv87bbzC5/bbANt2sCOO/pj221V4pKfVOySeo0bryvj2kybNp2ePXvmPJNILmm6o4hIyqjYRURSRsUuIpIyKnYRkZRRsYuIpIyKXUQkZVTsIiIpo2IXEUkZFbuISMqo2EVEUkbFLiKSMip2EZGUUbGLiKSMil1EJGUSKfYQwnkhBAshaI91EZHIsi72EEJb4EBgfvZxREQkW0mM2K8DBgJW2yeKiEjuZVXsIYRfAQvN7I2E8oiISJaC2cYH2iGEp4Ft1vOfLgEuBg40s6UhhHnAnmb22QZepwKoACgtLS2bPHlyRoGrqqpo3rx5Rs/NNzqW/JOW4wAdS77K5ljKy8tnmtmetX6imWX0AHYFPgXm1TzW4OfZt6ntuWVlZZapysrKjJ+bb3Qs+Sctx2GmY8lX2RwLMMPq0M8Zb2ZtZv8AWn39cW0jdhERqR+axy4ikjIZj9i/y8zaJ/VaIiKSOY3YRURSRsUuIpIyKnYRkZRRsYuIpIyKXUQkZWq98zQnbxrC/wEfZvj0rYC0zJXXseSftBwH6FjyVTbH0s7Mtq7tk6IUezZCCDOsLrfUFgAdS/5Jy3GAjiVf1cex6FSMiEjKqNhFRFKmEIt9TOwACdKx5J+0HAfoWPJVzo+l4M6xi4jIxhXiiF1ERDaiIIs9hHBMCOGtEEJ1CKHgrpSHEA4KIbwbQvgghHBR7DzZCCGMCyF8GkKYHTtLNkIIbUMIlSGEt2v+bvWPnSlTIYQmIYRXQghv1BzLlbEzZSOEUBJCeD2E8GjsLNkIIcwLIfwjhDArhDBOL7s3AAACZUlEQVQjl+9VkMUOzAaOBJ6LHWRThRBKgJuAXwKdgONDCJ3ipsrKHcBBsUMkYA1wnpl1AvYGzijg78sqYH8z6wJ0BQ4KIewdOVM2+gNzYodISLmZddV0x/Uwszlm9m7sHBnaC/jAzOaa2WpgMvCryJkyZmbPAUti58iWmS0ys9dqfr8cL5Ifxk2VmZrNdqpqPmxU8yjIi2khhDbAIcDY2FkKSUEWe4H7IfDRNz5eQIEWSFqFENoDuwMvx02SuZrTF7Pw7SufMrNCPZbrgYFAdewgCTDg6RDCzJo9oHMmsY02kraxTbTN7KH6ziPFIYTQHJgCDDCzZbHzZMrM1gJdQwibAw+EEDqbWUFdBwkh9AI+NbOZIYSesfMk4KdmtjCE0Ap4KoTwTs1PvInL22I3swNiZ8iRhUDbb3zcpubPJLIQQiO81Cea2f2x8yTBzL4IIVTi10EKqtiBHsBhIYSDgSZAyxDCnWZ2YuRcGTGzhTW/fhpCeAA/LZuTYtepmPr3KrBjCGG7EEJj4Djg4ciZil4IIQC3A3PM7NrYebIRQti6ZqROCOH7wM+Bd+Km2nRmNsjM2tRsu3kc8GyhlnoIoVkIocXXvwcOJIf/0BZksYcQjgghLAC6A38LITwRO1Ndmdka4EzgCfwC3T1m9lbcVJkLIdwFvAR0DCEsCCH8LnamDPUATgL2r5mONqtmpFiItgUqQwhv4gOJp8ysoKcKpkAp8EII4Q3gFeBvZjY1V2+mO09FRFKmIEfsIiKyYSp2EZGUUbGLiKSMil1EJGVU7CIiKaNiFxFJGRW7iEjKqNhFRFLm/wPAp9Nxz3rPOgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from sympy import *\n", "\n", "x = symbols('x')\n", "\n", "def func(x):\n", " return x**2-4*x+1\n", "\n", "x = np.linspace(-1, 5, 100) #0から2πまでの範囲を100分割したnumpy配列\n", "y = func(x)\n", "plt.plot(x, y, color = 'b')\n", "\n", "plt.plot(0, 0, \"o\", color = 'k')\n", "# plot([x1, x2], [y1, y2], color='k', linestyle='-', linewidth=2)\n", "plt.hlines(0, -1, 5, color='k', linestyle='-', linewidth=2)\n", "plt.vlines(0, -4, 6, color='k', linestyle='-', linewidth=2)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "もし，解析解が容易に求まるなら，その結果を使うほうがよい．\n", "pythonの解析解を求めるsolveは，sympyから呼び出して，" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-√3 + 2, √3 + 2]\n" ] } ], "source": [ "from sympy import *\n", "\n", "x = symbols('x')\n", "\n", "def func(x):\n", " return x**2-4*x+1\n", "\n", "pprint(solve(func(x), x))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "と即座に求めてくれる．数値解は以下の通り求められる．\n", "コメントを外してみてください．ちょっと注意が必要ということがわかるでしょうか？" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 0.26794919]\n", "[ 2.01500001]\n", "[ 0.26794919 3.73205081]\n", "[ 0.26794919 0.26794919]\n" ] } ], "source": [ "from scipy.optimize import fsolve\n", "def func(x):\n", " return x**2-4*x+1\n", "\n", "pprint(fsolve(func, 0))\n", "# pprint(fsolve(func, 2.0))\n", "pprint(fsolve(func, [0, 5]))\n", "pprint(fsolve(func, [0, 0.8]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 二分法とNewton法の原理\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 二分法(bisection) \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "二分法は領域の端$x_1, x_2$で関数値$f(x_1),f(x_2)$を求め，中間の値を次々に計算して，解を囲い込んでいく方法である．\n", "\n", "|$x_1$ | $x_2$ |$f(x_1)$ | $f(x_2)$ |\n", "|:----|:----|:----|:----|\n", "|0.0 | 0.8 |　　　　　 |　　　　　 |\n", "|　　　　| 　　　　| 　　　　| 　　　　 |\n", "|　　　　| 　　　　| 　　　　| 　　　　 |\n", "|　　　　| 　　　　| 　　　　| 　　　　 |\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGX+/vH3Jw1I6C10QlOkCol0QgKIwKpYcAVRXBQRC1Zsq19/uruuroqggiA27KhYQRApKVQlKCC9hioWUCQUKXl+f2TcRQTMMJOcSeZ+Xddczpl5Zs6dMcyd0805h4iIhJ8IrwOIiIg3VAAiImFKBSAiEqZUACIiYUoFICISplQAIiJhSgUgIhKmVAAiImFKBSAiEqaivA5wKpUrV3YJCQmn9dp9+/YRFxcX3EBBoFz+US7/KJd/imOuxYsX/+icq5Kvwc65kL0lJia605WWlnbary1IyuUf5fKPcvmnOOYCslw+v2O1CkhEJEypAEREwpQKQEQkTKkARETClApARCRMBaUAzOxlM/vezJaf5Hkzs2fMbL2ZLTOz1sGY7wm9+SYkJNCla1dISMibFhGRPwjWEsAEoOcpnu8FNPLdhgBjgzTf33vzTRgyBDZvxpyDzZvzplUCIiJ/EJQCcM5lArtPMaQP8JpvN9WFQHkzqx6Mef/O/ffD/v2/f2z//rzHRUTkd8wF6ZrAZpYATHHONTvBc1OAx5xzc33Ts4B7nHNZJxg7hLylBOLj4xMnTpyY7wxdunbN+8v/OM6MjNmz8/0+BSknJ4fSpUt7HeMPlMs/yuUf5fJPILlSU1MXO+eS8jU4v0eM/dkNSACWn+S5KUCnY6ZnAUl/9p5+Hwlct65z8IfboVq1/XufAlQcjzwsSMrlH+XyT3HMRQgeCbwdqH3MdC3fY8H1yCMQG/u7hw5El+DepH68tiCb3NzgLO2IiBQHhVUAnwADfXsDtQP2OOe+DfpcBgyA8ePJBnIB6tbl0HPj+PGCvjz48QqufuVLvt1zIOizFREpioK1G+jbwALgTDPbZmbXmtlQMxvqGzIV2AisB14AbgzGfE9owADqAZEA2dmUG/w3Jgw6h39d1Iys7J/oMTKTj77e/tuqKBGRsBWU00E75/r/yfMOuCkY8zodZsaV7erSqWFl7nxvKbe9s4TpK3byyMXNqRgX41UsERFPhdWRwAmV43j3+vbc07MxM1d9R4+Rmcxc+Z3XsUREPBFWBQAQGWHckNKAT27uRJUyJRj8WhbD31vKLwcPex1NRKRQhV0B/Oas6mX5+KaO3JzakA++2kavUXOYt/5Hr2OJiBSasC0AgJioCIafdybv39CBEtERDHjxC/7fx8vZf+iI19FERApcWBfAb1rVqcCnwzpzTcd6vLpgM72fnsPizac6s4WISNGnAvApFRPJgxc04e3r2nEk19F33AIenbqKg4ePeh1NRKRAqACO075BJT67LZl+59Th+cyNnP/sXJZu/dnrWCIiQacCOIHSJaJ49JLmvHpNG3IOHuGSsfN5cvoaDh3J9TqaiEjQqABOocsZVZh+ezIXt6rJ6LT1XDh6Lsu37/E6lohIUKgA/kS5UtE8eVlLXhyYxK59h7hozDxGzVzL4aNaGhCRok0FkE/dm8Qz4/Zkzm9RnVEz13HRmHms+vYXr2OJiJw2FYAfysfGMKpfK8Zdmch3vxzkwtFzeXbWOi0NiEiRpAI4DT2bVePz27vQs1l1RsxYyyXPzWfNzr1exxIR8YsK4DRVjIvh2f6tGDugNTt+PsD5z85h9GwtDYhI0aECCFCv5tWZcUcXzmtajSc/X8vFz81j9U5tGxCR0KcCCIKKcTGMvqI1Ywe0Zueeg1zwrLYNiEjoUwEEUa/m1fn89i708m0buGjMPFbu0NKAiIQmFUCQVYyL4Zn+v+0p9CsXjp7LyBlrdRSxiIQcFUAB6dmsGjNuT+aCljV4etY6Lhw9l2+26ShiEQkdKoACVCEuhpGXn82LA5PYve8QFz03j0lrD+kMoyISElQAhaB7k3hm3NGFS1rVZMrGw5z/7Fy+2vKT17FEJMypAApJuVLRPHFZS+5MLMH+X49w6dj5/HPKSg4c0tKAiHhDBVDImleJYvrtyQxoW4eX5m7ivFGZLNiwy+tYIhKGVAAeKFMymn9d1JyJQ9oRYdD/hYX8/cNv2HvwsNfRRCSMBKUAzKynma0xs/Vmdu8Jnk8xsz1mtsR3ezAY8y3q2tWvxLRbk7mucz0mfrmFHiMzSVv9vdexRCRMBFwAZhYJjAF6AU2A/mbW5ARD5zjnzvbd/hHofIuLUjGR3P+XJnxwY0fKlIxi0IRF3Dbxa3bvO+R1NBEp5oKxBNAGWO+c2+icOwRMBPoE4X3Dytm1yzN5WCdu7daIKcu+5dynMpi8dAfOOa+jiUgxZYF+wZhZX6Cnc26wb/oqoK1z7uZjxqQAHwDbgO3AcOfcipO83xBgCEB8fHzixIkT/c6UmpoKQFpamt+vLWg5OTmULl36lGO27s3l5eW/smlPLmdXiWRg0xgqlizYzTX5yeUF5fKPcvmnOOZKTU1d7JxLytdg51xAN6Av8OIx01cBo48bUxYo7bvfG1iXn/dOTEx0pwNweT9a6ElLS8vXuCNHc90LmRvcmQ9Mdc0e/My9sTDbHT2a63muwqZc/lEu/xTHXECWy+f3dzD+rNwO1D5mupbvsWNL5hfnXI7v/lQg2swqB2HexVZkhDG4c32m35ZM81rluP/D5fR/YSEbf8jxOpqIFBPBKIBFQCMzq2dmMUA/4JNjB5hZNTMz3/02vvlq5/d8qFspjjcHt+XxS1uw6ttf6Pn0HMakrdeppkUkYAEXgHPuCHAzMB1YBbzrnFthZkPNbKhvWF9guZktBZ4B+vkWVSQfzIy/nlObmXd0oVvjqjwxfQ0Xjp7Hsm0/ex1NRIqwqGC8iW+1ztTjHht3zP3RwOhgzCucVS1bkrFXJvLZ8p08+PFyLhozj2s61uOOHmcQGxOU/5UiEkZ0JHAR1LNZNWbc0YXLz6nDi3M30WNkJplrf/A6logUMSqAIqpcqWgevaQ57wxpR0xkBANf/pLb31miA8hEJN9UAEVc2/qVmHprZ27p2pApy3bQbUQ6H3y1TQeQicifUgEUAyWjI7mjx5lMGdaZhMpx3PHuUga+/CVbdu33OpqIhDAVQDFyZrUyTBragYcvbMrXW36mx6gMxmVs0C6jInJCKoBiJjLCuLpDAjPuSKZzoyo8Nm01F46ex9Kt2mVURH5PBVBMVS9XihcGJjHuykR27/uVi56bx0OfrCDn1yNeRxOREKECKOZ+22X0qnZ1eXVBNuc+lcH0FTu9jiUiIUAFEAbKlozmH32a8f4NHShXKprrX1/MkNey+HbPAa+jiYiHVABhpHWdCkwe1ol7ejYmc90PdB+RwSvzNnE0V7uMioQjFUCYiY6M4IaUBnx+WxcSEyry8OSVXPzcPLL3HPU6mogUMhVAmKpTKZZXB53Ds/1bsePngzy84CD/mLxSG4lFwogKIIyZGRe0rMGsO7uQUjuKV+Zv0kZikTCiAhDKlYrm6qYlfreRePCri9j2k44kFinOVADyX79tJP5778bMW7+Lc5/K5HkdSSxSbKkA5HeiIyMYktyAmXd2oWPDyjw6bTXnPzOXrOzdXkcTkSBTAcgJ1SxfihevTmL8VYnsPXiYvuMWcM+kZfyk002LFBsqADmlHk3zjiS+Prk+k77aRtcR6bybtZVcHTsgUuSpAORPxZWI4r7eZ/HpLZ1oUKU0d09axuXjF7Bm516vo4lIAFQAkm+Nq5Xl3evb8/ilLVj/fQ69n5nDv6euYp+OHRApklQA4peICOOv59Rm9p0pXJZYi/GZG+n+VAbTvvlWVyETKWJUAHJaKsTF8NilLXj/hg6Uj43hhje/4upXFpH94z6vo4lIPqkAJCCJdSsw+eaOPHh+E77a/BM9RmXy1Iy1HDyscwuJhLqgFICZ9TSzNWa23szuPcHzZmbP+J5fZmatgzFfCQ1RkRFc06kes+/sQs+m1Xhm1jrOHZnBrFXfeR1NRE4h4AIws0hgDNALaAL0N7Mmxw3rBTTy3YYAYwOdr4SeqmVL8kz/Vrx1XVtKREVy7atZDH41i627dUoJkVAUjCWANsB659xG59whYCLQ57gxfYDXXJ6FQHkzqx6EeUsI6tCgMlNv6cy9vRozf8OPdH8qg2dnrdNqIZEQExWE96gJbD1mehvQNh9jagLfBmH+J2VmBfn2kg+RZSpToeu1jDiSy2PvpLN71vMc3LjY61giIS0tLa1Q5hOMAggqMxtC3moi4uPjSU9P9zaQBOTo3h/58eP/kLN0OhW7X0/8ZQ+zf+0Cds96gaO/fO91PJGQlJOTUyjffcEogO1A7WOma/ke83cMAM658cB4gKSkJJeSknLawUJxv/T09HQC+ZkKSmHkOnQkl5fmbuKZ6EgqNOnIjSkNub5LfUpGR3qa63Qol3+Uyz+FlSsY2wAWAY3MrJ6ZxQD9gE+OG/MJMNC3N1A7YI9zrkBX/0joiYnKuxzlrDu70P2seEbOXEuPkZnaW0jEIwEXgHPuCHAzMB1YBbzrnFthZkPNbKhv2FRgI7AeeAG4MdD5StFVo3wpxgxozZuD2xITFcG1r2ZxzQQdRCZS2IKyDcA5N5W8L/ljHxt3zH0H3BSMeUnx0bFh3t5CE+Zv4umZ6+gxMpMhyfW5MbUBsTEht3lKpNjRkcDiqZiovAvQpA1P4fwW1Rmdtp7uIzL4dJnOLSRS0FQAEhKqli3JU5efzaSh7SkfG8NNb33FFS98wfa9uhylSEFRAUhISUqoyORhnfjXRc1Y+e0v/N/8Azz0yQr27D/sdTSRYkcFICEnMsK4sl1d0oenkFIritcWZJM6Ip23v9zCUV2JTCRoVAASsirExTCwaQkmD+tEgypx3PfBN/QZM5fFm3WBepFgUAFIyGtaoxzvXt+ep/udzY97D3Hp2AXcNvFrdu456HU0kSJNBSBFgpnR5+yazLqzCzelNmDq8p10HZHOmLT1OsmcyGlSAUiRElciirvOa8zM27vQqWFlnpi+hh4jM5m+Yqd2GxXxkwpAiqQ6lWIZPzCJN65tS8noCK5/fTFXvvQFa3bu9TqaSJGhApAirVOjvKOJ/9GnKcu3/0KvpzN58OPl/LTvkNfRREKeCkCKvKjICAa2TyB9eApXtavLm19sIeXJdF6Zt4nDR3UgmcjJqACk2KgQF8PDfZox7dbONK9Zjocnr6TX03NIX6PrDoiciApAip0z4svw+rVteGFgEkeO5vK3VxYx6JUvWf99jtfRREKKCkCKJTPj3CbxTL89mb/3bkxW9k/0HJXJw5NX8PN+bR8QARWAFHMloiLzzjZ6Vwp/Pac2r87PJuXJdCZo+4CICkDCQ+XSJfj3xc359JbONK1Rlocmr6TnqExmr/5Oxw9I2FIBSFg5q3pZ3ri2LS8MTCLXwTUTshj48pc6fkDCkgpAws5/tw/clsz/nd+EpVt/ptfTmfz9w2/4MedXr+OJFBoVgIStmKgIru1Uj4y7UhnYPoF3F20l5Yl0xqZv0PmFJCyoACTsVYiL4aELmzL99mTa1a/Efz5bTbcRGXyydIe2D0ixpgIQ8WlQpTQvXp3EW4PbUq5UNLe8/TUXPzdf1x+QYksFIHKcDg0rM3lYJx7v24IdPx/g0rELuPHNxWzetc/raCJBFeV1AJFQFBlh/DWpNue3qM74zI08n7GRGSu/Y2D7BIZ1beh1PJGgUAGInEJsTBS3dT+D/m3q8NTna3l53iYmLd5G77pG+yNHKREV6XVEkdMW0CogM6toZjPMbJ3vvxVOMi7bzL4xsyVmlhXIPEW8EF+2JP/p24Kpt3SmZe3yvL36EOc+lcmUZdpQLEVXoNsA7gVmOecaAbN80yeT6pw72zmXFOA8RTxzVvWyvHZNG4YnlSA2JpKb38rbULwoWxuKpegJtAD6AK/67r8KXBTg+4kUCc0qR/HpLZ15vG8Lvt1zgMvGLWDIa1ls+EFnHJWiI9ACiHfOfeu7vxOIP8k4B8w0s8VmNiTAeYqEhN82FKcPT2V4jzOYv2EXPUZm8sBH3/DDXh1RLKHP/mz9pZnNBKqd4Kn7gVedc+WPGfuTc+4P2wHMrKZzbruZVQVmAMOcc5knmd8QYAhAfHx84sSJE/P9w/wmNTUVgLS0NL9fW9BycnIoXbq01zH+QLn8c6Jcv/zq+HjDIdK3HiE6AnrVi+a8hGhKRpmnuUKBcvknkFypqamL872q3Tl32jdgDVDdd786sCYfr3kIGJ6f909MTHSng7wljtN6bUFLS0vzOsIJKZd/TpVr4w857oY3slzde6a4xH/OcK8vyHaHjhz1PJeXlMs/geQCslw+v8MDXQX0CXC17/7VwMfHDzCzODMr89t9oAewPMD5ioSsepXjeG5AIu/f0IGESrE88NFyzhuZyWfLd2qPIQkpgRbAY8C5ZrYO6O6bxsxqmNlU35h4YK6ZLQW+BD51zn0W4HxFQl5i3Qq8N7Q9469KxAyGvrGYS8dqjyEJHQEdCOac2wV0O8HjO4DevvsbgZaBzEekqDIzejStRtfGVXlv8TZGzVzLZeMW0P2sqtzdszFnxJfxOqKEMZ0LSKQQREVG0L9NHdKHp3LXeWfyxcbd9ByVyV3vLWXHzwe8jidhSgUgUohKxURyU2pDMu9O5ZqO9fh4yQ5Snkzn31NX6WL1UuhUACIeqBAXwwPnN2H28C5c0KIGL8zZSOfH0xiTtp79h454HU/ChApAxEO1KsQy4q8t+ezWZNrWq8gT09fQ5Yl0Xl+4mcNHc72OJ8WcCkAkBJxZrQwvXn0Ok4a2J6FSLP/30XK6P5V3VbLcXO06KgVDBSASQpISKvLu9e15+W9JlIqO5Ja3v+b8Z+eStuZ7HUMgQacCEAkxZkbXxvF8ektnRl1+Nnt/PcygVxZx+fMLydIxBBJEKgCREBUZYVzUqiaz7kjhn32asmnXPvqOW8CgV75kxY49XseTYkAFIBLiYqIiuKp9Ahl3pXB3zzNZvPkn/vLMXIa9/TUbdfppCYAKQKSIiI2J4saUhsy5pys3pzZk1qrvOHdkJve+v4ztOphMToOuCSxSxJQrFc3w887k6g4JjElbz1tfbOGDr7ZzRds6tCqhDcWSfyoAkSKqSpkSPHRhU65Lrs+zs9bx+sLNvGWO1W411yfXp3xsjNcRJcRpFZBIEVezfCkeu7QFM+/oQuuqkYzL2EDn/6Tx9Mx17D142Ot4EsJUACLFRL3KcQxtWZJpt3amfYNKjJy5luTH03g+YwMHDh31Op6EIBWASDHTuFpZxg9M4pObO9Kydnkenbaazo+n8fLcTRw8rCKQ/1EBiBRTLWqVZ8KgNkwa2p6GVeP4x5SVpDyRzhsLN3PoiM4zJCoAkWIvKaEiE4e0563BbalZoRQPfLSc1CfTeWfRFp1wLsypAETCRIeGlZk0tD0TBp1D5dIx3PP+N3QbkcGkxds4oiIISyoAkTBiZqScWZWPburIS1cnUaZkFMPfW8q5IzP56OvtHNWZR8OKCkAkDJkZ3c6KZ8qwTjx/VSIloiK47Z0l9BiZwcdLVAThQgUgEsbMjPOaVmPqLZ15bkBroiIiuHXiEs4blcknS3eoCIo5FYCIEBFh9G5enWm3dmbMFa2JMLjl7a/pOSqTySqCYksFICL/FRFh/KVFdT67NZln+7fCAcNUBMWWCkBE/iAiwrigZQ2m35bMM/1bAXlFoFVDxUtABWBml5nZCjPLNbOkU4zraWZrzGy9md0byDxFpPBERhgX+opg9BWt/rtqqMfIDO01VAwEugSwHLgEyDzZADOLBMYAvYAmQH8zaxLgfEWkEEVEGOe3qMFntyYz5oq8jcW3vbOEc5/K4IOvdBxBURVQATjnVjnn1vzJsDbAeufcRufcIWAi0CeQ+YqIN37bRjDt1s6MHdCaEtGR3PHuUro9lcG7WVt1ZHERUxjbAGoCW4+Z3uZ7TESKqIgIo1fz6ky9pRMvDMw7oOzuSctIfTKdt77Ywq9HdNK5osCcO/U6PDObCVQ7wVP3O+c+9o1JB4Y757JO8Pq+QE/n3GDf9FVAW+fczSeZ3xBgCEB8fHzixIkT8//T+KSmpgKQlpbm92sLWk5ODqVLl/Y6xh8ol3+U6/eccyz94SifbDjMxj25VCxp9K4XTXKtKGIiTZ+XnwLJlZqautg5d9Jtsr/jnAv4BqQDSSd5rj0w/Zjp+4D78vO+iYmJ7nQALu9HCz1paWleRzgh5fKPcp1Ybm6uy1jzves7dp6re88Ul/SvGW58xgY3bcZsT3OdjNef18kEkgvIcvn87i6MS0IuAhqZWT1gO9APuKIQ5isihczMSD6jCslnVGHhxl08O3sdj0xdReloWB+xjoEdEihbMtrrmOIT6G6gF5vZNvL+yv/UzKb7Hq9hZlMBnHNHgJuB6cAq4F3n3IrAYotIqGtXvxJvDm7H+zd0oEH5SJ78fC0dH5vNk9PXsHvfIa/jCQFeFN459yHw4Qke3wH0PmZ6KjA1kHmJSNGUWLcCtyeWpHKjVoxJW8+Y9PW8NHcTA9rW4brk+sSXLel1xLBVGKuARERoVrMcY69MZN13exmbvoFX5mfz2oLN9E2qxdDkBtSpFOt1xLCjU0GISKFqFF+Gpy4/m7Q7U+ibVItJWdtIHZHO7e8sYe13e72OF1ZUACLiiTqVYvn3xc2Zc08qgzokMH3FTnqMzOS617L4estPXscLC1oFJCKeii9bkgfOb8JNqQ2ZMD+bCfOzmbHyOzo0qMSNKQ3p2LASZuZ1zGJJSwAiEhIqxMVw+7lnMO/ertzf+yzWf5/DlS99QZ8x85j2zbfk6sRzQacCEJGQUrpEFNcl12fOPan8++Lm/HLgMDe8+RXdR2bw7qKtOs1EEKkARCQklYiK5Iq2dZh1ZwpjrmhNqehI7n5/GcmPpzE+cwM5vx7xOmKRpwIQkZAW6TsD6ZRhnXj92jY0qFKaf09dTYdHZ/HE9NX8sPdXryMWWdoILCJFgpnRuVEVOjeqwpKtP/N8xgaeS9/AC3M2cVliLa7rXJ+EynFexyxSVAAiUuScXbs8Y69MZOMPObwwZyPvZW3jrS+30KtZNa5PbkDL2uW9jlgkqABEpMiqX6U0j17SgtvPPYMJ87J5feFmpn6zk3b1K3J9cgNSzqyiXUhPQdsARKTIq1qmJHf3bMz8e7vywF/OYvOu/QyasIjzRmXyXtZWDh3RlcpORAUgIsVGmZLRDO5cn8y7U3nqry2JMOOuScvo/PhsxmVsYM+Bw15HDClaBSQixU50ZASXtK7Fxa1qkrnuR8ZnbuCxaasZPXs9/c6pzaBO9ahZvpTXMT2nAhCRYsvM6HJGFbqcUYXl2/fw4pyNvDI/m1fmZ/OX5tVpHRveB5VpFZCIhIVmNcsxql8rMu9O5ZqOCcxe/T0PLTjI5c8vYObK78LyVBMqABEJKzXLl+L+vzRh/n1dufzMGLbu3s/g17LoPjKDN7/YzIFD4bNUoAIQkbBUtmQ0vepFk3F3Kk/3O5u4mCju/3A5HR6bxYjP1/D93oNeRyxw2gYgImEtOjKCPmfX5MKWNfhy025emruJ0WnrGZexgQta1uDaTvVoWqOc1zELhApARIS8DcZt61eibf1KZP+4j1fmbeK9xdv44KvttK9fiWs61aNb46pERBSfA8u0CkhE5DgJleN4uE8zFtzbjXt7NSZ71z6uey2LriPSmTBvE/uKyZlIVQAiIidRLjaaoV0akHl3Ks/2b0X52BgemrySdo/O4pFPV7J1936vIwZEq4BERP5EdGQEF7SswQUta/DVlp94ZV42L8/L5qW5m+jRpBqDOibQpl7FInfeIRWAiIgfWtepQOs6FbivV2NeX7iZt7/cwmcrdtKkelkGdUzggpY1KBkd6XXMfAloFZCZXWZmK8ws18ySTjEu28y+MbMlZpYVyDxFREJBjfKluKdnYxbc241HL2nOkdxc7pq0jA6PzebJ6WvYuSf0dyMNdAlgOXAJ8Hw+xqY6534McH4iIiGlVEwk/dvUod85tVmwYRcvz8tmTHrebqQ9m+WtHmpdp0JIrh4KqACcc6uAkPzBREQKk5nRoWFlOjSszJZd+3ltQTbvZG1lyrJvaVazLFe3D73VQ4W1F5ADZprZYjMbUkjzFBHxRJ1KsTxwfhMW3teNf13UjF8P/2/10OOfrWb7zwe8jgiAOXfqEyCZ2Uyg2gmeut8597FvTDow3Dl3wvX7ZlbTObfdzKoCM4BhzrnMk4wdAgwBiI+PT5w4cWJ+f5b/Sk1NBSAtLc3v1xa0nJwcSpcu7XWMP1Au/yiXf8I9l3OOVbtzmbn5MF9/n3euodbxkXSrE81ZFSP+sBYlkFypqamLnXMn3Sb7h2CB3oB0ICmfYx8iryz+dGxiYqI7HeQtcZzWawtaWlqa1xFOSLn8o1z+Ua7/2bp7n3ts2ip39sPTXd17prhuI9Ldq/M3uV8OHHLujTecq1vX5Zo5V7du3rSfgCyXz+/uAl8FZGZxZlbmt/tAD/I2HouIhJ1aFWLz9h66rxtPXtaS2JhIHvx4Bf/ofz+Hrh0MmzdjzsHmzTBkCLz5ZoFlCWgjsJldDDwLVAE+NbMlzrnzzKwG8KJzrjcQD3zoW8SJAt5yzn0WYG4RkSKtZHQkfRNr0TexFl9v+Ym6idcR8+txu47u3w/33w8DBhRIhkD3AvoQ+PAEj+8AevvubwRaBjIfEZHirFWdCrBr54mf3LKlwOarcwGJiISCOnX8ezwIVAAiIqHgkUcgNvb3j8XG5j1eQFQAIiKhYMAAGD8e6tbFmUHdunnTBbT+H1QAIiKhY8AAyM4mY/ZsyM4u0C9/UAGIiIQtFYCISJhSAYiIhCkVgIhImFIBiIiEKRWAiEiYUgGIiIQpFYCISJhSAYiIhCkVgIhImFIBiIiEKRWAiEiYUgGIiIQpFYCISJhSAYiIhCkVgIhImFIBiIiEKRWAiEiYUgGIiIQpFYCISJgKqADM7AkzW21my8zsQzMrf5JxPc1sjZmtN7N7A5mniIgER6BLADOAZs5BMufSAAAFsElEQVS5FsBa4L7jB5hZJDAG6AU0AfqbWZMA5ysiIgEKqACcc5875474JhcCtU4wrA2w3jm30Tl3CJgI9AlkviIiErhgbgO4Bph2gsdrAluPmd7me0xERDxkzrlTDzCbCVQ7wVP3O+c+9o25H0gCLnHHvaGZ9QV6OucG+6avAto6524+yfyGAEMA4uPjEydOnOjfT+STk5ND6dKlT+u1BUm5/KNc/lEu/xTHXKmpqYudc0n5GuycC+gG/A1YAMSe5Pn2wPRjpu8D7svPeycmJrrTlZaWdtqvLUjK5R/l8o9y+ac45gKyXD6/vwPdC6gncDdwoXNu/0mGLQIamVk9M4sB+gGfBDJfEREJXKDbAEYDZYAZZrbEzMYBmFkNM5sK4PI2Et8MTAdWAe8651YEOF8REQlQVCAvds41PMnjO4Dex0xPBaYGMi8REQkuHQksIhKmVAAiImFKBSAiEqZUACIiYUoFICISpv70SGAvmdkPwObTfHll4McgxgkW5fKPcvlHufxTHHPVdc5Vyc/AkC6AQJhZlsvv4dCFSLn8o1z+US7/hHsurQISEQlTKgARkTBVnAtgvNcBTkK5/KNc/lEu/4R1rmK7DUBERE6tOC8BiIjIKRTpAvizi81bnmd8zy8zs9YhkquxmS0ws1/NbHhhZPIj2wDfZ/WNmc03s5YhkquPL9cSM8sys06hkOuYceeY2RHfBZA8z2VmKWa2x/d5LTGzB0Mh1zHZlpjZCjPLCIVcZnbXMZ/VcjM7amYVQyBXOTObbGZLfZ/XoKAGyO+FA0LtBkQCG4D6QAywFGhy3Jje5F2m0oB2wBchkqsqcA7wCDA8xD6zDkAF3/1eIfSZleZ/qyxbAKtDIdcx42aTd8bbvqGQC0gBphTW75YfucoDK4E6vumqoZDruPEXALNDIRfwd+A/vvtVgN1ATLAyFOUlgPxcbL4P8JrLsxAob2bVvc7lnPveObcIOFzAWU4n23zn3E++yYVArRDJleN8/wqAOKAwNl7l53cMYBjwPvB9IWTyJ1dhy0+uK4APnHNbIO/fQojkOlZ/4O0QyeWAMmZm5P0RtBs4EqwARbkA8nOxeS8uSO/FPPPL32zXkrcEVdDylcvMLjaz1cCnwDWhkMvMagIXA2MLIU++c/l08K02m2ZmTUMk1xlABTNLN7PFZjYwRHIBYGaxQE/yCj0Uco0GzgJ2AN8AtzrncoMVIKALwkjxZWap5BVAoaxrzw/n3IfAh2aWDPwT6O5xJIBRwD3Oudy8P9JCxlfkrWbJMbPewEdAI48zQd53TiLQDSgFLDCzhc65td7G+q8LgHnOud1eB/E5D1gCdAUakHf1xTnOuV+C8eZFeQlgO1D7mOlavsf8HeNFLq/kK5uZtQBeBPo453aFSq7fOOcygfpmVjkEciUBE80sG+gLPGdmF3mdyzn3i3Mux3d/KhAdIp/XNmC6c26fc+5HIBMo6B0N/Pn96kfhrP6B/OUaRN4qM+ecWw9sAhoHLUFBb+gowA0oUcBGoB7/24DS9Lgxf+H3G4G/DIVcx4x9iMLdCJyfz6wOsB7oEGK5GvK/jcCtyfuHYl7nOm78BApnI3B+Pq9qx3xebYAtofB5kbc6Y5ZvbCywHGjmdS7fuHLkrWOPK+j/h358XmOBh3z3432/95WDlaHIrgJyzh0xs98uNh8JvOycW2FmQ33PjyNvr4ze5H2h7SevTT3PZWbVgCygLJBrZreRt/U/KIt1gWQDHgQqkfeXLMARV8AnpcpnrkuBgWZ2GDgAXO58/yo8zlXo8pmrL3CDmR0h7/PqFwqfl3NulZl9BiwDcoEXnXPLvc7lG3ox8Llzbl9B5vEz1z+BCWb2DXl/yN7j8pacgkJHAouIhKmivA1AREQCoAIQEQlTKgARkTClAhARCVMqABGRMKUCEBEJUyoAEZEwpQIQEQlT/x+C8XC0B/lpAwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "def func(x):\n", " return x**2-4*x+1\n", "\n", "x = np.linspace(0, 0.8, 100) #0から2πまでの範囲を100分割したnumpy配列\n", "y = func(x)\n", "plt.plot(x, y)\n", "\n", "plt.plot(0, func(0), \"o\", color = 'r')\n", "plt.plot(0.8, func(0.8), \"o\", color = 'r')\n", "# plot([x1, x2], [y1, y2], color='k', linestyle='-', linewidth=2)\n", "plt.hlines(0, 0, 0.8, color='k', linestyle='-', linewidth=2)\n", "plt.vlines(0, -2, 1, color='k', linestyle='-', linewidth=2)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Newton法(あるいはNewton-Raphson法) \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Newton法は最初の点$x_1$から接線をひき，それが$x$軸(y=0)と交わった点を新たな点$x_2$とする．さらにそこでの接線を求めて...\n", "\n", "という操作を繰り返しながら解を求める方法である．関数の微分をdf(x)とすると，これらの間には\n", "\n", "|　　　$x_{i+1} = x_i + \\ldots$　　　|\n", "|:----|\n", "|　　　　　　　|\n", "という関係が成り立つ．\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2⋅x - 4\n", "-2.0⋅x\n", "0 -2.00000000000000\n", "1.00000000000000 -3.00000000000000\n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from sympy import *\n", "\n", "x = symbols('x')\n", "\n", "def func(x):\n", " return x**2-4*x+1\n", "\n", "def df(x):\n", " return diff(func(x), x)\n", "\n", "pprint(df(x))\n", "\n", "x1 = 1.0\n", "df(x).subs(x, x1)*(x-x1)+func(x1)\n", "\n", "def line_f(x, x1):\n", " return df(x).subs(x, x1)*(x-x1)+func(x1)\n", "\n", "\n", "pprint(line_f(x, 1.0))\n", "x0 = 0.0\n", "x1 = 1.0\n", "\n", "y0 = line_f(x, x1).subs(x, x0)\n", "y1 = line_f(x, x1).subs(x, x1)\n", "print(y0, y1)\n", "\n", "yy0 = line_f(x, x0).subs(x, x0)\n", "yy1 = line_f(x, x0).subs(x, x1)\n", "print(yy0, yy1)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VMXXwPHvpEMCgqSBlFBCb4HQWwKINOmKGMSG+FOxoHQUQcVCExGVooAFRZQuIAokdJCASC+hI0hvoUPm/WMSg75gQrLZu3tzPs+zT7Kbze4ZgyeTc8/MKK01Qggh7MPD6gCEEEI4liR2IYSwGUnsQghhM5LYhRDCZiSxCyGEzUhiF0IIm5HELoQQNiOJXQghbEYSuxBC2IyXFW8aGBiow8LCnPqeFy9exN/f36nv6Sx2HhvYe3wyNvdlxfjWr19/UmsdlNbzLEnsYWFhxMfHO/U94+LiiIqKcup7Ooudxwb2Hp+MzX1ZMT6l1IH0PE9KMUIIYTOS2IUQwmYksQshhM1IYhdCCJuRxC6EEDYjiV0IIWxGErsQQtiMJHYhhLAZSxYoOZtSCgA531UIkR3IjF0IIWxGErsQQtiMJHYhhLAZSexCCGEzktiFEMJmJLELIYTNSGIXQgibkcQuhBA2I4ldCCFsxq0Su6wcFUKItLlVYp+y9iDPfbOeQ6cvWR2KEEK4LLdK7NdvJhG38wSNRy5l5K+7uHztptUhCSGEy3GrxP5knaIsfq0BTcqFMnrxbhqNiGPuH0ekRCOEELdwq8QOUCBPDj7uFMH33WpyT04fXvzudzqOX8PWI+esDk0IIVyC2yX2FDWK5eOnF+sypG15dh+7wIMfr6DfjM2cSrxqdWhCCGEphyR2pdREpdRxpdQWR7xeenl6KGJqFCGuZzRP1C7KD/GHiBoex+fL93LtRpIzQxFCCJfhqBn7ZKCpg17rrt2T05uBD5bl51fqEVE4L+/M207TUcuI3XEcpkxhH3ATICwMpkyxKkwhhHAKhyR2rfUy4LQjXiszSgTn4ssnqzHxiUgAZvYaytWnuhJG8kAPHIBu3SS5CyFszW1r7HeilKJh6RB+fqU+Q9Z9h++1K/98wqVLMGCANcEJIYQTKEe1CiqlwoCftNbl7/D1bkA3gJCQkKpTp051yPv+lwYNG6JuMz6tFEuXLMny93eWxMREAgICrA4jy9h5fDI292XF+KKjo9drrSPTep7TDrPWWo8HxgNERkbqqKiorH/TwoVN+eVfTvvnISm0DNGlgv8+6NqdxcXF4ZT/nhax8/hkbO7Llcdnu1LMPwwZAjlz/uOhm75+jG3xP75+4zO+junJjkOnLApOCCGyhqPaHb8DVgOllFKHlVJPO+J1My0mBsaPZz+QBFCkCJ5ffE6vb97hwVa1CF8bi2dEBJPe+IwTF6T/XQhhDw4pxWitOznidbJETAxFO3cGQO/fD4AP0O6xBzjbLor5739O0PSpRFOE7nUK80SjMvh5e1oXrxBCZJK9SzFpyOPvy6Nvv0DZpfOoVSKIuk+0YU7UQ8xfukX2nxFCuK1sndhTFAsKYMLj1bg0dx45PTTVm9fl/WffI36/5a35Qghx1ySx36J6tVI0XzqDDRN/ZFOOEJ78aBGje49h/8mLVocmhBDpJon9Xzw8FE06NmbisCfoWzYnrb54n/01o/jk0zmcuXjN6vCEECJNktjvIIePJzHPtCTnrm1crNuAh3s/TptBMxm/bA9XrssBH0II1yWJPQ3B+XLTYvJwzmzeSdGyRbnSdwCftu/BnHX7SUqSC6xCCNcjiT2dShYNZvKT1ak34Hmid66mTLP6DH5xJKv3yAInIYRrkcR+lyJa1KfStrWcGDCYfPsT6DRhDa989DO7jl2wOjQhhAAksWeIh6cHtXs8SbdZY3ijTn7eGNCJlS1iGDRpOX+du5L2CwghRBaSxJ4Jft6ePP1gFTy2bKF8cE66d2/Fy698xtCfd3D+ynWrwxNCZFNO293RzvKG3Ue1+VM5umIdhROus+K7BRyYcokqzz5K55qF8fWSLQqEEM4jid2B8tetxrC6sDfgDDleeI7tq2bxZJsXeLjz/bSqVAAPD/ffIlgI4fqkFJMFinVoQf5DCYQ9/CBvTh9Kj6kbaDl6OXE7j8seNEKILCeJPav4+FDsvYGEJ2xm1CMRvPVJD5Z0H8hjY1ey8dBZq6MTQtiYlGKymIeXJ60jCnJ92kQKPPM8V96Moe+K58nXrDE9HyhF8SD7Hh0mhLCGzNidxDuiEgXWrSD/Z6N4KKIAq7cf4dnXv6XfjE0cPXfZ6vCEEDYiM3ZnUooc7dvyENB4+Wp8Wj7OD7/Xp9XKR2nbqALPNShOXn8fq6MUQrg5mbFbJG+9Wvgn7KRDuSCWfPEcM37eQP2hsXy8eDcXr96wOjwhhBuTxG6loCACJk0g1x8b+HZAK3oeXsG6cd9Rf2gsE1fs4+oN2UVSCHH3pBTjCooUoSRQsn0tOnZ/ie3bF9Lz4GN8saIkLzcKp12V+/DylN/BQoj0kWzhSpo3x2/HNiKebM/XV9cTGOBD/2kbaDJqGT9tOiLbBAsh0kVm7K7GxwdefZUCwKwzZ7hSrgKf1+zAy8ca8ul9een5QEmiSwWjlKxiFULcnszYXZjKm5ccCxfQ/fwWfp/ei+D9O3lqcjztP1vFqoSTVocnhHBRMmN3dRUqoH79ldxz5zKhWnXmr97Nt8t38ejBs9Quno/XmpS0OkIhhIuRGbs7UApatcI7fyitPU4x9etezNvzI3/t+5P2n61m5PorbD58zuoohRAuQhK7u2nTBrVtG+Xy+bJo1kB6Nwlnz9mbPDhmBd2+imfHX+etjlAIYTFJ7O4oOBjGjsVjzWqejw5nwe/j+PCev1i95xRNRy3nhW83kHBcjuoTIruSxO7OcucG4FytGrSd+D4bfvuI10t4ELvjOE0+XEaP7zey7+RFi4MUQjibJHZ3pxSn6tSBLVvwbtyIrv5nWd4riu5VAvl5y180HrmU16b9wYFTkuCFyC4ksduFry+89hp06kS+7Zt4tXsrfsuzkydrFOSnTUdoOGIpvX/8g0OnL1kdqRAii0lit6Pq1WHRInItmMvr/Tux4omydKlVhFkbjxA9PI6+0zdJghfCxhyS2JVSTZVSO5VSCUqpvo54TZFJFSrAr7/CmDEEFS/Mm75/srJ1ATrXLMKM3/+UBC+EjWU6sSulPIFPgGZAWaCTUqpsZl9XOIBSEB1tPh48SFDzxgxaMoHlz1Tm0RqFmbHBJPg+P27i4ClJ8ELYhSNm7NWBBK31Xq31NWAq0NoBryscqWtX2LYNrl4lpE8P3mpdnqW9o4ipUZiZG/8kekQcvX6Qi6xC2IEjthS4Dzh0y/3DQA0HvK4D+QIrUGoEMB9YAVyzNiQLeQAB33/PT8AaYG9APnLXaM/3Vx9g2roDXNway7nV07hx5ojFkQrhumJjY60O4Y6ctleMUqob0A0gJCSEuLg4Z701Jok/h6kWDQFWA68CbYH1wEEnxmK9JOA8MBxTQ9uZeIoXF4/n4JofyF2jPQGVm+JfLpqL25aaBH/6sLUBC+GCEhMTnZzH7oLWOlM3oBaw8Jb7/YB+//U9VatW1c4EaDNUIynJ3J54QuvAQK3LlNH6o49Sv+ZuYmNjM/7NV65oPWyY1rt2aX3ypNanT+tj5y/rIfO26dKvL9BhfX/SL0xZr3ccPe+weO9Wpsbn4mRs7suK8QHxOh152RE19nVAuFKqqFLKB3gEmOOA180ySpnbpElw7BhMngzly5uv1asHrVrBZ5/B/v1WRukkvr7QsyeEh8OCBVC6NMFfT6R/k3BW9Inmfw2KE7vjOA+MWsazX8ez5U/ZbEwIV5fpxK61vgF0BxYC24FpWuutmX1dZ/HwMG3fDRua+7NmwSOPwKpV0Lu3eWz6dFi0CK5etS5Op+jc2bRITp8OtWqRL4cXfZqWZmXfhrzUKJxVe07R8uMVPDV5HRsOnrE6WiHEHTikxq61no+5Kun2AgPh0UfNLcWxYzB8uGkqadUKvv4arl0zhx3ZTsWK5rfY7t3g6QkjR5KnRQtevb8UXesV5atV+/lixT7afbqKOiXy0T06nJrF7pUTnYRwIbLyNB2efx5Wr4Y9e0zXYMpjZcuaVfyLF0NSkrUxOpRSULIk6OQzVuvWhR49yH05ke4Nw1nRpyH9m5dm17FEOk1YQ4exq4ndcTzlGosQwmKS2O9CYCA0aGA+Hz8evvwS7rkHRowwj82bB2PHwoED1sXoUErBq6/C1q1w+bK5GAH4eym61S/O8t7RvN26HH+du8KTk9fR8uMVzN98lJty6LYQlpLEnkEeHlCtGgwcCPPnm/v+/rByJURGmtn82bNw/rwNavPJ+7/TowesW2fKNb/8gp+3J4/VCiO2ZxTDOlTk8rWbPD9lA/d/uJQf4g9x/aad/owRwn1IYnegqChTfz92DKZMgTx5YNo0kxdbtza58Zy7N5VERsK778ILL0DLlnDiBD5eHjwUWYhfX23AmEcj8PPypNePm2gwNJbJK/dx+dpNq6MWIluRxJ4FPDwgIsJ83rWrqc137Ghm81evwtq1qbV5t5vNK2V+S23ZAg8+aH57bdoEZ87g6aFoWbEA816qy6Qnq1EgTw4Gzd1G3Q+W8ElsAucuX7c6eiGyBUnsTpDSafP112b2HhJiavMDBpj7a9fCpUtuVpv39YVnnwVvb5g7F0qXhk8+gRs3UEoRXSqYH5+rzbRna1Gh4D0MW7iTOu8v4b0F2zl+/orV0Qtha5LYLRAWZmrza9aY2XylSqaVslo1U5vv2RN27bI6yrswYEBq/3vnzv/4UvWi9zL5yerMe6ku0aWDmbBsL3U/iKXfjM1ybJ8QWUQSu8UCA8HPz5Su//rLNJ7kygWJiXD0qKl6jBsHB119O5uKFU1t6ZNPTAdN166wc+ffXy5X4B4+7hRBbM8oOkQWZPqGwzQcEcfzU9az6fBZCwMXwn4ksbuQlFWwb74JVaqYBN+xI6xYAVWrmm0OtIZly8wCKZejFOTLZxY2lS4NdeqYTpozqatUi+Tz5922FVjRJ5rnGhRn+e6TtBqzkkcnrGHZrhPSCy+EA0hid2EBAam1+WPH4IknTPtk794QFARt2sDChVZHeRs+PqaetG2buXiwcydcuQI3bvz9lOBcfvRuWppVfRvSr1lpEo4n0mXib7QYvYLZG//khrRKCpFhktjdhIcH5MhhLrquWQMJCfDQQ2ZynJRkNi/r2dNUQ1xmNh8cbOpINWua306VKsEvv/zjKbn8vHm2QXGW94nmg/YVuHrjJi9P3UjU8DgmrdzHpWs37vDiQog7kcTupoKCICYGGjc290eMMKWb/v3huefMY/Pnu1CnTdeupv/9+efNhYOb/+xt9/XypGO1wvzaowETukQSmtuPwXO3Ufv9JUzffY0TF9ytL1QI60hit4Fba/Nr18KECebxOXNSV8EOG2Yes6yEndL/vnWrSe6enjB79j/q7wAeHor7y4bw43O1mf5cLaqH3ctPe65T54Ml9JuxiT0nEi0agBDuQxK7DXkk/1THjk3ttClXzjzWsqWpzVvWaePrCw88YH7DLF5sLrJ++uk/6u8pqha5l/FdInmvXg46VC3I9A1/0mjEUrp+uY61e0/JhVYh7kASu815eprZfPPm5v6XX8LDD5tOm2eeMY/Nn29BbV4pGD06tf/9gw/u+NRQfw/ebVuBVcn7wq8/cIaO49fQ+pOVzPnjiFxoFeJfJLFnM7eugk3pqDlyxKwxCgoyF2S1vu0EOmuk7P/eq5fZpqBVq3/0v/8j9gBfXr2/JKv6NuKdNuW5cOUGL333Ow2GxfH58r1cuCJbFggBktgF5rpmyirYZ54xk+levUz5pmdPWLIkixO9UqZFMjwc6tdP7X8/f/62T8/h40nnmkVY/GoDxj9Wlfvy5uCdedup/d4S3vlpG4fPXMrCYIVwfZLYxd8CA6FJE/P58OHmTNhcueDtt01LZWxsFtfmU85f3bbNXChQCo4eRd28/e6QHh6KJuVCmfZsLeZ0r0N06WAmrdpP/aGxvDBlA+sPyPF9InuSxC5uK6U2/+abJqH7+Jj9vpYvN6tgy5UzJZxLl7Jgh8rg4NT+zWHDiOza1dTi/0PFgnkY3SmC5b2jeaZ+MZbtPkH7z1bR5pOVzJU6vMhmJLGLdKtbF775xnTaTJoEoaGmpTJlv/lx4/5f92LmjRjBvqefNs35KQ36/6FAnhz0a1aGNf0aMbhVOc5eusaL3/1O/aGxjF26h3OXpA4v7E8Su7hrKbN5Dw945BGzCvbhh81s/vRpcw3UYatgleJk3bqm//2ll8zCpqFD0/wN4u/rxeO1w1j8WhQTukRSJJ8/7y/YQc33FvP6rM0kHJd+eGFfkthFpqWsgv3mGyhe3Gx7kCtXaqfNkiUmwWeqNu/rC2XKmD1n9u79z/73W3kmL3j6rltN5r9UjxYV8zNt3WEaj1zK4xN/I27ncZLkjFZhM5LYhcMVKmRq8ymdNtWrm48ptfmU66MZ4u9vVl79+iv8+KP5syCdC5XKFsjN8IcqsapfQ169vyTbjp7niUnraPzhUr5avZ+LV2VfGmEPkthFlgoMNLtUlimTWpvPlQtOnjTdjO3aZbDTJmX/9yZNYOJEc0zfHfrf/19MAb681CiclX0aMqpjZXL5ejFw9lZqvruYt+ZuY78cACLcnCR24TS3dtrUrw9eXtC+fWqnTcp+NmvWpLM2r5S5de4MDRqYq7uvvprupnsfLw/aRNzHrBfqMP252kSXDuar1fuJHhHHU5PXsXTXCSnTCLckiV1YJmfO1Nr8X3+ZpperV+GVV1L3m587Nx0vlNL/vnWrKfJ7ecHvv6c7wSulqFokL6M7RbCyb0NebBjOpsPneHzibzQeuZTJK/fJqlbhViSxC5fg6WlKNr6+/9xvXmtza948HZ02wcHwwgvmGwYPvu3+72kJye3Hq/eXZGXfaEZ1rEzuHN4MmruNmu8u5o1ZW9h97ELmBytEFpPELlxSSqdNq1bm/qBBJvH37596XvaSJXeozSsFM2fCkCFmi+CUGs9d8PXy/LtMM/uFOjQtn5/v4w9x/4fL6DR+DT9vOSqLnoTLksQuXJ5SpjY/aJDZb/7bb83js2endtoMGXKbb2rTxpRnnnoKDh0yU/4MrKCqVCgPIx6uxOq+DendtBQHT1/if99soN7QWEYv3s3xC1cyPUYhHEkSu3A7Xl7m40cfpXbalCljHouJ+dd+876+5oDtHDngwgXT//7ZZxna1SxfgC/PR5VgWe9oxj9WlRLBAYz8dRe131tC9283yB7xwmVkKrErpR5SSm1VSiUppSIdFZQQ6ZXSadOunbk/apSpzS9fnlqjj42FJZsCufbxOFNznz8fTpwwi50y8p7Jm499/XQNYntG8XjtMJbtOkHH8Wt4YNQyvlq9Xy62Cktldsa+BWgHLHNALEJk2q2rYNeuNRWZgwehXz/ztbaDKnF9xlxuBuc3rZGtWsGuXRl+v6KB/rzRsixr+zdmaPuK+Hp5MnD2Vmq8u5h+Mzaz9cg5B45OiPTJVGLXWm/XWqdvVYgQFnn8cZPkExJMud3bG956CyLiPuSn8/W5Vr0ON98YlKn3yOHjycPVCjGnu7nY2qJCfmb+fpgWo1fQ5pOV/BB/iMvXbr/9sBCOJjV2kW0EBZkFqgADB8LYSb7ER/WkbfhWrlSrx29rNQufm8XBvRnfWkApRaVCeRj2UCXW9mvMwJZlSbx6g14/bqLGu4sYNGertEyKLKfSutijlFoEhN7mSwO01rOTnxMH9NRax//H63QDugGEhIRUnTp1akZjvmvR0dEAxMbGOu09nSkxMZGAgACrw8gyzhrfrg2eVHunH77nzvBu4Ps0GxlEUNBVlAJv74xfFNVas+tMEksOXif+2E1uaiiZ14OoQt6U9r/CvffY82cn/y4dLzo6er3WOs3rmWkm9vRIT2K/VWRkpI6PT9dTHUIpBWDbjoW4uDiioqKsDiPLOHV8WnNz5myuv9wTzzGjWejZnJgYiI6GZs3MFgiBgRl/+VOJV/lx/WG+/e0gB05dwt8bOlYvSqfqhQgPyeW4cbgA+XfpeEqpdCV2KcUIcSul8GzXBr+ErXi3fICW12ZwtFMPYpqfYflyOHwY9u83Z8KmbEd8N/IF+PJsg+LEvhbFlK41KJfPk6/X7Of+D5fR4bNVTF9/WGrxItMy2+7YVil1GKgFzFNKLXRMWEJYzNfX9FLWqUPOpIs89EZpvqn9KZUr3MTHx+wenNJp89NP5kzYQ4fS//IeHoo6JQJ5vrIfq/s1ol+z0py6eI3XfviD6u8uYuDsLWw7cvvDvIVIi1dmvllrPROY6aBYhHA9ISEwfrzZmmDiRFCKAt4nGDQoiEGDTDu8t7dJ6lWrmuMCmzWDLl2gQoX0vUVg8iy+W/1irNl7mu/XHWTqukN8tfoAFQveQ8dqhWhVqQC5/LyzdKjCPqQUI0R6VK4Mo0eb8wA7dDD977t3ExQEefJAkSJw7Bh88YWZzR8+bMo0HTua3wvpmc0rpahVPB+jHongt/6NePPBsly9nsSAmVuoPmQxPX/4g3X7T9v2WpFwHEnsQtytX36BevWgVi2TtZN5ekKNGmZPm2bNzKrXBx+EpUshIsLsQw+wYUPatfk8OX14sk5Rfn6lHrNeqEObiAIs2HyUh8auptHIpYxdukf2qBF3lKlSjBDZkq+vuXrapQskJsLZszBtmln95OX1j6d17mxuN2+apyYlwcsvw6ZNptOmSpVA/quxQilF5UJ5qFwoD6+3KMu8zUeZtu4Q7y/YwbCFO4kuFczDkQWJLh2Mt6fM04QhiV2IjAoJMbcDB2DqVPj4Y/jwQ2jc+P891dPTHPINZh+bEydg4ULYtcsk406doGBBM9OvWxd8fP7/2/n7evFwZCEejizEnhOJ/BB/mOkbDrNo+zECA3xoG3EfD0UWoqTN2ibF3ZNf8UJkVpEi5gSQt9+GF180BfU06uBBQWYm37DhcQB69EjttEnZg37lyjvX5osHBdC3WWlW9W3I510iqVokL5NW7qfJh8toPWYFX685wLlLshFZdiUzdiEcIWX/91atzAXWnj1N3WXgQHN1NQ3Vq6fuOX89OR/Pnm0acUJDzfXaQYP+//d5e3rQuGwIjcuGcDLxKrM3HuGH+EO8MWsLb/+0jSZlQ+hQtSD1woPw9FAOHbJwXTJjF8KRPJL/l+rVCy5ehFKl4PPP7+olvJO7GocOTe20KVHCPPbCC9C27e07bQIDfHm6blEWvFyPn16sy6PVC7My4SRPTFpHrfcW89787eySfWqyBUnsQmSFkBBz2scvv5gCO5jTnO5SSqdNynGAgwaZbQ2WLoUmTcwfBatXmz3nUzptlFKUv+8eBrUqx9r+jRnbuQoVC+bhixX7aPLhMh78eAWTV+7jVOJVx4xVuBwpxQiRlSpVMrfERFOqKVMGhg+HkiUz9HIptfnOnU0ZXymzxcGHH5pt5aOj4euvIWdO88eDj5cHTcvnp2n5/H+XamZsOMygudt4Z952okoF077KfTQsE4yvl6djxy4sI4ldCGcICIAtW8wip9q1zcGtTZpk6iWT97ajUydzO37czOT9/WHkSHNkYLNm5lavXmqp5um6Rdnx13lmbviTmb//yaLtx8jt50XLSgVoF3EfVYvk/XvjPOGeJLEL4Sy39r/nygUrV1Jg9mzT3+iV+f8Vg4PNcYAAr7xiXnb+fOjbF2bNgtOnYdUqk+hLF8pNv+a56d20NCsTTjJjw2FmbDjMt2sPUujeHLStfB+tI+6jeJB9t921M6mxC+FsISGmVpInD0FxcWZZ6qJFDn2LlNr84MHw229QoIDptklZBVu+vPkDIummomZYEKMeiSD+9fsZ8VAlwvL583FsAo1GLKX1mBVMWrmPExekHu9OZMYuhFXKleOPkSOJOnsWhg2DBg3M495Zs9lXRARMmWJWwcbHQ1iY6ZVv0yZlv3kv2rYtSPuqBfnr3BXm/nGEWRv/ZHByPb5OiUDaVC5Ak3KhBPhK6nBl8tMRwkpKmf7Ftm3N/QYNIDIS3ngjXf3vGZEymweIioLdu80q2AULoGxZ8/jIYX40a1aMmf8rxoEzF5i18U9mbzzCq9P+wM97M43KhNC6UgEalAqSi64uSBK7EK5k2jST1EuVgk8+MSuTstitnTZgLsLmyAF9+phOm3HjctGrY2keLV+Ko9fPMHvjEeZtPsq8TUfJ7edF8wr5aVWpADWK5ZNFUC5CErsQruTW/d+vXjW3VatMrcRJgoNNbX7wYJPklYKTJ6FSJUX+/PfSrNm9jOtYlkT/k8z54whz/jjC1HWHCM7lS4uK+XmwUgEiCmXNXxsifSSxC+GKKlc2H3ftgmeeMTWSESMgPNypYQQHp35+7BisW2dKNvv3etC+fTCT3wvm+RpJ+Jc4zppjh5my9iCTVu6nYN4cVMpzg+CS5ymTP5e0TzqZJHYhXFnJkmbFakr/e3y82XTMAp6eULOmuYHpsrn/fliwwIOFA0N5/PFQ4odcZ8KsU/x+/hAL9h1n3ujlFAvyp2WF/LSsVEB2nnQSSexCuLqU/vdnn4XcuWHUKPPYM884pP89o7y9/7nf/NmzkMvXmwWfhxIfH0qZcieIbA4nc+3h49gERi9JoGRIAC0rFqBFxfzSI5+FpI9dCHeRO7f5GBUF33+fJf3vGeXpCfnymXr8okWmgtQw6iQR+YP4rltN6h17gLKH6nLtYCAjF+6m0YilNB21jI8X72bPiUSrw7cdSexCuJvKlc2uX2+9BXFx5rGzZy0N6d+Cg+H++4/x1FPmfvf/eVGtxD1cWlmO4GVNGdiyLNf/uocPph/8O8mPXrybhOOS5B1BSjFCuKNb+98vXoRy5eCRR7K0/z0zUvabHzwYrl71wNe3KAd+hd+/1+TMc5394ccZcfQPRv66i5IhATQrn5/mFfJTMiRALrxmgMzYhXB3/v7mhOzz56F0aVi2zOqI/pOvr/k4eDCcOK74cYoPr7UvyNr+jSh/pA5bv6zAkJFXaPTWbzQasZShP+9g8+FLI4e1AAAOHUlEQVRz6DROpRKpZMYuhB2EhMCECbBxozk8dedOcxLHbc5fdSX/7LTxY+JwPxYuhJlz7mHZDE3wm+sZPe0kI788S7Hyl2leKYSm5UOpUjgvHrIY6o4ksQthJyn97zt2mC6a8uXN/u9O7n/PqOBgeOwxeOwxz+T95mvwdZ7rDH4nifWzPdlU6CTjHognOMiDB8qH0LRcKDWL5cPHS4oPt5L/GkLYUd26sG0b1KkDrVunHqTqRlJK64918iZhqy+H9nsxqm8go58oT64D4Yx8NowHHztPqafX0/2bjczffJSLV29YG7SLkMQuhF35+kLv3rBpk2k6f/hhGDsWbrhn8gsOhq5PetG2SgEWjyvCwhk5aVU1hCsry7F4bSLPjNlBsXbb6ThyI1N/O5ittxqWUowQdpeyiKl/f3MCxyefwJgxqdsEuyFPT2hQ15MGdc0ipxs3a/P9L+d5/32Y8UZOZua8QmCLddSs5kHDUiE0rRhCieDssyBKErsQ2UVK//usWabvPSkJDhyAokWtjizTvDw9iGmWh5hmcOOG5ocF19l3PT8L1pzjpdcK41fkFIUr/km7Noo2tQOpWiSvrXeilFKMENlJSv9769awfTtUqwavveZyC5wyw8tL0elBf/q3K87yoVXYuPkGj3X05OqhQCbPOUeHj38jf/39dBiYwKz1R7hwxf2uP6RFErsQ2VW5cuZ8vJT+9/37rY4oS1QMz8G4t4LYtyYfO7+szNC2lSkWkoP5XwTRrk4QJTpsJ+bzNYyae4ADpy5aHa5DZCqxK6WGKaV2KKU2KaVmKqVcb8mbEOLOQkNN/3tcnNk1cvp0WLzY6qiyTC4/bx6pF8qa6aFcOJSb+csTeeqRHBw9eZ3XHgklvPRNCkcf5PmP9rJqz0mu30yyOuQMyeyM/VegvNa6IrAL6Jf5kIQQTle6tCnT+PmZXSPbtDFn5tmYp4eiabW8vNc5nCV965Gw7yZ93kkkT4AH039N5NEJawmteZiorgcYO+8Ixy9cMd84ZQqEhdGgYUNzcOyUKZaO43Yyldi11r9orVN6p9YABTMfkhDCMi1amP732rVh0CDzWJJ7zlrvVtHgnLz9bAE2zS3I3ullGRtTlVq1YPM6X154KJDi9Y8xrMtAtjz5AdcP/InS2lx87tbN5ZK7I7tingK+d+DrCSGs4Odn+t8BLl2CKlWgRw94+mlL9393Jn9fL5pWCKXpaNBas/nweeavS+LxZ8bz9PUJeHKTWSQfQH7pEjf79cMzJsbaoG+h0tpYRym1CAi9zZcGaK1nJz9nABAJtNN3eEGlVDegG0BISEjVqVOnZibuuxKdfF5kbGys097TmRITEwkIsG+Prp3H5w5jC0hIoMSYMXhduMCOPn1ILFkyXd/nDmO7Ww0aNkRpzXW88CZ1oVcSiugh8yifz5MKQZ6UyOOBVxa0U0ZHR6/XWkem9bw0E3uaL6DUE8CzQCOt9aX0fE9kZKSOj4/P1PvejZRtP+26O1xcXBxRUVFWh5Fl7Dw+txmb1jBzpumFz5kTLlxIc/8Ztxnb3QgLM+WXfzkfXICug39kw8Ez3EjS+Pt4Uqt4IA1KBlIvPIiwQH+HvL1SKl2JPVN/VymlmgK9gQbpTepCCDekFLRrZz6fPx+6dIEnnoDXX3fJ/d+zzJAhpqZ+6ZZ0lzMnuUcOZVpMLS5cuc6qPadYtusEy3afYNH2YwAUujcH9cKDqB8eSJ0SgeTy887SMDNbMBsD+AK/Js+K12it/5fpqIQQrqt5c9P//sYbULWqWejk42N1VM6RUkcfMAB98CCqcGGT7JMfz+XnzQPlQnmgXChaaw6cusSy3SdYvvskczYe4du1B/m8SySNy4ZkaZiZSuxa6xKOCkQI4UZS+t9PnjRJvX9/aNTI3OwuJgZiYliaRqlJKUVYoD9hgf50qRXG9ZtJbDx0lnIFcmd5iLLyVAiRcYGB5mNkpOl/b93a9v3vGeXt6UG1sHvJ6ZP1nUWS2IUQmdeuXWr/+8KFAHhcvmxxUNmXJHYhhGP4+UGfPtC9O2zbRs2YGBg3Dm7etDqybEcSuxDC8cqWZdMHH8C330JEBBw8aHVE2Ur2WEYmhHC6xPBws7nY/PlQoACsWgVBQW5z/qo7kxm7ECLrKGX2n/HyMgds16oFvXrBuXNWR2ZrktiFEM7x1FOm//3sWejUyepobE0SuxDCeVL632fOhMuXzWInG+//bhVJ7EII5/P1NV00Xbum7v++d6/VUdmGJHYhhDVS9p9J6X8/dQoSE211/qpVJLELIayVsv97tWrw88/mNKdx4+DGjbS/V9yWJHYhhOvo0MG0R377rUn0165ZHZFbksQuhHAtVaqY/veJE80GYxMnQkKC1VG5FUnsQgjXo5RZsQpw+jTUrCn973dBErsQwrX17Gn638+cgeHDzWM2PQ3NUSSxCyFcX2gofP45vPWWWcFapQosWWJ1VC5LErsQwn0oBaVKmSP5unY1/e9//WV1VC5HErsQwr0oBe3bm/73Bg0gIMBcXJX6+98ksQsh3JOfH/ToYRL7jBlmJi/7vwOS2IUQdtC7NyxYYPrf27WzOhrLSWIXQthDRITpf//sMzNr79Ej256/KoldCGEfSplDPW7ehJCQbLv/uyR2IYT9+PhA376p/e8rV5q9Z7JJ/V0SuxDCvlL635s3hx9+MP3vsbFWR5Xl5MxTIUT28MgjZib/9NNQqRJMnWr2hbchmbELIbKHW/vfY2JMUl+82Jb1d0nsQojsxc/PbA8MMGeO2f99/Hhb1d8lsQshsq+PPoJ582DKFBgwwOpoHEYSuxAie0vZ/33gQDhwAB5+2O33f5fELoQQSkHOnKb3vUqV1P3fz5+3OrIMyVRiV0q9rZTapJTaqJT6RSlVwFGBCSGE0/n5pfa/X7wIly+bgz7crP6e2Rn7MK11Ra11ZeAnYKADYhJCCGuFhsKnn5oZ/PDhbrf/e6b62LXWt/6d4g/IsSZCCHsZMsQk9qefhuhocwari8v0AiWl1BCgC3AOiM50REII4UqUMu2RLVuaEo3W8NVXUKSI1ZHdkdJpnB2olFoEhN7mSwO01rNveV4/wE9r/eYdXqcb0A0gJCSk6tSpUzMc9N2Kjja/b2JtupQ4MTGRgIAAq8PIMnYen4zN/Xhcu0be+HgOVKzo9PFFR0ev11pHpvW8NBN7eimlCgPztdbl03puZGSkjo+Pd8j7podSCgBHjdXVxMXFERUVZXUYWcbO45OxuS8rxqeUSldiz2xXTPgtd1sDOzLzekIIITIvszX295VSpYAk4ADwv8yHJIQQIjMy2xXT3lGBCCGEcAxZeSqEEDYjiV0IIWxGErsQQtiMJHYhhLAZSexCCGEzktiFEMJmJLELIYTNSGIXQgibkcQuhBA2I4ldCCFsRhK7EELYjCR2IYSwGUnsQghhM5LYhRDCZiSxCyGEzWT6MGt3oLUmLi7O6jCEEMIpZMYuhBA2I4ldCCFsRhK7EELYjCR2IYSwGUnsQghhM5LYhRDCZiSxCyGEzUhiF0IIm5HELoQQNqO01s5/U6VOAAec/LaBwEknv6ez2HlsYO/xydjclxXjK6K1DkrrSZYkdisopeK11pFWx5EV7Dw2sPf4ZGzuy5XHJ6UYIYSwGUnsQghhM9kpsY+3OoAsZOexgb3HJ2NzXy47vmxTYxdCiOwiO83YhRAiW7BdYldKNVVK7VRKJSil+t7m60opNTr565uUUlWsiDMj0jG2mOQxbVZKrVJKVbIizoxIa2y3PK+aUuqGUqqDM+PLrPSMTykVpZTaqJTaqpRa6uwYMyod/y7vUUrNVUr9kTy2J62IMyOUUhOVUseVUlvu8HXXzCdaa9vcAE9gD1AM8AH+AMr+6znNgQWAAmoCa62O24Fjqw3kTf68mZ3GdsvzlgDzgQ5Wx+3gn10eYBtQOPl+sNVxO3Bs/YEPkj8PAk4DPlbHns7x1QeqAFvu8HWXzCd2m7FXBxK01nu11teAqUDrfz2nNfCVNtYAeZRS+Z0daAakOTat9Sqt9Znku2uAgk6OMaPS83MDeBGYDhx3ZnAOkJ7xPQrM0FofBNBau8sY0zM2DeRSSikgAJPYbzg3zIzRWi/DxHsnLplP7JbY7wMO3XL/cPJjd/scV3S3cT+NmUm4gzTHppS6D2gLfObEuBwlPT+7kkBepVScUmq9UqqL06LLnPSMbQxQBjgCbAZe1lonOSe8LOeS+SRbHGad3SilojGJva7VsTjQKKCP1jrJTPxsxwuoCjQCcgCrlVJrtNa7rA3LIR4ANgINgeLAr0qp5Vrr89aGZV92S+x/AoVuuV8w+bG7fY4rSlfcSqmKwOdAM631KSfFllnpGVskMDU5qQcCzZVSN7TWs5wTYqakZ3yHgVNa64vARaXUMqAS4OqJPT1jexJ4X5uidIJSah9QGvjNOSFmKZfMJ3YrxawDwpVSRZVSPsAjwJx/PWcO0CX5anZN4JzW+qizA82ANMemlCoMzAAec7OZXppj01oX1VqHaa3DgB+B590kqUP6/l3OBuoqpbyUUjmBGsB2J8eZEekZ20HMXyIopUKAUsBep0aZdVwyn9hqxq61vqGU6g4sxFytn6i13qqU+l/y18diOiqaAwnAJcxswuWlc2wDgXzAp8kz2xvaRTcpulU6x+a20jM+rfV2pdTPwCYgCfhca33bFjtXks6f3dvAZKXUZkz3SB+ttVvs+qiU+g6IAgKVUoeBNwFvcO18IitPhRDCZuxWihFCiGxPErsQQtiMJHYhhLAZSexCCGEzktiFEMJmJLELIYTNSGIXQgibkcQuhBA28395qq1yUKCj7AAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(x0-0.05, x1+0.05, 100)\n", "y = func(x)\n", "plt.plot(x, y)\n", "\n", "plt.plot(x0, func(x0), \"o\", color = 'r')\n", "plt.plot(x1, func(x1), \"o\", color = 'r')\n", "# plot([x1, x2], [y1, y2], color='k', linestyle='-', linewidth=2)\n", "plt.hlines(0, x0, x1, color='k', linestyle='-', linewidth=2)\n", "plt.vlines(0, -3.5,1.5, color='k', linestyle='-', linewidth=2)\n", "\n", "plt.plot([x0, x1], [y0, y1], color='b', linestyle='--', linewidth=1)\n", "plt.plot([x0, x1], [yy0, yy1], color='r', linestyle='--', linewidth=1)\n", "\n", "plt.grid()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "|$x_1$ |$f(x_1)$ | $df(x_1)$ |\n", "|:----|:----|:----|\n", "|1.0　| 　　　　 | 　　　　 |\n", "|　　　　　　　| 　　　　　　　 | 　　　　　　　 |\n", "|　　　　　　　| 　　　　　　　 | 　　　　　　　 |\n", "|　　　　　　　| 　　　　　　　 | 　　　　　　　 |\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 二分法とNewton法のコード\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 二分法(bisection) \n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x1 x2 f1 f2 \n", "0.000 0.800 1.000 -1.560\n", "0.000 0.400 1.000 -0.440\n", "0.200 0.400 0.240 -0.440\n", "0.200 0.300 0.240 -0.110\n", "0.250 0.300 0.062 -0.110\n", "0.250 0.275 0.062 -0.024\n" ] } ], "source": [ "x1, x2 = 0.0, 0.8\n", "f1, f2 = func(x1), func(x2)\n", "print('%-6s %-6s %-6s %-6s' % ('x1','x2','f1','f2'))\n", "print('%-6.3f %-6.3f %-6.3f %-6.3f' % (x1,x2,f1,f2))\n", "for i in range(0, 5):\n", " x = (x1 + x2)/2\n", " f = func(x)\n", " if (f*f1>=0.0):\n", " x1, f1 = x, f\n", " else:\n", " x2, f2 = x, f\n", " print('%-6.3f %-6.3f %-6.3f %-6.3f' % (x1,x2,f1,f2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Newton法(あるいはNewton-Raphson法) \n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from sympy import *\n", "\n", "x = symbols('x')\n", "def func(x):\n", " return x**2-4*x+1\n", "def df(x):\n", " return diff(func(x), x)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000000000 -2.0000000000000000000000000\n", "0.0000000000 1.0000000000000000000000000\n", "0.2500000000 0.0625000000000000000000000\n", "0.2678571429 0.0003188775510204081378197\n", "0.2679491900 0.0000000084726737969074341\n", "0.2679491924 0.0000000000000000059821834\n" ] } ], "source": [ "x1 = 1.0\n", "f1 = func(x1)\n", "print('%-15.10f %-24.25f' % (x1,f1))\n", "for i in range(0, 5):\n", " x1 = x1 - f1 / df(x).subs(x,x1)\n", " f1 =func(x1)\n", " print('%-15.10f %-24.25f' % (x1,f1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 収束性と安定性" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "実際のコードの出力からも分かる通り，解の収束の速さは2つの手法で極端に違う．2分法では一回の操作で解の区間が半分になる．このように繰り返しごとに誤差幅が前回の誤差幅の定数($<1$)倍になる方法は1次収束(linear convergence)するという．Newton法では関数・初期値が素直な場合($f^{\\prime}(x) <> 0$)に，収束が誤差の2乗に比例する2次収束を示す．以下はその導出をMapleで示した．\n", "\n", "\n", "```maple\n", "> restart; ff:=subs(xi-x[f]=ei,series(f(xi),xi=x[f],4));\n", "```\n", "\n", "$$\n", "{\\it ff}\\, := \\,f \\left( x_{{f}} \\right) +D \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei}+\\frac{1}{2}\\, D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2} +\\frac{1}{6}\\, \n", "D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{3}+O \\left( {{\\it ei}}^{4} \\right)\n", "$$\n", "```maple\n", "> dff:=subs({0=x[f],x=ei},series(diff(f(x),x),x,3));\n", "```\n", "$$\n", "{\\it dff}\\, := \\,D \\left( f \\right) \\left( x_{{f}} \\right) + \n", "D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei}+\n", "\\frac{1}{2}\\, D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2} +O \\left( {{\\it ei}}^{3} \\right)\n", "$$\n", "```maple\n", "> ei1:=ei-ff/dff;\n", "```\n", "$$\n", "{\\it ei1}\\, := \\,{\\it ei}-{\\frac {f \\left( x_{{f}} \\right) +D \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei}+\\frac{1}{2}\\, D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2}+\\frac{1}{6}\\, D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{3}+O \\left( {{\\it ei}}^{4} \\right) }{D \\left( f \\right) \\left( x_{{f}} \\right) + D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei} +\\frac{1}{2}\\, D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2}+O \\left( {{\\it ei}}^{3} \\right) }}\n", "$$\n", "```maple\n", "> ei2:=simplify(convert(ei1,polynom));\n", "```\n", "$$\n", "{\\it ei2}\\, := \\,\\frac{1}{3}\\,\\frac {3\\, D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2}+2\\, D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{3}\n", "-6\\,f \\left( x_{{f}} \\right) }{2\\,D \\left( f \\right) \\left( x_{{f}} \\right) +2\\, D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei}+ D^{ \\left( 3 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2}\n", "}\n", "$$\n", "```maple\n", "> ei3:=series(ei2,ei,3);\n", "```\n", "$$\n", "{\\it ei3}\\, := \\,-{\\frac {f \\left( x_{{f}} \\right) }{D \\left( f \\right) \\left( x_{{f}} \\right) }}+{\\frac {f \\left( x_{{f}} \\right) \\left( D^{ \\left( 2 \\right) } \\right) \\left( f \\right) \\left( x_{{f}} \\right) {\\it ei}}{ \\left( D \\left( f \\right) \\left( x_{{f}} \\right) \\right) ^{2}}}+ \\notag \\\\\n", "\\frac{1}{6}\\, \\frac{ 3\\, \\left( D^{ \\left( 2 \\right) } \\right) \\left( f \\right) \\left( x_{{f}} \\right) +3\\,{\\frac {f \\left( x_{{f}} \\right) \\left( D^{ \\left( 3 \\right) } \\right) \\left( f \\right) \\left( x_{{f}} \\right) }{D \\left( f \\right) \\left( x_{{f}} \\right) }}-6\\,{\\frac {f \\left( x_{{f}} \\right) \\left( \\left( D^{ \\left( 2 \\right) } \\right) \\left( f \\right) \\left( x_{{f}} \\right) \\right) ^{2}}{ \\left( D \\left( f \\right) \\left( x_{{f}} \\right) \\right) ^{2}}}}\n", "{ \\left( D \\left( f \\right) \\left( x_{{f}} \\right) \\right)}{{\\it ei}}^{2} +O \\left( {{\\it ei}}^{3} \\right)\n", "$$\n", "```maple\n", "> subs(f(x[f])=0,ei3);\n", "```\n", "$$\n", "\\frac{1}{2}\\,{\\frac { D^{ \\left( 2 \\right) } \\left( f \\right) \\left( x_{{f}} \\right) {{\\it ei}}^{2}}{D \\left( f \\right) \\left( x_{{f}} \\right) }}+O \\left( {{\\it ei}}^{3} \\right)\n", "$$\n", "注意すべきは，この収束性には一回の計算時間の差は入っていないことである．Newton法で解析的に微分が求まらない場合，数値的に求めるという手法がとられるが，これにかかる計算時間はばかにできない．二分法を改良した割線法(secant method)がより速い場合がある(NumRecipe9章参照)．\n", "\n", "二分法では，収束は遅いが，正負の関数値の間に連続関数では必ず解が存在するという意味で解が保証されている．しかし，Newton法では，収束は速いが，必ずしも素直に解に収束するとは限らない．解を確実に囲い込む，あるいは解に近い値を初期値に選ぶ手法が種々考案されている．解が安定であるかどうかは，問題，解法，初期値に大きく依存する．収束性と安定性のコントロールが数値計算のツボとなる．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 収束判定条件\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "どこまで値が解に近づけば計算を打ち切るかを決める条件を収束判定条件と呼ぶ．以下のような条件がある．\n", "\n", "\n", "|手法|判定条件|解説\n", "|:----|:----|:----|\n", "|$\\varepsilon$(イプシロン，epsilon)法 |\n", "|$\\delta$(デルタ，delta)法 |\n", "|占部法 | $\\left|f(x_{i+1})\\right| > \\left|f(x_i)\\right|$ | 数値計算の際の丸め誤差までも含めて判定する条件\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $\\epsilon, \\delta$を説明するための図 \n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VPX9xvH3ByI7FRCMsSBLpVaqBknIhhI2ZVFBARWKiCyNCLgVKQjYKohazE8t1KLIqoK4ICqKIksAIWSlynIUpQoKahUBbSpr+P7+INCEBkgyk9xZntc5c5g7c78zz/dcM48z985cc84hIiJyXCWvA4iISGBRMYiISBEqBhERKULFICIiRagYRESkCBWDiIgUoWIQEZEiVAwiIlKEikFERIqI8DpAWdSvX981adKkTGP/85//ULNmTf8GCnCac3jQnMODL3POzc3d7ZxrcKb1grIYmjRpQk5OTpnGrlq1inbt2vk3UIDTnMOD5hwefJmzme0oyXr6KElERIpQMYiISBF+KQYz62JmW81sm5mNKeb+35jZejM7aGb3lWasiIhULJ+LwcwqA08DXYEWQF8za3HSanuAu4DUMowVEZEK5I93DHHANufc5865Q8ACoEfhFZxz3znnsoHDpR0rIiIVyx9HJf0S+KrQ8k4g3t9jzSwFSAGIjIxk1apVpQ4KkJeXV+axwUpzDg+ac3ioiDkHzeGqzrnpwHSA2NhYV9bDtXR4W3jQnMOD5lw+/PFR0i6gUaHlhgW3lffYUsvIyODll19GpzMVETk1fxRDNtDczJqaWRWgD/BWBYwttRdeeIFnnnmG22+/ncOHT97dISIi4IePkpxzR8xsBLAUqAzMcs5tMbOhBfc/Y2bnATnAL4CjZnYP0MI591NxY33NdCpTp07lxx9/5LnnnmP79u28+uqrnH322eX1dCIiQckv+xicc0uAJSfd9kyh699y7GOiEo0tL5UqVWLIkCF07NiRlJQU2rRpwzvvvEPjxo0r4ulFRIJCWH7zeeDAgSxdupSdO3cSHx9Pdna215FERAJGWBYDQIcOHVi/fj3Vq1cnOTmZ119/3etIIiIBIWyLAeDiiy8mMzOT6OhoevfuTWpqqo5YEpGwF9bFAHDuueeycuVKevfuzahRoxg6dKiOWBKRsBY0X3ArT9WrV2fBggVceOGFPProo3zxxRc6YklEwlbYv2M4rlKlSjzyyCPMnDmTtLQ02rRpw44dJTqnhYhISFExnGTQoEEsXbqUXbt2ER8fT1ZWlteRREQqlIqhGB06dCA9PZ0aNWqQnJzMa6+95nUkEZEKo2I4hYsvvpiMjAwuv/xybrzxRv7yl7/oiCURCQsqhtM4fsTSzTffzJgxYxgyZIiOWBKRkKejks6gWrVqzJ8/n+bNm/Pwww+zfft2XnvtNerWret1NBGRcqF3DCVQqVIlJk6cyNy5c/nggw9ISkri888/9zqWiEi5UDGUwq233sry5cv57rvviI+PZ926dV5HEhHxOxVDKbVt25aMjAzq1q1Lhw4dmD9/vteRJIBs3LiRoUOHcskll1C7dm0qV65M7dq1ueyyyxg1ahTff/+91xFFzkjFUAbNmzcnIyODxMRE+vXrx0MPPaQjloS5c+cSExPDs88+y5YtW3DOUbduXQ4fPsymTZtITU1l8+bNXscUOSMVQxnVq1eP999/nwEDBvDggw/Sv39/Dh486HUs8ciPP/7I8OHDOXLkCH/84x/59ttvycvLY/fu3ezfv59//vOfTJo0iZYtW3odVeSMdFSSD6pUqcLs2bNp3rw548ePZ/v27SxatIgGDRp4HU0qWE5ODv/5z3+oVKkSY8eOLfI7W2ZGs2bNGDt2rIcJRUpO7xh8ZGaMGzeOl19+mdzcXBISEvj444+9jiUV7Fe/+hUREREcPXqUuLg4nnrqKTZu3KiPGCUoqRj85KabbmLVqlXk5eWRmJjIihUrvI4kFahJkybMnDkTOPZO8oknniA6Opr69eszbNgwvv32W48TipScisGP4uPjyczMpFGjRnTp0oXnnnvO60hSQTIzM5kwYQJTp05l06ZNfPnll3z99df07NmTadOmER0dzdatW72OKVIiKgY/a9KkCevWraNTp06kpKQwatQo8vPzvY4l5SgzM5OOHTsyYMAARowYceL2qKgonnvuObp37853333HXXfd5WFKkZJTMZSDX/ziFyxevJgRI0aQmppKr169yMvL8zqWlIN9+/bRu3dv6tWrx+jRo4tdZ/DgwQCsXLmS/fv3V2Q8kTJRMZSTiIgIpk6dypQpU1i8eDFt27Zl586dXscSP3v88cfZuXMnAwYMoEqVKsWuc/wotSNHjrBv376KjCdSJiqGcnbnnXeyePFiPvvsM+Lj48nNzfU6kvjRnDlzAOjYseMp19m9ezdw7Ag2/fiiBAMVQwXo1q0b6enpRERE0LZtWxYtWuR1JPGDL774gq+//hqASy+99JTrpaenA8f2P1WrVq1Cson4QsVQQS699FKysrK49NJL6dWrF5MnT9Yx7kGu8C/snu6dwJtvvglAp06dyj2TiD+oGCpQZGQkaWlp3HjjjYwePZrBgwdz6NAhr2NJGRU+2uyHH34odp3Fixef+MLjoEGDKiSXiK/8Ugxm1sXMtprZNjMbU8z9ZmZTCu7faGatCt233cw2mdmHZpbjjzyBrHr16rz00ks88MADzJ49m86dO7Nnzx6vY0kZNG3a9MT14+8KCtu9ezd33HEHANdeey0JCQkVlk3EFz4Xg5lVBp4GugItgL5m1uKk1boCzQsuKcC0k+5v75xr6ZyL9TVPMKhUqRITJkzgxRdfJD09nYSEBD799FOvY0kpNW/enOjoaABGjx7N+++/D4BzjuXLl9OmTRt27dpF48aNmT59updRRUrFH+8Y4oBtzrnPnXOHgAVAj5PW6QE8747JAOqYWZQfnjuo9evXj7S0NPbt20dCQgJpaWleR5JSmj59OjVq1GDPnj107tyZ2rVrU7NmTa666io+/fRTWrZsycqVK4mKCvv/3CWI+KMYfgl8VWh5Z8FtJV3HAcvNLNfMUvyQJ6gkJSWRmZlJVFQUV199NTNmzPA6kpRCXFwcGzZsoF+/fkRFRXHw4EFq1KhBx44dmTlzJtnZ2TRr1szrmCKlYr4eGWNmvYEuzrkhBcv9gXjn3IhC67wNPOacW1uwvAIY7ZzLMbNfOud2mdm5wDLgTufcmmKeJ4VjH0MRGRkZs2DBgjLlzcvLo1atWmUaW57y8vKYMGEC2dnZ3HTTTaSkpFC5cmW/PXYgzrk8ac7hQXMunfbt2+eW6CN755xPFyARWFpo+X7g/pPWeRboW2h5KxBVzGM9CNx3pueMiYlxZZWWllbmseXt8OHDbsSIEQ5w1113nfv3v//tl8cN5DmXF805PGjOpQPkuBK8rvvjo6RsoLmZNTWzKkAf4K2T1nkLuLXg6KQE4Efn3DdmVtPMagOYWU3gaiBsz314/Gc0/va3v/HOO+9wxRVX8NVXX515oIiIH/lcDM65I8AIYCnwMfCKc26LmQ01s6EFqy0BPge2Ac8BwwpujwTWmtlHQBbwjnPuPV8zBbvhw4fzzjvv8MUXXxAXF0dWVpbXkUQkjPjl1J7OuSUce/EvfNszha47YHgx4z4Hov2RIdR06dKF9evXc+2115KcnMycOXO4+eabvY4lImFA33wOYC1atCAzM5OYmBj69OnDhAkT9DMaIlLuVAwBrkGDBqxYsYL+/fvz5z//mX79+nHgwAGvY4lICFMxBIGqVasyd+5cHnnkEV566SXat2+vcwiLSLlRMQQJM+P+++9n4cKFbNy4kbi4ODZu3Oh1LBEJQSqGINOzZ08++OADjh49SlJSEosXL/Y6koiEGBVDEGrVqhVZWVlcfPHF9OjRg9TUVO2UFhG/UTEEqfPPP5/Vq1fTu3dvRo0axZAhQ3RuBxHxC798j0G8UaNGDRYsWMBvfvMbJk6cyD//+U8WLlzIOeec43U0EQliescQ5I6f22H+/PlkZGQQHx9/4oxhIiJloWIIEX379mXVqlXk5eWRkJDA0qVLvY4kIkFKxRBCEhISyMrKomnTpnTr1o33br0V17gxyR06QJMmMG+e1xFFJAhoH0OIueCCC1i7di3PJidz5QsvYMfv2LEDUgrOg9Svn1fxRCQIqBhCUK1atfjD7t0YcA9P8iEtj93xMzC42rHftw1x+/a1pE4dr1NUrHCb87ZtUL16DJ995nWS0KNiCFF2qvM4HNTvLEloyMuDgwereB0jJKkYQtUFF8COHTzFvUVu3n/uuVRf9S+PQlWcVas+pF27dl7HqFDhNud27WDfvv1AVa+jhBztfA5VkyZBjRpFbtpvxu+//56nn37ao1AiEgxUDKGqXz+YPh0aN8aZHft3+nR+vOYaRowYwfDhwzl8+LDXKUUkAKkYQlm/frB9O6tXroTt26kxZAhvvPEGo0aN4u9//ztdu3Zl7969XqcUkQCjYggzlStXZvLkycyePZs1a9aQkJDAp59+6nUsEQkgKoYwddttt7Fy5Ur27NlDfHw8y5cv9zqSiAQIFUMYu+KKK8jOzqZhw4Z06dKFv//9715HEpEAoGIIc02aNCE9PZ1u3boxfPhwhg0bpp3SImFOxSDUrl2bRYsWMWrUKKZNm0bXrl3Zs2eP17FExCMqBgGK7pT+4IMPiI+P55NPPvE6loh4QMUgRRzfKf3TTz+RkJDA+++/73UkEalgKgb5H23atCErK4vGjRvTtWtXpkyZonNKi4QRFYMUq3Hjxqxbt47u3btz9913c/vtt+uc0iJhwi/FYGZdzGyrmW0zszHF3G9mNqXg/o1m1qqkY8U7tWrVYuHChYwbN47nnnuOq666it27d3sdS0TKmc/FYGaVgaeBrkALoK+ZtThpta5A84JLCjCtFGPFQ5UqVeLhhx9m3rx5ZGZmEhcXx+bNm72OJSLlyB/vGOKAbc65z51zh4AFQI+T1ukBPO+OyQDqmFlUCcdKAPjd737HmjVrOHDgAImJiSxevNjrSCJSTvxxPoZfAoXPCrMTiC/BOr8s4Vi/ueeee1i1ahV1wuk0V8C+ffv8NufGjRuzefNmunfvTrNmzWjUqJFfHtff/DnnYBFuc/7ww6c4cuQI7drd53WUClW/fv1yP+9G0Jyox8xSOPYxFJGRkaxatarUj7Fz507y8/PZt2+fn9MFNn/PuWnTpnz11Vd8/vnn7N27l0aNGmFmZx5YgbSdQ9+RI0dwzoXVnAHOPvvsMr3+lYpzzqcLkAgsLbR8P3D/Ses8C/QttLwViCrJ2OIuMTExrqzS0tLKPDZYlcecjx496iZMmOAAFx8f77755hu/P4cvtJ1DX3Kyc9HRe72OUeF82c5AjivB67o/9jFkA83NrKmZVQH6AG+dtM5bwK0FRyclAD86574p4VgJQGbGAw88wMKFC9m0aROtW7dmw4YNXscSET/wuRicc0eAEcBS4GPgFefcFjMbamZDC1ZbAnwObAOeA4adbqyvmaTi9OzZk/T0dCpVqsQVV1zBK6+84nUkEfGRX/YxOOeWcOzFv/BtzxS67oDhJR0rwSU6Oprs7Gx69uzJzTffzObNm3nwwQepVEnfnxQJRvrLFb8499xzWbFiBQMHDmTixInceOON5OXleR1LRMpAxSB+U7VqVWbOnMkTTzzBG2+8QZs2bdixY4fXsUSklFQM4ldmxr333suSJUvYsWMHrVu3Zu3atV7HEpFSUDFIuejcuTOZmZnUrVuXDh06MHPmTK8jiUgJqRik3Fx00UVkZGTQvn17hgwZwt13382RI0e8jiUiZ6BikHJVt25d3nnnHf7whz8wZcoUnTZUJAioGKTcRURE8H//93/MmjWLNWvWEB8fz8cff+x1LBE5BRWDVJiBAweSlpZ24rShS5bo6ysigUjFIBUqKSmJ7OxsfvWrX3HttdcyefJknTZUJMCoGKTCXXDBBaxdu5Ybb7yR0aNH079/f/bv3+91LBEpoGIQT9SoUYMFCxYwceJE5s2bR3JyMrt27fI6loigYhAPmRnjx49n0aJFfPzxx7Ru3ZrMzEyvY4mEPRWDeO76669n/fr1VKtWjeTkZF544QWvI4mENRWDBIRLLrmErKwskpKSuPXWWxk1ahT5+flexxIJSyoGCRj169dn6dKlDB8+nNTUVK699tqwO22jSCBQMUhAOeuss/jb3/7Gs88+y/Lly4mPj2fr1q1exxIJKyoGCUgpKSmsWLGCvXv3EhcXx7vvvut1JJGwoWKQgNW2bVuys7Np2rQp11xzDY8//ri+DCdSAVQMEtAaN27MunXr6N27N3/84x+59dZb9WU4kXKmYpCAV7NmTV5++WUefvhhXnzxRdq2bcvOnTu9jiUSslQMEhTMjHHjxvHmm2/yySef0Lp1a9avX+91LJGQpGKQoNK9e3cyMjKoWbMm7dq1Y9asWV5HEgk5KgYJOr/97W/JysoiOTmZwYMHc9ddd3H48GGvY4mEDBWDBKV69eqxZMkS7r33XqZOnUrnzp3ZvXu317FEQoKKQYJWREQETzzxBHPmzCE9PZ24uDg2bdrkdSyRoKdikKA3YMAAVq9ezYEDB0hMTOT111/3OpJIUFMxSEiIj48nJyeHSy65hF69ejFr1iyOHj3qdSyRoKRikJBx/vnns2rVKgYOHMgLL7xAz549+emnn7yOJRJ0fCoGM6tnZsvM7LOCf+ueYr0uZrbVzLaZ2ZhCtz9oZrvM7MOCSzdf8ohUq1aNmTNnMmLECN5++20SExPZtm2b17FEgoqv7xjGACucc82BFQXLRZhZZeBpoCvQAuhrZi0KrfKkc65lwWWJj3lEMDN69erF+++/z7fffkvr1q15//33vY4lEjR8LYYewNyC63OB64tZJw7Y5pz73Dl3CFhQME6kXHXo0IHs7GwaNWpE165dSU1N1Y/wiZSA+fKHYmb7nHN1Cq4bsPf4cqF1egNdnHNDCpb7A/HOuRFm9iAwEPgRyAFGOuf2nuK5UoAUgMjIyJgFCxaUKXNeXh61atUq09hgFe5z3r9/P4899hhr1qyhU6dO3HfffVStWtXjhP4Xbtv5nntakp+fz9Sp4XWIsi/buX379rnOudgzruicO+0FWA5sLubSA9h30rp7ixnfG5hRaLk/8LeC65FAZY69c5kEzDpTHuccMTExrqzS0tLKPDZYac7OHT161E2cONGZmYuJiXFffvmlN8HKUbht5+Rk56Kj9/r8ONu3b3djx451LVu2dOecc46rWrWqa9iwoUtOTnYjR450P//8s+9h/ciX7QzkuBK8xkaUoDg6neo+M/uXmUU5574xsyjgu2JW2wU0KrTcsOA2nHP/KvRYzwFvnymPSFmYGePHjyc6Opp+/foRGxvLa6+9xpVXXul1NPHQ66+/zi233HLip9xr1qxJtWrV2LVrFzt37uSjjz4iNTXV45QVz9d9DG8BAwquDwDeLGadbKC5mTU1sypAn4JxFJTJcTdw7J2ISLm57rrryMrKok6dOnTo0IFnnnnG60jikV27dnHLLbdw4MABxo4dy5dffkleXh779u3j0KFDbNiwgTlz5ngd0xO+FsNjwFVm9hnQqWAZMzvfzJYAOOeOACOApcDHwCvOuS0F4yeb2SYz2wi0B+71MY/IGf3mN78hMzOTq6++mjvuuIOUlBQOHTrkdSypYG+//Tb79++ne/fuTJo0iUaN/vvBRkREBJdffjk9eoTncTJn/CjpdJxzPwAdi7n9a6BboeUlwP8ciuqc6+/L84uUVZ06dXjrrbd44IEHePTRR9myZQsLFy7kvPPO8zqaVJAjR44AkJmZyZo1a2jbtq3HiQKHvvksYaty5co88sgjvPzyy3z44YfExsaSlZXldSypIH369OHXv/413377LcnJyVSrVo3zzjuPCy+80OtonlMxSNi76aabSE9P56yzzqJt27bMnTv3zIMk6J1zzjm88sorREdHA3Dw4EH+9a9/Ua9ePY+TeU/FIAJER0eTnZ1NmzZtuO2227j77rt18p8Q5pxj9OjRxMTEcOmll5KVlcVPP/2Ec07vGlExiJxQv359li5dyj333MOUKVN08p8QlpqayuTJk7njjjt44YUXaN26NbVr1/Y6VsBQMYgUEhERwZNPPsnzzz9Peno6sbGx/OMf//A6lvjZX//6VwCGDh3qcZLApGIQKUb//v1Zu3Yt+fn5tGnThpdeesnrSOInBw4cYNeuXQB8//33HqcJTCoGkVOIjY0lJyeH2NhYfve73zFq1KgThzhK8KpWrRpRUce+W3v77bezbNmyE/uTDh48yCeffMKkSZN49dVXvYzpKRWDyGlERkayfPlyhg0bRmpqKl27duWHH37wOpb46OGHHwbg008/5eqrr6ZatWrUq1eP6tWrc/HFFzN+/HjOPvtsj1N6R8UgcgZVqlTh6aefZsaMGaxZs4bWrVvz0UcfeR1LfDBo0CCWLVtGz549adiwIRERERw8eJDGjRtzzTXX8Pjjj5OUlOR1TM/49M1nkXAyePBgfvvb39KrVy+SkpKYPXs2N910k9expIw6depEp06n/I3QsKZ3DCKlkJCQQE5ODi1btuTmm29mzJgx5Ofnex1LxK9UDCKlFBUVRVpaGrfffjt/+ctfuOaaa9izZ4/XsUT8RsUgUgZVqlThmWeeYfr06axcuZLWrVuzaVN4nUlMQpeKQcQHv//971m9ejX79+8nMTGR1157zetIIj5TMYj4KDExkZycHC677DJuvPFG7r//fu13kKCmYhDxg/PPP5+0tDRSUlJ47LHHtN9BgpqKQcRPqlatyrPPPsuzzz6r/Q4S1FQMIn6WkpJyYr9DQkICr7zyiteRREpFxSBSDhITE8nNzT3xfYfRo0drv4MEDRWDSDk5/n2HoUOHMnnyZP3OkgQNFYNIOapSpQrTpk1jxowZrF69mtjYWP3OkgQ8FYNIBRg8eDBr1qzh8OHDJCYm6vwOEtBUDCIVJD4+ntzc3BPndxg5cqTO7yABScUgUoEiIyNZsWIFI0aM4IknnqBz5846i5gEHBWDSAU766yzmDp1KnPmzGHdunXExsaSm5vrdSyRE1QMIh4ZMGAA69atA6BNmzbMnTvX40Qix6gYRDwUExNDTk4Obdq04bbbbmPEiBEcOnTI61gS5nwqBjOrZ2bLzOyzgn/rnmK9WWb2nZltLst4kVDWoEEDli5dyn333cfTTz9Nx44d+eabb7yOJWHM13cMY4AVzrnmwIqC5eLMAbr4MF4kpEVERPD444/z0ksvsWHDBmJiYli/fr3XsSRM+VoMPYDjH4zOBa4vbiXn3BqguJ+aLNF4kXDRp08f1q9fT/Xq1UlOTmbatGk457yOJWHG12KIdM4df8/7LRBZweNFQs5ll11GTk4OnTp1YtiwYQwePJgDBw54HUvCSMSZVjCz5cB5xdw1rvCCc86ZWZn/1+ZM480sBUiBY8eCr1q1qkzPk5eXV+axwUpzDk4jR46kfv36zJ49m/T0dB566CEiI0/9/06hMOfS2LevJfn5+WE1Z6ig7eycK/MF2ApEFVyPAraeZt0mwOayji98iYmJcWWVlpZW5rHBSnMObm+88YarXbu2q1+/vluxYsUp1wulOZdEcrJz0dF7vY5R4XzZzkCOK8FrrK8fJb0FDCi4PgB4s4LHi4S8Hj16kJ2dTYMGDbjqqqtITU3VfgcpV74Ww2PAVWb2GdCpYBkzO9/MlhxfycxeAtYDF5nZTjMbfLrxIlLURRddRGZmJjfccAOjRo2iT58+5OXleR1LQtQZ9zGcjnPuB6BjMbd/DXQrtNy3NONF5H/Vrl2bV199lcmTJzN27Fi2bNnCokWLaN68udfRJMTom88iQcTMGD16NO+99x7ffPMNsbGxLF682OtYEmJUDCJB6KqrriI3N5cLL7yQ7t278+c//5mjR496HUtChIpBJEg1adKEtWvXcttttzFhwgTGjh3L3r17vY4lIUDFIBLEqlevzqxZs5g2bdqJkwBt3LjR61gS5FQMIkHOzBg6dCh//etfOXDgAAkJCcybN8/rWBLEVAwiIaJFixZs2LCB1q1bc8stt3DXXXdx+PBhr2NJEFIxiISQyMhIli9fzj333MPUqVPp0KGDfsJbSk3FIBJizjrrLJ588knmz59/4ie8j58pTqQkVAwiIapv375kZGRQs2ZN2rVrx9SpU/VTGlIiKgaREHbppZeSnZ1N165dueuuu+jfvz8///yz17EkwKkYREJcnTp1eOONN5g4cSLz588nMTGRbdu2eR1LApiKQSQMVKpUifHjx7NkyRK++uorYmNjefvtt72OJQFKxSASRrp06UJubi7NmjXjuuuu409/+hP5+flex5IAo2IQCTNNmzZl3bp13HbbbUycOJFrrrmGH374wetYEkBUDCJhqPBPaaxcuZLY2Fg2bNjgdSwJECoGkTB1/Kc0PvjgA44cOUJSUhKzZs3yOpYEABWDSJiLj49nw4YNXHHFFQwePJjbb7+dgwcPeh1LPKRiEBEaNGjAe++9x5gxY5g+fTpXXnklX375pdexxCMqBhEBICIigkcffZRFixaxdetWWrVqxbJly7yOJR5QMYhIEddffz3Z2dlERUXRuXNnJk2apLPDhRkVg4j8j1//+tdkZGTQt29fxo8fz/XXX8++ffu8jiUVRMUgIsWqWbMmL774IlOmTOHdd98lNjaWjz76yOtYUgFUDCJySmbGnXfeyerVq9m/fz+JiYk8//zzXseScqZiEJEzSkpKYsOGDcTHxzNgwACGDRumQ1pDmIpBREokMjKSZcuWMWrUKKZNm0bbtm356quvvI4l5UDFICIlFhERweTJk1m4cCEff/wxrVq1Yvny5V7HEj9TMYhIqfXs2ZPs7GwiIyO5+uqrdUhriPGpGMysnpktM7PPCv6te4r1ZpnZd2a2+aTbHzSzXWb2YcGlmy95RKTiXHTRRWRmZtKnTx/Gjx9Pjx492Lt3r9exxA98fccwBljhnGsOrChYLs4coMsp7nvSOdey4LLExzwiUoFq1qzJvHnzmDp1KkuXLiWqMte9AAAHb0lEQVQmJoZ//OMfXscSH/laDD2AuQXX5wLXF7eSc24NsMfH5xKRAGRmjBgxgtWrV3Po0CGSkpKYPXu217HEB74WQ6Rz7puC698CkWV4jDvNbGPBx03FfhQlIoEvMTGRDRs20KZNGwYNGsTvf/97Dhw44HUsKQNzzp1+BbPlwHnF3DUOmOucq1No3b3OuVPtZ2gCvO2cu6TQbZHAbsABE4Eo59ygU4xPAVIAIiMjYxYsWHDa3KeSl5dHrVq1yjQ2WGnO4SFQ5pyfn8/s2bOZN28ezZs356GHHiIqKsrvz3PPPS3Jz89n6tRNfn/sQObLdm7fvn2ucy72jCs658p8AbZy7MUcIArYepp1mwCby3p/4UtMTIwrq7S0tDKPDVaac3gItDkvXrzY1alTx9WpU8ctXrzY74+fnOxcdPRevz9uoPNlOwM5rgSvsb5+lPQWMKDg+gDgzdIMNrPC/xtxA7D5VOuKSHC59tpryc3NpWnTplx33XWMGzeO/Px8r2NJCfhaDI8BV5nZZ0CngmXM7HwzO3GEkZm9BKwHLjKznWY2uOCuyWa2ycw2Au2Be33MIyIBpFmzZqxbt47BgwfzyCOP0LlzZ7777juvY8kZRPgy2Dn3A9CxmNu/BroVWu57ivH9fXl+EQl81atXZ8aMGSQlJTF8+HBatWrFq6++SmJiotfR5BT0zWcRqRCDBg0iPT2dqlWr0rZtW6ZMmXJ8/6IEGBWDiFSYyy+/nNzcXLp168bdd99Nnz59+Pe//+11LDmJikFEKlSdOnVYtGgRjz32GK+99hpxcXFs2bLF61hSiIpBRCpcpUqVGD16NCtWrGDv3r3ExcUxb948r2NJARWDiHimXbt2bNiwgZiYGG655RadAChAqBhExFPnn38+K1as4L777mPatGlceeWV7Nixw+tYYU3FICKeO+uss3j88cd5/fXX2bp1K61ateLdd9/1OlbYUjGISMC44YYbyMnJoWHDhnTr1o0HHnhA35b2gIpBRAJK8+bNycjIYODAgTz88MP6trQHVAwiEnCqV6/OrFmzmDlzJuvWrePyyy9n3bp1XscKGyoGEQlYgwYNYv369VSvXp127drxxBNP6NvSFUDFICIBrWXLluTm5nLdddcxcuRIevfuzc8zZkBGBnU++hCaNAF9B8KvVAwiEvDOPvtsFi5cSGpqKtUXLYKUFDhYcHa4HTuOLasc/EbFICJBwcwYOXIkM849lxonf5z0888wbpw3wUKQTz+7LSJS0aoVHKHUkg+L3vHllx6kCU0qBhEJLhdcADt28NTJ5/W64AJv8oQgfZQkIsFl0iSoUaPobTVqHLtd/ELFICLBpV8/mD4dGjfGmUHjxseW+/XzOlnIUDGISPDp1w+2b2f1ypWwfbtKwc9UDCIiUoSKQUREilAxiIhIESoGEREpQsUgIiJFWDD+UqGZfQ+U9dx/9YHdfowTDDTn8KA5hwdf5tzYOdfgTCsFZTH4wsxynHOxXueoSJpzeNCcw0NFzFkfJYmISBEqBhERKSIci2G61wE8oDmHB805PJT7nMNuH4OIiJxeOL5jEBGR0wipYjCzLma21cy2mdmYYu7vZ2YbzWyTmaWbWXRJxwYiH+e7veD2D80sp2KTl10J5tyjYM4fmlmOmV1R0rGBysc5h+R2LrReazM7Yma9Szs20Pg4Z/9uZ+dcSFyAysA/gWZAFeAjoMVJ6yQBdQuudwUySzo20C6+zLdgeTtQ3+t5lMOca/Hfj0gvAz4J1m3s65xDeTsXWm8lsAToHerb+VRzLo/tHErvGOKAbc65z51zh4AFQI/CKzjn0p1zewsWM4CGJR0bgHyZb7AqyZzzXMFfClATcCUdG6B8mXOwKum2uhNYCHxXhrGBxpc5+10oFcMvga8KLe8suO1UBgPvlnFsIPBlvnDsxWO5meWaWUo55CsPJZqzmd1gZp8A7wCDSjM2APkyZwjR7WxmvwRuAKaVdmyA8mXO4OftHJbnfDaz9hx7obziTOuGglPM9wrn3C4zOxdYZmafOOfWeJPQv5xzi4BFZtYWmAh08jhSuTvNnEN1Oz8FjHbOHTUzr7NUlNPN2a/bOZSKYRfQqNByw4LbijCzy4AZQFfn3A+lGRtgfJkvzrldBf9+Z2aLOPZWNtBfMEq1nZxza8ysmZnVL+3YAFLmOTvndofwdo4FFhS8QNYHupnZkRKODURlnrNz7g2/b2evd7r4cedNBPA50JT/7rz57UnrXABsA5JKOzbQLj7OtyZQu9D1dKCL13Py05wv5L87YlsV/HFZMG5jP8w5ZLfzSevP4b87n0N2O59mzn7fziHzjsE5d8TMRgBLObbnfpZzbouZDS24/xngT8A5wN8LWveIcy72VGM9mUgJ+TJfIJJjHzvAsf8g5zvn3vNgGqVSwjn3Am41s8PAfuBmd+wvJui2Mfg2ZzML5e1cqrEVkdsXvsyZcvh71jefRUSkiFA6KklERPxAxSAiIkWoGEREpAgVg4iIFKFiEBGRIlQMIiJShIpBRESKUDGIiEgR/w+FuF7VLurYMgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "def func(x):\n", " return 0.4*(x**2-4*x+1)\n", "x1=0.25\n", "x0=0.4\n", "x = np.linspace(0.2, 0.4, 100)\n", "y = func(x)\n", "plt.plot(x, y, color = 'k')\n", "plt.plot(x1, func(x1), \"o\", color = 'r')\n", "plt.plot(x0, func(x0), \"o\", color = 'r')\n", "plt.plot([0.2,0.45],[0,0], color = 'k')\n", "plt.plot([x1,x0],[func(x1),func(x1)], color = 'b')\n", "plt.plot([x0,x0],[func(x0),func(x1)], color = 'b')\n", "\n", "plt.text(0.41, -0.07, r'$\\epsilon$', size='24') \n", "plt.text(0.32, 0.05, r'$\\delta$', size='24') \n", "\n", "\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2変数関数の場合\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "2変数の関数では，解を求める一般的な手法は無い．この様子は実際に2変数の関数で構成される面の様子をみれば納得されよう．\n", "```maple\n", "> restart;\n", "> f:=(x,y)->4*x+2*y-6*x*y; g:=(x,y)->10*x-2*y+1;\n", "```\n", "$$\n", "f\\, := \\,( {x,y} )\\mapsto 4\\,x+2\\,y-6\\,xy \\notag \\\\\n", "g\\, := \\,( {x,y} )\\mapsto 10\\,x-2\\,y+1 \\notag\n", "$$\n", "```maple\n", "> p1:=plot3d({f(x,y)},x=-2..2,y=-2..2,color=red):\n", " p2:=plot3d({g(x,y)},x=-2..2,y=-2..2,color=blue):\n", " p3:=plot3d({0},x=-2..2,y=-2..2,color=gray):\n", " with(plots):\n", " display([p1,p2,p3],axes=boxed,orientation=[-150,70]);\n", "```\n", "\n", "\n", "解のある程度近くからは，Newton法で効率良く求められる．\n", "```maple\n", "> fsolve({f(x,y)=0,g(x,y)=0},{x,y});\n", "```\n", "$$\n", "\\left\\{ x=- 0.07540291160,y= 0.1229854420 \\right\\}\n", "$$\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": true }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('