{ "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-21\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": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiNUlEQVR4nO3dfZzVc/7/8ce7qbQ1pa00pLZIIlFMpFo0ftbXRSy5FrtYRpuLkkQuU0p0gZVFIu0qQ3IZKmpysYRKKLlsk1IskaZsqXn//niNLTXVdM7nzPt8znneb7dz05mZ8/m8PjPjOe/z/rwvnPceERGJryqhCxARkeQoyEVEYk5BLiIScwpyEZGYU5CLiMRc1RAnbdCggW/WrFlCr129ejW1atWKtqBAdC3pJ1OuA3Qt6SqZa5k9e/a33vtdNv94kCBv1qwZs2bNSui1M2bMoHPnztEWFIiuJf1kynWAriVdJXMtzrkvyvu4ulZERGJOQS4iEnMKchGRmFOQi4jEnIJcRCTmIgly51xd59wTzrmPnHMLnHMdojiuiIhsX1TDD+8CJnvvT3XOVQdqRnRcERHZjqRb5M65nYHDgQcBvPfrvPc/JHvc8kybBuPH/y4VhxYRSak1a6BXL/jqqxqRH9slux65c64tMAr4EGgDzAZ6eu9Xb/Z1hUAhQF5eXn5RUdEOn+u++/ZkwoTGPProWzRsuDapukMrKCgAoLi4OHAl0SgpKSE3Nzd0GUnLlOsAXUu6ee653RgxoiVDhrxB+/brEjpGQUHBbO99uy0+4b1P6gG0A9YD7cue3wUM3NZr8vPzfSIWLfK+SpVSf801Cb08rQDevv2Zobi4OHQJkciU6/Be15JOSku9328/7w880Pvp04sTPg4wy5eTqVHc7FwCLPHev1X2/AngoAiOu4WmTaFTp28ZNQp++ikVZxARid706TB/Plx+OTgX/fGTDnLv/XLgS+dcy7IP/T+smyUlunZdwooVMG5cqs4gIhKtu+6CXXaBM89MzfGjGkd+GTDOOfc+0BYYHNFxt9CmzUratLFvjLYbFZF09/nnMGkSdO8ONaK/zwlEFOTe+7ne+3be+wO89yd577+P4rjlcQ569oR582DGjFSdRUQkGiNHQk6OBXmqxHJm51lnQYMG1ioXEUlXq1bBQw/B6adDo0apO08sg7xGDbj4Ynj2WXvbIiKSjsaMgR9/tJucqRTLIAfo0QOqVoW//S10JSIiW9qwwXoNOnSA9u1Te67YBnmjRnDGGfa2ZeXK0NWIiPzapEmwcCFccUXqzxXbIAf7BpWUwOjRoSsREfm1O+6wuS8nn5z6c8U6yA86CA4/3LpX1q8PXY2IiHn3XXjlFbjsMusCTrVYBzlYq3zxYnjqqdCViIiYO+6A3Fy48MLKOV/sg/yEE2DPPe0bJyIS2rJlUFQEF1wAO+9cOeeMfZDn5NgEoTffhJkzQ1cjItlu5Ejr6k31kMNNxT7Iwf7y1a0Lw4eHrkREstnq1XDvvXaDs3nzyjtvRgR5bq5NEHrySRvuIyISwpgx8P33cOWVlXvejAhysLvDVarAnXeGrkREstGGDXav7tBDoWPHyj13xgT57rvD2WfbBKHvU7Zkl4hI+Z55xnoE+vSp/HNnTJAD9O5tfVT33x+6EhHJNsOG2Qi6k06q/HNnVJC3aQNHHWUThNYltiWeiMgOe/NNe/TqZSPpKltGBTnY25ply2D8+NCViEi2GDoUfvtbOP/8MOfPuCA/+mg44AD7xpaWhq5GRDLdJ5/A00/biqy5uWFqyLggdw6uugo+/BBefDF0NSKS6YYPh+rVbeRcKBkX5GDL2zZpArffHroSEclky5fD2LFw3nmQlxeujowM8mrVbATLq69q2r6IpM7dd9vAisqeALS5jAxysFXHfvtb6ysXEYnaqlXw979D167QokXYWjI2yHNz7ebDU0/ZzQgRkSiNHg0//GD35ELL2CAHu/mw005qlYtItNatgxEj4IgjUr8fZ0VkdJDn5dnKiGPHwtKloasRkUwxbhwsWQL9+oWuxGR0kINNECot1cYTIhKNDRvgttvgwANt3ko6yPgg32MPOPNMuO8+WLEidDUiEnfPPAMffwzXXGPzVtJBxgc5wNVX22Ja99wTuhIRiTPvYcgQ2GsvOOWU0NVslBVBvv/+0KUL3HWXBbqISCKmT4d33oG+fcMsjrU1WRHkYG+DvvsOHnggdCUiEleDB8Nuu8Gf/hS6kl/LmiDv1MmGCg0dCmvXhq5GROJm5kxrkffpY8Oa00nWBDnAddfBV1/ZcEQRkR0xaBDUr2/7A6ebyILcOZfjnHvXOTcpqmNG7aij4OCDbejQ+vWhqxGRuJg7FyZNso0jatUKXc2WomyR9wQWRHi8yDlnrfKFC6GoKHQ1IhIXgwdDnTpw6aWhKylfJEHunGsMHA+MjuJ4qXTCCTaKZfBgbTwhItv30UfwxBNwySVQt27oasrnvPfJH8S5J4BbgdpAH+99l3K+phAoBMjLy8svSrBJXFJSQm6S23BMn96QgQNb0b//PI444tukjpWogoICAIqLi4OcP2pR/FzSQaZcB+haonLrrfvwyiu7UFQ0k7p1f076eMlcS0FBwWzvfbstPuG9T+oBdAH+XvbvzsCk7b0mPz/fJ6q4uDjh1/5i/Xrv997b+zZtvC8tTfpwCQG8ffszQxQ/l3SQKdfhva4lCp995n1Ojve9ekV3zGSuBZjly8nUKLpWOgEnOucWAUXAkc65RyI4bsrk5MD118N778Fzz4WuRkTS1a23QtWqNgEonSUd5N77ft77xt77ZsCZwHTv/TlJV5ZiZ50FzZvDgAE27VZEZFOLFtlQ5cJCmwSUzrJqHPmmqlaFa6+F2bO1SbOIbGnIEKhSJf1b4xBxkHvvZ/hybnSmq3PPhaZN1SoXkV/78kt46CH4y1+gcePQ1Wxf1rbIwTZpvvZaeOstmDo1dDUiki6GDLH/XnNN2DoqKquDHOC886BJE+jfX61yEbHW+OjRlg2/+13oaiom64O8enWb7TlzJkyZEroaEQnt1lutUXfddaErqbisD3KA88+3vvKbblKrXCSbLV5srfELLrBMiAsFORtb5W+/rREsItls8GD777XXhq1jRynIy5x3HjRrpla5SLb64gsbqXLhhfHpG/+FgrxMtWpwww0waxY8/3zoakSksg0aZCukxq01DgryXzn3XJvteeONWhlRJJt8/jmMGWOzOOMwbnxzCvJNVKtmXSvvvgtPPhm6GhGpLAMGbJxXEkcK8s2cfTbsu6+1yjdsCF2NiKTaggXwyCO2aUS6r6myNQryzeTk2F/nBQtg/PjQ1YhIqt10E9SsGY81VbZGQV6Orl2hbVub7flz8uvIi0iamjsXJkyAK66ABg1CV5M4BXk5qlSBgQNtb88xY0JXIyKpcuONtn1b796hK0mOgnwrjj8eDj3Uull++il0NSIStTfesI1lrroqfffirCgF+VY4Z2suLF0K99wTuhoRiZL30K8f5OVBz56hq0megnwbOneGo4+2QF+5MnQ1IhKVKVPg1VdtEmCtWqGrSZ6CfDsGD4YVK2DYsNCViEgUSkutNb7HHnDRRaGriYaCfDvy8+G00+COO+Drr0NXIyLJmjDBRqsMGGAL5mUCBXkFDBwI//0v3HJL6EpEJBk//wzXXw/7728bsGcKBXkFtGxp6xPff7+tySAi8fTAA/DZZ9ZlmpMTuproKMgrqH9/qFo1XruGiMhGq1bBzTfD4Yfb8OJMoiCvoEaNbNLAY4/BO++ErkZEdtTw4fDNN3D77Ta8OJMoyHdA3742jffqq7X5hEicLF9uI89OPRXatw9dTfQU5DugTh0bd1pcDJMnh65GRCrq5pth7dqNW7llGgX5DureHfbc01rnWuZWJP199JHd5CwshBYtQleTGgryHVS9us30nDcPHn44dDUisj1XX23L1N50U+hKUkdBnoDTTrMFta6/HkpKQlcjIlszYwY8+6zN5GzYMHQ1qaMgT4Bzdgd8+XIYOjR0NSJSntJSuPJKaNIEevUKXU1qKcgT1LGjtcyHDrUVEkUkvYwbB3Pm2A3O3/wmdDWppSBPwpAhsH69jWQRkfSxZo1tpJyfb/vwZrqkg9w518Q5V+yc+9A5N985lwGr+1bMnnvC5ZfbTc85c0JXIyK/GD4cliyx/1bJguZqFJe4HrjSe98KOBS4xDnXKoLjxsL110P9+tYHp0lCIuEtXWrvlk85BY44InQ1lSPpIPfeL/Pezyn79ypgAbB7sseNi7p1bVXE116DiRNDVyMi/fpZl2c2DURwPsJmpHOuGfAq0Np7/+NmnysECgHy8vLyi4qKEjpHSUkJubm5SVYarQ0boLCwHWvW5DB27DtUr1663dcUFBQAUFxcnOryKkU6/lwSkSnXAdl5LQsW1KZHj3zOPvsLLrro35VQ2Y5L5udSUFAw23vfbotPeO8jeQC5wGyg6/a+Nj8/3yequLg44dem0vTp3oP3gwdX7OsBb9/+zJCuP5cdlSnX4X32XUtpqfcdOnifl+f9jz+mvqZEJfNzAWb5cjI1ktsAzrlqwERgnPf+ySiOGTcFBXDSSTBokIYjioQwfjy8+aYNN6xdO3Q1lSuKUSsOeBBY4L0fkXxJ8TV8uPXNXX116EpEssuqVXDVVdCuHZx3XuhqKl8ULfJOwLnAkc65uWWP4yI4buzsuaf9Mo0bB6+/HroakewxaBAsWwYjR2bHcMPNRTFq5XXvvfPeH+C9b1v2eCGK4uKoXz+bEnzZZVodUaQyfPopjBhhLfFMXGu8IrLwb1dq1axpC9jPnWtLZ4pIavXqZVPwhwwJXUk4CvIUOO006NzZ9vf89tvQ1YhkrueegxdesCVq8/JCVxOOgjwFnIO774aVK229BxGJ3po1tkRGq1bWlZnNFOQp0ro19OwJo0fDW2+FrkYk89x6KyxaBPfcA9Wqha4mLAV5CvXvD7vtBj166ManSJQ+/RRuvx26dbNuzGynIE+h2rXtbvqcOXDffaGrEckM3sOll0KNGjawQBTkKXf66XDUUXbj8+uvQ1cjEn8TJ8LUqTBwIOy6a+hq0oOCPMWcs0kKP/0EvXuHrkYk3lautBucBx5oXZZiFOSVoGVLmyg0fry1JEQkMddfb3vl3n8/VK0aupr0oSCvJNdcAy1aWCvip59CVyMSP++8YyNULrkEDj44dDXpRUFeSWrUsBuen39u60KISMVt2OAoLLRRYLfcErqa9KM3J5XoyCPh3HNt2JSIVNzEibszdy5MmAA77xy6mvSjFnklGz4c6tQJXYVIfCxcCA89tAddutg+nLIlBXkl22UXuPPO0FWIxIP30L075OR47r3XRoHJlhTkAXTrtvHfixeHq0Mk3f3jH/DSS3DRRQtp3Dh0NelLQR7Apq2K7t2t1SEiv/b113DFFdCpE5x44lehy0lrCvLAXnwRHnkkdBUi6eeyy2D1alvXPxt3/dkR+vYE1qGDrZK4bFnoSkTSx8SJNkLlpptg331DV5P+FOSBPfSQravco4e6WETANmPp0QMOOsj2wJXtU5AHts8+tvjP00/DY4+FrkYkvMsvh++/hzFjtM54RSnI00Dv3nDIIbY05zffhK5GJJxnnoFHH7U1VQ44IHQ18aEgTwM5Odb6KCmBiy9WF4tkp//8BwoLoW1bW2ROKk5BniZatbI1JJ5+Gv75z9DViFQu7+Gvf4UffrCx4+pS2TEK8jRyxRVw2GE27OrLL0NXI1J5Hn3URqoMGAD77x+6mvhRkKeRX7pYNmyACy6A0tLQFYmk3tKltjRthw7Qp0/oauJJQZ5mmje3fQhfftnWXhbJZKWl1mhZtw7GjrXGjOw4BXkauvhiOO446NsX5s8PXY1I6tx9t+2aNXy4bbwiiVGQpyHnbKJQ7dq2wNbataErEonevHlw9dVwwgnWeJHEKcjTVF6ehfl778ENN4SuRiRaa9daI2XnnWH0aC1PmywFeRrr0sVWRxw2DKZPD12NSHSuvRbef99u7jdsGLqa+FOQp7lhw6BlSzjnHJswIRJ3L74II0bYSJXjjgtdTWaIJMidc8c45z52zn3mnLsmimNmsnHjxv3v382aNfvV883VqgVFRbBiBZx/vmZ9SrwtWwZ//rNNvx82LHQ1mSPpIHfO5QD3AMcCrYCznHOtkj1upho3bhyFhYX/e/7FF19QWFi4zTBv08Z+6Z9/Hu66qzKqFIleaaltPr56tTVOatQIXVEG8d4n9QA6AFM2ed4P6Let1+Tn5/tEAXrooYcesX0UFxcnk3+zysvUqiRvd2DTCeVLgPabf5FzrhAoBMjLy2PGjBkRnFpEJF5KSkqiz79ttZwr8gBOBUZv8vxcYOS2XpNMizyZv2bpoGnTpuX+lW7atGmFXv/6697n5Hh/yinel5amttYdEfefyy8y5Tq8T59r+fpr7xs18r5FC+9XrkzsGOlyLVFIRYs8ipudS4EmmzxvXPYxKcegQYOoWbPmrz5Ws2ZNBg0aVKHXd+oEQ4bYAkN3352KCkWis2GDjbj67jvbuq1OndAVZaYogvwdoIVzbg/nXHXgTODZCI6bkbp168aoUaP+97xp06aMGjWKbt26VfgYV15ps+H69IG3305FlSLRGDwYXnrJGh1t2oSuJnMlHeTe+/XApcAUYAHwuPdeK4Rsw6ahvWjRoh0KcbBZcA8/DI0awamnany5pKcpU2zz5HPOgQsvDF1NZotkHLn3/gXv/d7e++be+4r1EUhS6tWDJ5+0reHOOgvWrw9dkchGixbB2WdD69Zw332agp9qmtkZYwcdBPfeC9Om2R6HIungp5+ga1frH3/ySZvUJqkVxfBDCej882HmTLjtNtvAuWvX0BVJNvMeevSAd9+F556DvfYKXVF2UIs8A/ztb9C+PfzpT/DBB6GrkWw2cqTdv7nhBlv0TSqHgjwD7LSTvYWtUwf++Ecb6iVS2aZNs31nTzwR+vcPXU12UZBniEaNLMyXLoXTT9fNT6lcCxfa713LlvDPf0IVJUul0rc7gxx6KIwaZWuX9+4duhrJFj/+aO8EvYdnn9WknxB0szPD/PnPtmD/iBGwzz5240kkVdavt+GvCxbA5Mm2ebhUPgV5Brr9dvjkE7j8cvsf6//+L3RFkqn69IEXXrBhsEcdFbqa7KWulQyUkwPjx8N++1m/5XzNs5UUuPdeWx+/Vy/bklDCUZBnqNq1YdIkqFkTjj8eli8PXZFkkhdegMsusyGG2uknPAV5BmvSxML8228tzFetCl2RZIJZs+C002wRrPHj7R2ghKUgz3D5+fD44/Dee9bN8vPPoSuSOFu40BoFDRva1oO1a4euSEBBnhWOO84WLpo8GS6+GG3gLAn59ls49lhrDLz4Iuy6a+iK5BcatZIlLrwQvvwSBgyw1tSQIaErkjhZtcoaBIsXw9SpNrRV0oeCPIv072/L3t52GzRoYEPHRLZn7Vo4+WSYM8dmDx92WOiKZHMK8izinC1qtGIFXHUV1K9vqyeKbM0vW7VNmwZjx9o6KpJ+FORZJifH1sL44QfrbqlVy26CimyutNR+R554AoYPt9U1JT3pZmcWql7d3iJ37Ajdutn6GCKb8h4uvdSWpO3fX2v3pDsFeZaqVcuGjx14oI0JnjIldEWSLry3+yf33gt9+8KNN4auSLZHQZ7F6tSxIYn77gsnnWS7nUt28x769bNF1y67zEY3ab/N9Kcgz3L16lmA77233ciaOjV0RRKK93D11TaqqXt3uPNOhXhcKMiFXXaxUQktW1qYq5sl+3hvI5mGDoVLLoG//12bQ8SJflQC2LjyadOsm+XEE+GZZ0JXJJWltNSWPB4+3LpT7r5bLfG4UZDL/9Svb2Heti2ccootiCSZbf16uOACm1/Qu7ctS6sQjx8FufxKvXrw8ss2e++cc+D++0NXJKmybh2ceaZN9Ln5ZluOViEeTwpy2ULt2rbe9LHH2k2vQYO00Fam+fFHW8Vw4kQboXLjjQrxOFOQS7l+8xt46imbMHT99TY5ZMOG0FVJFJYvh86dobgYxoyBK64IXZEkS1P0ZauqV4d//AMaNbLRDMuXwyOPWMhLPH36qe3h+vXX8Nxz9q5L4k8tctmmKlVsM+c777QWekGBhYDEzyuvwKGH2pK0xcUK8UyiIJcK6dnT+lPffx/at4d580JXJDti7Fj4wx8gLw/eegsOOSR0RRKlpILcOTfUOfeRc+5959xTzrm6EdUlaejkk+G112y0Q8eOtlaLpLcNG2zK/XnnweGHwxtvwJ57hq5KopZsi/wloLX3/gDgE6Bf8iVJOsvPtxbdXnvBCSfALbfYhBJJP6tWVaVLF1sv5eKLbXu2unVDVyWpkFSQe++neu/Xlz2dCTROviRJd02awOuvw9lnww032OqJa9ZoK/V0Mm8edO+ez7RpNhfgvvugWrXQVUmqOB/RAGHn3HPAY977R7by+UKgECAvLy+/qKgoofOUlJSQm5ubcJ3poqCgAIDi4uLAlSTOe5gwoTH339+c3XZbzc03L6B589Why0pKJvx+TZ6cx5137k3Nmj8zYMCHtG79Y+iSkpYJP5dfJHMtBQUFs7337bb4hPd+mw/gZWBeOY8/bvI11wFPUfaHYXuP/Px8n6ji4uKEX5tOAG/f/vibMcP7+vX/62vU8P7BB70vLQ1dUeLi/Pu1erX3F1zgPXjfubP3Eyf+K3RJkYnzz2VzyVwLMMuXk6nb7Vrx3h/lvW9dzuMZAOfceUAXoFvZiSTLHHEEjBo1i9//Hv7yF5va/8MPoavKLu+/byNRxoyx7q6XX4Z69daFLksqSbKjVo4B+gIneu/XRFOSxFG9ej8zeTIMHAiPPQZt2sCrr4auKvOVltoU+4MPhu++s41CBgywvVkleyQ7amUkUBt4yTk31zl3XwQ1SUzl5Nh0/n/9y26sde5sGxX897+hK8tMixbB0UfDlVfa5J4PPrDnkn2SHbWyl/e+ife+bdmje1SFSXy1bw9z51o3y+2327K4b74ZuqrMUVpqGz+0bm1DQUeNslm3DRqErkxC0cxOSYncXHjgAdtt6KefoFMnmx36Y/wHUAS1YAEceaTt4tOxow0zvOgirVyY7RTkklJHH21v+f/6V9t5Zt994fHHtSzujlqzBq691u49vPfexj+STZuGrkzSgYJcUq5OHbjnHpg5E3bdFc44Y2PAy7Z5b3/4WrWCW2+1SVgffwwXXqhWuGykIJdKc8gh8Pbb1jKfPdv6zi++WKspbs3bb9tOTWecATvvDDNmwMMPQ8OGoSuTdKMgl0qVk2ObVHz2mW30+9BD0Lw5XHcdfP996OrSw/z5tmdq+/b2fXrgAZgzx8bri5RHQS5B1Ktna5zPnw9dusDgwbYq38CB2RvoH34I554L++8PL70EN90En3xi3SgaFy7boiCXoPbeG4qKbLjiYYfZ3pG/+x307QtffRW6usrx9tu2RPB++9ma7336wL//Df372/0Fke1RkEtaaNMGnn3WAr1LFxg+HJo1sz1DZ87MvFEu69bB+PHQoYN1ocyYYVPrFy+2sff164euUOJEQS5ppU0bePRR61Lo0QMmTbKwa9fORr6sWBG6wuR8/DFcc40NG+zWDb791rqYFi+2qfWa1COJUJBLWmre3AJuyRIL8A0b7CbpbrvB6afDE0/Y2Oo4WLrURup07Aj77APDhtnaKC+8YMHesyfUrh26SomzqqELENmW2rWtZd6jh3W7jBljLfYJE6BmTVtj5PjjbWf4Ro1CV2u8t9UIJ0+27qI33rCPt24NQ4fa6pC77hq2RsksCnKJjbZt4a67rP/8tdesVf7003aDEOCAA2yI3u9/b4/KCvbSUps6/9prtnNScfHGG7Vt29pInFNPtda4SCooyCV2qlaFggJ7jBxpM0RffBGmToUHH7RuDLBW74EHWpi2bGn7jO61l02oSWRW5M8/w5df2tjuzz6z4YLvvmtT5lev3njOww+HY45Jr3cJktkU5BJrzllL/IADbMncn3+2cH3jDfvv3Lk2Jnv9+o2vqVYN8vIsdOvWtQW+cnPt48uWtWTsWFi7FkpKYNUqW+d7+XK7Mbnp6Jnate3m7AUXwEEH2buA5s01dV4qn4JcMkq1arYUwCGHbPzYunXwxRfw6afWkl62zIJ5+XJYudKWCFi1ym6orl37W3baCWrU2Bjwe+xhNyp33RUaN4YWLeyx224KbUkPCnLJeNWrbwzf7ZkxYyadO3dOeU0iUdLwQxGRmFOQi4jEnIJcRCTmFOQiIjGnIBcRiTkFuYhIzCnIRURiTkEuIhJzCnIRkZhTkIuIxJyCXEQk5hTkIiIxpyAXEYk5BbmISMwpyEVEYi6SIHfOXemc8865BlEcT0REKi7pIHfONQGOBhYnX46IiOyoKFrkdwB9Ab+9LxQRkeg57xPPX+fcH4Ejvfc9nXOLgHbe+2+38rWFQCFAXl5eflFRUULnLCkpITc3N8GK04uuJf1kynWAriVdJXMtBQUFs7337bb4hPd+mw/gZWBeOY8/Am8BO5d93SKgwfaO570nPz/fJ6q4uDjh16YbXUv6yZTr8F7Xkq6SuRZgli8nU7e7+bL3/qjyPu6c2x/YA3jP2VbijYE5zrlDvPfLd/APjYiIJGi7Qb413vsPgIa/PN9e14qIiKSGxpGLiMRcwi3yzXnvm0V1LBERqTi1yEVEYk5BLiIScwpyEZGYU5CLiMRcUjM7Ez6pc/8Bvkjw5Q2ATBniqGtJP5lyHaBrSVfJXEtT7/0um38wSJAnwzk3y5c3RTWGdC3pJ1OuA3Qt6SoV16KuFRGRmFOQi4jEXByDfFToAiKka0k/mXIdoGtJV5FfS+z6yEVE5Nfi2CIXEZFNKMhFRGIulkHunDvNOTffOVfqnIvdkCTn3DHOuY+dc585564JXU8ynHMPOee+cc7NC11LMpxzTZxzxc65D8t+t3qGrilRzrkazrm3nXPvlV3LzaFrSoZzLsc5965zblLoWpLhnFvknPvAOTfXOTcrymPHMsixHYq6Aq+GLmRHOedygHuAY4FWwFnOuVZhq0rKw8AxoYuIwHrgSu99K+BQ4JIY/1zWYlswtgHaAsc45w4NW1JSegILQhcRkQLvfVuNIwe89wu89x+HriNBhwCfee8Xeu/XAUXYtnmx5L1/FVgRuo5kee+Xee/nlP17FRYcu4etKjFlu4KVlD2tVvaI5agG51xj4HhgdOha0lksgzzmdge+3OT5EmIaGJnKOdcMOBDbkzaWyroj5gLfAC957+N6LXcCfYHSwHVEwQNTnXOzyzajj0xkG0tEzTn3MrBrOZ+6znv/TGXXI9nBOZcLTAR6ee9/DF1Porz3G4C2zrm6wFPOudbe+1jdx3DOdQG+8d7Pds51DlxOFH7vvV/qnGsIvOSc+6jsHW3S0jbIt7bpcwZYCjTZ5Hnjso9JYM65aliIj/PePxm6nih4739wzhVj9zFiFeRAJ+BE59xxQA2gjnPuEe/9OYHrSoj3fmnZf79xzj2FdbNGEuTqWql87wAtnHN7OOeqA2cCzwauKes55xzwILDAez8idD3JcM7tUtYSxzn3G+APwEdBi0qA976f975x2TaSZwLT4xrizrlazrnav/wbOJoI/7DGMsidcyc755YAHYDnnXNTQtdUUd779cClwBTshtrj3vv5YatKnHPuUeBNoKVzbolz7i+ha0pQJ+Bc4Miy4WFzy1qCcbQbUOycex9rOLzkvY/10L0MkAe87px7D3gbeN57Pzmqg2uKvohIzMWyRS4iIhspyEVEYk5BLiIScwpyEZGYU5CLiMScglxEJOYU5CIiMff/AchuZFY26exuAAAAAElFTkSuQmCC", "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": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[2 - √3, √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": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.26794919]\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": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAla0lEQVR4nO3dd3hUZf7+8fcnjSodQyc0RTokSg9EEAEVUHEFUWyIugqKupbVrz9119W1AApYsIEVkbUAgkhJiDQlICC9hqZYEJDQy/P7I+NuxFAmM8mZzNyv65qLOTNn5rkZ9M7kOc2cc4iISPiL8jqAiIgUDBW+iEiEUOGLiEQIFb6ISIRQ4YuIRIgYrwOcSoUKFVxCQkKeXrtv3z5KlCgR3EBBoFz+US7/KJd/wjHXokWLfnHOVcz1SedcyN4SExNdXqWmpub5tflJufyjXP5RLv+EYy4gw52kUzWlIyISIVT4IiIRQoUvIhIhVPgiIhFChS8iEiGCUvhm9qaZ/WRmy0/yvJnZi2a23syWmVmLYIybq/feg4QEOlx4ISQkZC+LiEjQvuGPAbqe4vluQD3fbSDwcpDG/aP33oOBA2HzZsw52Lw5e1mlLyISnMJ3zqUDv55ilZ7A277dRBcAZcyscjDG/oOHH4b9+//42P792Y+LiEQ4c0E6H76ZJQCTnXONcnluMvC0c26Ob3km8IBzLiOXdQeS/VsA8fHxiePGjTvjDB0uvBBzDvMt//43c2bMnjXLn79OvsnKyqJkyZJex/gT5fKPcvlHufwTSK6UlJRFzrmkXJ882RFZ/t6ABGD5SZ6bDLTLsTwTSDrde/p9pG3Nms5l97wDnPPdDler7t/75KNwPLIvPymXf5TLP+GYixA40nY7UD3HcjXfY8H15JNQvPgfHjoQW4QHk/rw9vxMjh/X1b1EJHIVVOFPBPr79tZpBexxzv0Q9FH69YPRo/+3XLMmh196hV8u682jn63g+re+4Yc9B4I+rIhIYRCs3TI/AOYD55rZNjO72cxuM7PbfKtMATYC64HXgL8GY9xc9ev3v/uZmZQecANjbjyff/ZqREbmLroMS+fTb7f/PrUkIhIxgnJ6ZOdc39M874A7gjFWXpgZ17aqSbu6Fbj3o6Xc/eESpq3YwZOXN6ZciTivYomIFKiIOtI2oUIJxt/amge61mfGqh/pMiydGSt/9DqWiEiBiKjCB4iOMm7vWIeJd7aj4llFGPB2Bvd9tJTfDh7xOpqISL6KuML/3XmVS/HZHW25M6UuHy/eRrfhXzF3/S9exxIRyTcRW/gAcTFR3Hfxufzn9jYUiY2i3+tf8/8+W87+w0e9jiYiEnQRXfi/a16jLJ8Pas9NbWsxdv5mur/wFYs2n+pMESIihY8K36dYXDSPXtaAD25pxdHjjt6vzOepKas4eOSY19FERIJChX+C1nXK88XdyfQ5vwavpm/k0hFzWLp1t9exREQCpsLPRckiMTx1RWPG3nQBWQePcsXL83hu2hoOHz3udTQRkTxT4Z9Ch3MqMm1IMpc3r8rI1PX0GDmH5dv3eB1LRCRPVPinUbpYLM9d1ZTX+yexc99heo2ay/AZazlyTN/2RaRwUeGfoc4N4pk+JJlLm1Rm+Ix19Bo1l1U//OZ1LBGRM6bC90OZ4nEM79OcV65N5MffDtJj5BxGzFynb/siUiio8POga6NKfDmkA10bVeb56Wu54qV5rNmx1+tYIiKnpMLPo3Il4hjRtzkv92vB97sPcOmIrxg5S9/2RSR0qfAD1K1xZabf04GLG1biuS/XcvlLc1m9Q3P7IhJ6VPhBUK5EHCOvacHL/VqwY89BLhuhuX0RCT0q/CDq1rgyXw7pQDff3H6vUXNZ+b2+7YtIaFDhB1m5EnG82Pf3PXkO0WPkHIZNX6ujdEXEcyr8fNK1USWmD0nmsqZVeGHmOnqMnMN323SUroh4R4Wfj8qWiGPY1c14vX8Sv+47TK+X5jJh7WGdgVNEPKHCLwCdG8Qz/Z4OXNG8KpM3HuHSEXNYvGWX17FEJMKo8AtI6WKxPHtVU+5NLML+Q0e58uV5/GPySg4c1rd9ESkYKvwC1rhiDNOGJNOvZQ3emLOJi4enM3/DTq9jiUgEUOF74KyisfyzV2PGDWxFlEHf1xbw90++Y+/BI15HE5EwFpTCN7OuZrbGzNab2YO5PH+Dmf1sZkt8twHBGLewa1W7PFPvSuaW9rUY980WugxLJ3X1T17HEpEwFXDhm1k0MAroBjQA+ppZg1xW/dA518x3ez3QccNFsbhoHr6kAR//tS1nFY3hxjELuXvct/y677DX0UQkzATjG/4FwHrn3Ebn3GFgHNAzCO8bUZpVL8OkQe24q1M9Ji/7gYuGzmbS0u9xznkdTUTChAVaKGbWG+jqnBvgW74OaOmcuzPHOjcATwE/A2uBIc65rSd5v4HAQID4+PjEcePG+Z0pJSUFgNTUVL9fm9+ysrIoWbLkKdfZuvc4by4/xKY9x2lWMZr+DeMoVzR/N7ecSS4vKJd/lMs/4ZgrJSVlkXMuKdcnnXMB3YDewOs5lq8DRp6wTnmgiO/+rcCsM3nvxMRElxeAy/6rhZ7U1NQzWu/osePutfQN7txHprhGj37h3l2Q6Y4dO+55roKmXP5RLv+EYy4gw52kU4PxtXE7UD3HcjXfYzl/qOx0zh3yLb4OJAZh3LAWHWUMaF+baXcn07haaR7+ZDl9X1vAxp+zvI4mIoVUMAp/IVDPzGqZWRzQB5iYcwUzq5xjsQewKgjjRoSa5Uvw3oCWPHNlE1b98BtdX/iKUanrdeplEfFbwIXvnDsK3AlMI7vIxzvnVpjZE2bWw7faYDNbYWZLgcHADYGOG0nMjL+cX50Z93SgU/2zeXbaGnqMnMuybbu9jiYihUhMMN7EOTcFmHLCY4/muP8Q8FAwxopkZ5cqysvXJvLF8h08+tlyeo2ay01ta3FPl3MoHheUf0oRCWM60rYQ6tqoEtPv6cDV59fg9Tmb6DIsnfS1P3sdS0RCnAq/kCpdLJanrmjMhwNbERcdRf83v2HIh0t0wJaInJQKv5BrWbs8U+5qz+AL6zJ52fd0ej6Njxdv0wFbIvInKvwwUDQ2mnu6nMvkQe1JqFCCe8Yvpf+b37Bl536vo4lICFHhh5FzK53FhNva8HiPhny7ZTddhs/mldkbtAuniAAq/LATHWVc3yaB6fck075eRZ6eupoeI+eydOtur6OJiMdU+GGqculivNY/iVeuTeTXfYfo9dJcHpu4gqxDR72OJiIeUeGHud934byuVU3Gzs/koqGzmbZih9exRMQDKvwIUKpoLE/0bMR/bm9D6WKx3PrOIga+ncEPew54HU1ECpAKP4K0qFGWSYPa8UDX+qSv+5nOz8/mrbmbOHZcu3CKRAIVfoSJjY7i9o51+PLuDiQmlOPxSSu5/KW5ZO455nU0EclnKvwIVaN8ccbeeD4j+jbn+90HeXz+QZ6YtFIbdUXCmAo/gpkZlzWtwsx7O9CxegxvzdukjboiYUyFL5QuFsv1DYv8YaPugLEL2bZLR+qKhBMVvvzX7xt1/969PnPX7+Sioem8qiN1RcKGCl/+IDY6ioHJdZhxbwfa1q3AU1NXc+mLc8jI/NXraCISIBW+5KpqmWK8fn0So69LZO/BI/R+ZT4PTFjGLp1+WaTQUuHLKXVpmH2k7q3JtZmweBsXPp/G+IytHNe++yKFjgpfTqtEkRge6n4enw9uR52KJbl/wjKuHj2fNTv2eh1NRPygwpczVr9SKcbf2ppnrmzC+p+y6P7iV/xryir2ad99kUJBhS9+iYoy/nJ+dWbd25GrEqsxOn0jnYfOZup3P+gqWyIhToUveVK2RBxPX9mE/9zehjLF47j9vcVc/9ZCMn/Z53U0ETkJFb4EJLFmWSbd2ZZHL23A4s276DI8naHT13LwiM7NIxJqVPgSsJjoKG5qV4tZ93aga8NKvDhzHRcNm83MVT96HU1EcghK4ZtZVzNbY2brzezBXJ4vYmYf+p7/2swSgjGuhJazSxXlxb7Nef+WlhSJiebmsRkMGJvB1l91igaRUBBw4ZtZNDAK6AY0APqaWYMTVrsZ2OWcqwsMA/4d6LgSutrUqcCUwe15sFt95m34hc5DZzNi5jpN84h4LCYI73EBsN45txHAzMYBPYGVOdbpCTzmuz8BGGlm5vJ5tw4zy8+3Fz8MfhIGex1CpJBITU3Nl/cNRuFXBbbmWN4GtDzZOs65o2a2BygP/HLim5nZQGAgQHx8PGlpaUGIKCJSeGRlZeVL9wWj8IPKOTcaGA2QlJTkOnbsGMh7BSlV8KSlpRHI3ym/FESuw0eP88acTbw4cx3HneOvHetya4faFI2N9jRXXiiXf5TLP/mVKxgbbbcD1XMsV/M9lus6ZhYDlAZ2BmFsKUTiYrIvrzjz3g50Pi+eYTPW0mVYuvbmESkgwSj8hUA9M6tlZnFAH2DiCetMBK733e8NzMrv+XsJXVXKFGNUvxa8N6AlcTFR3Dw2g5vG6KAtkfwWcOE7544CdwLTgFXAeOfcCjN7wsx6+FZ7AyhvZuuBe4A/7bopkadt3ey9ef7evT5fb9xJl2HpPDdtDfsP69w8IvkhKHP4zrkpwJQTHns0x/2DwFXBGEvCS1xM9gVXejWrytNTVzMydT0fL97Gw5c0oHvjSl7HEwkrOtJWQsLZpYoy9OpmTLitNWWKx3HH+4u55rWv2b5Xl1cUCRYVvoSUpIRyTBrUjn/2asTKH37j/+Yd4LGJK9iz/4jX0UQKPRW+hJzoKOPaVjVJu68jHavF8Pb8TFKeT+ODb7ZwTFfaEskzFb6ErLIl4ujfsAiTBrWjTsUSPPTxd/QcNYdFm3VBdZG8UOFLyGtYpTTjb23NC32a8cvew1z58nzuHvctO/Yc9DqaSKGiwpdCwczo2awqM+/twB0pdZiyfAcXPp/GqNT1OimbyBlS4UuhUqJIDH+7uD4zhnSgXd0KPDttDV2GpTNtxY6QPJWGSChR4UuhVKN8cUb3T+Ldm1tSNDaKW99ZxLVvfM2aHXu9jiYSslT4Uqi1q5d9tO4TPRuyfPtvdHshnUc/W86ufYe9jiYSclT4UujFREfRv3UCafd15LpWNXnv6y10fC6Nt+Zu4sgxHbgl8jsVvoSNsiXieLxnI6be1Z7GVUvz+KSVdHvhK9LW/OR1NJGQoMKXsHNO/Fm8c/MFvNY/iaPHjnPDWwu58a1vWP9TltfRRDylwpewZGZc1CCeaUOS+Xv3+mRk7qLr8HQen7SC3fs1vy+RSYUvYa1ITDQDk+uQ+reO/OX86oydl0nH59IYo/l9iUAqfIkIFUoW4V+XN+bzwe1pWKUUj01aSdfh6cxa/aP235eIocKXiHJe5VK8e3NLXuufxHEHN43JoP+b32j/fYkIKnyJOP+d3787mf+7tAFLt+6m2wvp/P2T7/gl65DX8UTyjQpfIlZcTBQ3t6vF7L+l0L91AuMXbqXjs2m8nLZB5+eRsKTCl4hXtkQcj/VoyLQhybSqXZ5/f7GaTs/PZuLS7zW/L2FFhS/iU6diSV6/Pon3B7SkdLFYBn/wLZe/NE/n35ewocIXOUGbuhWYNKgdz/Ruwve7D3Dly/P563uL2Lxzn9fRRAIS43UAkVAUHWX8Jak6lzapzOj0jbw6eyPTV/5I/9YJDLqwrtfxRPJEhS9yCsXjYri78zn0vaAGQ79cy5tzNzFh0Ta61zRaHz1GkZhoryOKnDFN6YicgfhSRfl37yZMGdyeptXL8MHqw1w0NJ3Jy7RhVwqPgArfzMqZ2XQzW+f7s+xJ1jtmZkt8t4mBjCnipfMql+Ltmy7gvqQiFI+L5s73szfsLszUhl0JfYF+w38QmOmcqwfM9C3n5oBzrpnv1iPAMUU816hCDJ8Pbs8zvZvww54DXPXKfAa+ncGGn3VGTgldgRZ+T2Cs7/5YoFeA7ydSaPy+YTftvhTu63IO8zbspMuwdB759Dt+3qsjdiX0WCDzj2a22zlXxnffgF2/L5+w3lFgCXAUeNo59+kp3nMgMBAgPj4+cdy4cX7nSklJASA1NdXv1+a3rKwsSpYs6XWMP1Eu/+SW67dDjs82HCZt61Fio6BbrVguToilaIx5misUKJd/AsmVkpKyyDmXlOuTzrlT3oAZwPJcbj2B3Sesu+sk71HV92dtIBOoc7pxnXMkJia6vABc9l8t9KSmpnodIVfK5Z9T5dr4c5a7/d0MV/OByS7xH9PdO/Mz3eGjxzzP5SXl8k8guYAMd5JOPe1umc65zid7zsx+NLPKzrkfzKwykOu15Jxz231/bjSzNKA5sOF0Y4sURrUqlOClfoks2ryLp6as4pFPl/PmnE3c37U+FzeMJ/uXYZGCF+gc/kTget/964HPTlzBzMqaWRHf/QpAW2BlgOOKhLzEmmX56LbWjL4uETO47d1FXPmy9ugR7wRa+E8DF5nZOqCzbxkzSzKz133rnAdkmNlSIJXsOXwVvkQEM6NLw0pMuzuZp65ozPbd2Xv0DBi7kLU/6hz8UrACOtLWObcT6JTL4xnAAN/9eUDjQMYRKexioqPoe0ENejWryptzN/FK2ga6Dk/nyhbVGHLROVQpU8zriBIBdKStSAEqFhfNHSl1Sb8/hZva1uKzJd/T8bk0/jVllS6uLvlOhS/igbIl4njk0gbMuq8DlzWpwmtfbaT9M6mMSl3P/sNHvY4nYUqFL+KhamWL8/xfmvLFXcm0rFWOZ6etocOzabyzYDNHjh33Op6EGRW+SAg4t9JZvH79+Uy4rTUJ5Yvzf58up/PQ7KtuHT+uk7NJcKjwRUJIUkI5xt/amjdvSKJYbDSDP/iWS0fMIXXNTzorpwRMhS8SYsyMC+vH8/ng9gy/uhl7Dx3hxrcWcvWrC8jQPvwSABW+SIiKjjJ6Na/KzHs68o+eDdm0cx+9X5nPjW99w4rv93gdTwohFb5IiIuLieK61gnM/ltH7u96Los27+KSF+cw6INv2ajTMYsfVPgihUTxuBj+2rEuXz1wIXem1GXmqh+5aFg6D/5nGdt3H/A6nhQCuqatSCFTulgs9118Lte3SWBU6nre/3oLHy/ezjUta9C8iDbsysmp8EUKqYpnFeGxHg25Jbk2I2au450Fm3nfHKvdam5Nrk2Z4nFeR5QQoykdkUKuapliPH1lE2bc04EWZ0fzyuwNtP93Ki/MWMfeg0e8jichRIUvEiZqVSjBbU2LMvWu9rSuU55hM9aS/Ewqr87ewIHDx7yOJyFAhS8SZupXKsXo/klMvLMtTauX4ampq2n/TCpvztnEwSMq/kimwhcJU02qlWHMjRcw4bbW1D27BE9MXknHZ9N4d8FmDh/VeXoikQpfJMwlJZRj3MDWvD+gJVXLFuORT5eT8lwaHy7cohO0RRgVvkiEaFO3AhNua82YG8+nQsk4HvjPd3R6fjYTFm3jqIo/IqjwRSKImdHx3LP59I62vHF9EmcVjeG+j5Zy0bB0Pv12O8d0Zs6wpsIXiUBmRqfz4pk8qB2vXpdIkZgo7v5wCV2GzeazJSr+cKXCF4lgZsbFDSsxZXB7XurXgpioKO4at4SLh6czcen3Kv4wo8IXEaKijO6NKzP1rvaMuqYFUQaDP/iWrsPTmaTiDxsqfBH5r6go45ImlfnirmRG9G2OAwap+MOGCl9E/iQqyrisaRWm3Z3Mi32bA9nFr6mewk2FLyInFR1l9PAV/8hrmv93qqfLsNnaq6cQCqjwzewqM1thZsfNLOkU63U1szVmtt7MHgxkTBEpeFFRxqVNqvDFXcmMuiZ74+7dHy7hoqGz+Xix9uMvLAL9hr8cuAJIP9kKZhYNjAK6AQ2AvmbWIMBxRcQDv8/xT72rPS/3a0GR2GjuGb+UTkNnMz5jq47cDXEBFb5zbpVzbs1pVrsAWO+c2+icOwyMA3oGMq6IeCsqyujWuDJTBrfjtf7ZB3DdP2EZKc+l8f7XWzh0VCdpC0XmXOBzcGaWBtznnMvI5bneQFfn3ADf8nVAS+fcnSd5r4HAQID4+PjEcePG+Z0nJSUFgNTUVL9fm9+ysrIoWbKk1zH+RLn8o1x/5Jxj6c/HmLjhCBv3HKdcUaN7rViSq8UQF236vPwUSK6UlJRFzrlcp9hPe8UrM5sBVMrlqYedc5/lKdEpOOdGA6MBkpKSXMeOHfP8XoG8Nr+kpaUplx+Uyz9e5koB7nKOr9b9wohZ63h31S6mbTMGtq9N9aKb9Xn5Ib9ynbbwnXOdAxxjO1A9x3I132MiEmbMjORzKpJ8TkUWbNzJiFnreHLKKkrGwvqodfRvk0CporFex4xYBbFb5kKgnpnVMrM4oA8wsQDGFREPtapdnvcGtOI/t7ehTplonvtyLW2fnsVz09bw677DXseLSIHulnm5mW0DWgOfm9k03+NVzGwKgHPuKHAnMA1YBYx3zq0ILLaIFBaJNcsyJLEokwe1o13dCoxKW0/bp2fxz8kr+fG3g17HiyinndI5FefcJ8AnuTz+PdA9x/IUYEogY4lI4daoamlevjaRdT/u5eW0Dbw1L5O352+md1I1bkuuQ43yxb2OGPZ0pK2IFKh68Wcx9OpmpN7bkd5J1ZiQsY2U59MY8uES1v641+t4YU2FLyKeqFG+OP+6vDFfPZDCjW0SmLZiB12GpXPL2xl8u2WX1/HCUkBTOiIigYovVZRHLm3AHSl1GTMvkzHzMpm+8kfa1CnPXzvWpW3d8piZ1zHDgr7hi0hIKFsijiEXncPcBy/k4e7nsf6nLK5942t6jprL1O9+4LhO1BYwFb6IhJSSRWK4Jbk2Xz2Qwr8ub8xvB45w+3uL6TxsNuMXbtVpGwKgwheRkFQkJpprWtZg5r0dGXVNC4rFRnP/f5aR/Ewqo9M3kHXoqNcRCx0VvoiEtGjfGTonD2rHOzdfQJ2KJfnXlNW0eWomz05bzc97D3kdsdDQRlsRKRTMjPb1KtK+XkWWbN3Nq7M38FLaBl77ahNXJVbjlva1SahQwuuYIU2FLyKFTrPqZXj52kQ2/pzFa19t5KOMbbz/zRa6NarErcl1aFq9jNcRQ5IKX0QKrdoVS/LUFU0YctE5jJmbyTsLNjPlux20ql2OW5Pr0PHcitqlMwfN4YtIoXf2WUW5v2t95j14IY9cch6bd+7nxjELuXh4Oh9lbOXwUV2JC1T4IhJGzioay4D2tUm/P4Whf2lKlBl/m7CM9s/M4pXZG9hz4IjXET2lKR0RCTux0VFc0aIalzevSvq6XxidvoGnp65m5Kz19Dm/Oje2q0XVMsW8jlngVPgiErbMjA7nVKTDORVZvn0Pr3+1kbfmZfLWvEwuaVyZFsUj6yAuTemISERoVLU0w/s0J/3+FG5qm8Cs1T/x2PyDXP3qfGas/DEiTt2gwheRiFK1TDEevqQB8x66kKvPjWPrr/sZ8HYGnYfN5r2vN3PgcPh+61fhi0hEKlU0lm61Ypl9fwov9GlGibgYHv5kOW2ensnzX67hp73hdzUuzeGLSESLjY6iZ7Oq9GhahW82/cobczYxMnU9r8zewGVNq3Bzu1o0rFLa65hBocIXESF7A2/L2uVpWbs8mb/s4625m/ho0TY+Xryd1rXLc1O7WnSqfzZRUYX3QC5N6YiInCChQgke79mI+Q924sFu9cncuY9b3s7gwufTGDN3E/sK6Zk6VfgiIidRungst3WoQ/r9KYzo25wyxeN4bNJKWj01kyc/X8nWX/d7HdEvmtIRETmN2OgoLmtahcuaVmHxll28NTeTN+dm8sacTXRpUIkb2yZwQa1yIX/eHhW+iIgfWtQoS4saZXmoW33eWbCZD77ZwhcrdtCgcilubJvAZU2rUDQ22uuYudKUjohIHlQpU4wHutZn/oOdeOqKxhw9fpy/TVhGm6dn8dy0NezYE3q7dQZU+GZ2lZmtMLPjZpZ0ivUyzew7M1tiZhmBjCkiEkqKxUXT94IaTLs7mfcHtKRFjbKMSltPu3/P4s73F7No8684FxpH8QY6pbMcuAJ49QzWTXHO/RLgeCIiIcnMaFO3Am3qVmDLzv28PT+TDzO2MnnZDzSqWorrW3s/3RPQN3zn3Crn3JpghRERCQc1yhfnkUsbsOChTvyzVyMOHfnfdM8zX6xm++4DnuSyYPyqYWZpwH3OuVyna8xsE7ALcMCrzrnRp3ivgcBAgPj4+MRx48b5nSclJQWA1NRUv1+b37KysihZsqTXMf5EufyjXP6J9FzOOVb9epwZm4/w7U/Z5+ppER9NpxqxnFcu6k979wSSKyUlZZFzLvcpdufcKW/ADLKnbk689cyxThqQdIr3qOr782xgKZB8unGdcyQmJrq8IPsHS55em99SU1O9jpAr5fKPcvlHuf5n66/73NNTV7lmj09zNR+Y7Do9n+bGztvkfjtw2Ll333WuZk133My5mjWzl/0EZLiTdOpp5/Cdc53z9GPmj++x3ffnT2b2CXABkB7o+4qIFDbVyhbnga71uatTPSYv+4G352fy6Gcr+O6Zl3hy6gjiDh3EADZvhoEDs1/Ur19Qxs73/fDNrAQQ5Zzb67vfBXgiv8cVEQllRWOj6Z1Yjd6J1fh2yy5qJt5C3KETduXcvx8efjhohR/obpmXm9k2oDXwuZlN8z1excym+FaLB+aY2VLgG+Bz59wXgYwrIhJOmtcoS7mdO3J/csuWoI0T0Dd859wnwCe5PP490N13fyPQNJBxRETCXo0a2dM4uT0eJDrSVkQkFDz5JBQv/sfHihfPfjxIVPgiIqGgXz8YPRpq1sSZQc2a2ctBmr8HFb6ISOjo1w8yM5k9axZkZga17EGFLyISMVT4IiIRQoUvIhIhVPgiIhFChS8iEiFU+CIiEUKFLyISIVT4IiIRQoUvIhIhVPgiIhFChS8iEiFU+CIiEUKFLyISIVT4IiIRQoUvIhIhVPgiIhFChS8iEiFU+CIiEUKFLyISIVT4IiIRQoUvIhIhAip8M3vWzFab2TIz+8TMypxkva5mtsbM1pvZg4GMKSIieRPoN/zpQCPnXBNgLfDQiSuYWTQwCugGNAD6mlmDAMcVERE/BVT4zrkvnXNHfYsLgGq5rHYBsN45t9E5dxgYB/QMZFwREfGfOeeC80Zmk4APnXPvnvB4b6Crc26Ab/k6oKVz7s6TvM9AYCBAfHx84rhx4/KUJysri5IlS+bptflJufyjXP5RLv+EY66UlJRFzrmkXJ90zp3yBswAludy65ljnYeBT/D9ADnh9b2B13MsXweMPN24zjkSExNdXqWmpub5tflJufyjXP5RLv+EYy4gw52kU2NO99PCOdf5VM+b2Q3ApUAn32An2g5Uz7FczfeYiIgUoED30ukK3A/0cM7tP8lqC4F6ZlbLzOKAPsDEQMYVERH/BbqXzkjgLGC6mS0xs1cAzKyKmU0BcNkbde8EpgGrgPHOuRUBjisiIn467ZTOqTjn6p7k8e+B7jmWpwBTAhlLREQCoyNtRUQihApfRCRCqPBFRCKECl9EJEIE7Ujb/GBmPwOb8/jyCsAvQYwTLMrlH+Xyj3L5Jxxz1XTOVcztiZAu/ECYWYY72eHFHlIu/yiXf5TLP5GWS1M6IiIRQoUvIhIhwrnwR3sd4CSUyz/K5R/l8k9E5QrbOXwREfmjcP6GLyIiOajwRUQiRKEv/NNdIN3MipjZh77nvzazhBDJlWxmi83sqO+qYAXiDHLdY2YrfRemn2lmNUMk121m9p3vrKxzCuq6yKfLlWO9K83MmVmB7OJ3Bp/XDWb2s+/zWmJmA0Ihl2+dv/j+G1thZu+HQi4zG5bjs1prZrtDJFcNM0s1s299/092z+19ztjJroxSGG5ANLABqA3EAUuBBies81fgFd/9PmRfhjEUciUATYC3gd4h9HmlAMV9928Poc+rVI77PYAvQiGXb72zgHSyr+ucFAq5gBs4wyvLFXCuesC3QFnf8tmhkOuE9QcBb4ZCLrI33t7uu98AyAxkzML+Df9MLpDeExjruz8B6GRm5nUu51ymc24ZcDyfs/ibK9X972I2J7swvRe5fsuxWAIoiL0NzuS/L4B/AP8GDhZAJn9yFbQzyXULMMo5twvAOfdTiOTKqS/wQYjkckAp3/3SwPeBDFjYC78qsDXH8jbfY7mu47IvxrIHKB8Cubzgb66bgan5mijbGeUyszvMbAPwDDA4FHKZWQugunPu8wLIc8a5fK70TQNMMLPquTzvRa5zgHPMbK6ZLfBdNS8UcgHgm8KsBcwKkVyPAdea2TayrykyKJABC3vhSz4xs2uBJOBZr7P8zjk3yjlXB3gAeMTrPGYWBQwF7vU6Sy4mAQnOuSbAdP73W67XYsie1ulI9jfp18ysjJeBTtAHmOCcO+Z1EJ++wBjnXDWyLyr1ju+/uzwp7IV/JhdI/+86ZhZD9q9FO0MglxfOKJeZdQYeJvtaxYdCJVcO44Be+RnI53S5zgIaAWlmlgm0AiYWwIbb035ezrmdOf7tXgcS8znTGeUi+1vsROfcEefcJmAt2T8AvM71uz4UzHQOnFmum4HxAM65+UBRsk+sljf5vWEinzd6xAAbyf4V7PeNHg1PWOcO/rjRdnwo5Mqx7hgKbqPtmXxezcnekFQvxP4d6+W4fxmQEQq5Tlg/jYLZaHsmn1flHPcvBxaESK6uwFjf/QpkT2mU9zqXb736QCa+A1JD5POaCtzgu38e2XP4ec6X73+pAvjQupP9LWED8LDvsSfI/nYK2T8RPwLWA98AtUMk1/lkf9vZR/ZvHCtCJNcM4Edgie82MURyvQCs8GVKPVXxFmSuE9YtkMI/w8/rKd/ntdT3edUPkVxG9jTYSuA7oE8o5PItPwY8XRB5/Pi8GgBzff+OS4AugYynUyuIiESIwj6HLyIiZ0iFLyISIVT4IiIRQoUvIhIhVPgiIhFChS8iEiFU+CIiEeL/A6qCQqKZTACKAAAAAElFTkSuQmCC", "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": 5, "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": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuL0lEQVR4nO3dd3hUxdfA8e+kQwKCpIEgoYTeAgHpJIBIk14NdkBF9CcISFEEFQtNRFCKgr6KIkpvokBCCUUCIr2EjnSkhQ6Z949JDCqakGz25m7O53n2SfbuZvfMk+Rkcu6ZuUprjRBCCPtyszoAIYQQGSOJXAghbE4SuRBC2JwkciGEsDlJ5EIIYXMeVrypv7+/DgkJcep7Xr58GV9fX6e+p7O48tjAtccnY7MvK8a3cePGM1rrgL8ftySRh4SEEBcX59T3jImJISIiwqnv6SyuPDZw7fHJ2OzLivEppQ7d7biUVoQQwuYkkQshhM1JIhdCCJuTRC6EEDYniVwIIWxOErkQQticJHIhhLA5SeRCCGFzliwIcjalFACy97oQwhXJjFwIIWxOErkQQticJHIhhLA5SeRCCGFzksiFEMLmJJELIYTNSSIXQgibk0QuhBA2J4lcCCFszlaJXFZmCiHEP9kqkU9bf5gXvt7IkT+uWB2KEEJkGbZK5DdvJxKz+zQNR69g9M97uHrjttUhCSGE5WyVyJ+uVYRlr9ajUdlgxi7bS4NRMcz/7ZiUXIQQ2ZqtEjlAgTw5+LhzGN91r859Ob146dtf6ThpHduPXbA6NCGEsITtEnmyh4rmY8FLtRnWuhx7T17i0Y9XM2DWVs4mXLc6NCGEcCrbJnIAdzdF1EOFiekTyVM1i/B93BEiRsbw2ar93LiVaHV4QgjhFA5J5EqpKUqpU0qpbY54vXt1X05PBj9ahh9fqUPYg3l5Z+FOGo9ZSfSuUzBtWsoTQ0L+el8IIVyAo2bkXwCNHfRa6VY8MBdfPl2VKU+FAzC773CuP9M15QmHDkH37pLMhRAuxSGJXGu9EvjDEa+VUUop6pcK4sdX6jJsw7d437gGQHjyE65cgUGDLItPCCEcTTmqdU8pFQIs0FqX+5fHuwPdAYKCgqpMnz7dIe/7X+rVr4/SmnzAWqBE0nGtFCuWL8/093eWhIQE/Pz8rA4j07jy+GRs9mXF+CIjIzdqrcP/8YDW2iE3IATYlpbnVqlSRTtF4cJagwa0J2iddDvjl1cv23lCJyYmOieOTBYdHW11CJnKlccnY7MvK8YHxOm75FRbd62katgwyJkTgJtJh257+zCh2fN89canfBXVh11HzloXnxBCOIBrJ/KoKJg0KeV+4cK4f/4Zfb9+h0db1CB0fTTuYWFMfeNTTl+S/nMhhD05qv3wW0wZuqRS6qhS6llHvK5DREWlfH7wIERF4eXhRpvHH6H0lrVseOE1AmZOJ3L4Mib8uI1rN2X/FiGEvTiqa6Wz1jq/1tpTa11Qa/25I143s+Xx9eaxt1+kzIqF1CgeQO2nWjEvoj2LVmyT/VuEELbh2qWVNCoa4MfkJ6tyZf5CcrppqjWtzfvPvUfcwSzRUSmEEP9JEvkdqlUtSdMVs9g05Qe25Aji6Y+WMrbfOA6euWx1aEII8a8kkf+Nm5uiUceGTBnxFP3L5KTF5+9zsHoE4z+Zx7nLN6wOTwgh/kES+b/I4eVOVLfm5Nyzg8u169Gh35O0GjKbSSv3yQlRIUSWIok8FYH5ctPsi5Gc27qbImWKcK3/ID5p24t5Gw6SmCgnRIUQ1pNEnkYligTyxdPVqDOoB5G711K6SV2GvjSatftkQZEQwlqSyO9RWLO6VNyxntODhpLvYDydJ6/jlY9+ZM/JS1aHJoTIpiSRp4Obuxs1ez1N9znjeKNWft4Y1JnYZlEMmbqKExeuWR2eECKbkUSeAT6e7jz7aGXctm2jXGBOevZswf9e+ZThP+7i4rWbqb+AEEI4gIfVAbiCvCEPUHXRdI6v3sCD8TdZ/e1iDk27QuXnHqNL9Qfx9nC3OkQhhAuTRO5A+WtXZURt2O93jhwvvsDONXN4utWLdOjyMC0qFsDNTVkdohDCBUlpJRMUbdeM/EfiCenwKG/OHE6v6ZtoPnYVMbtPyR4uQgiHk0SeWby8KPreYELjtzKmUxhvje/F8p6DeXxCLJuPnLc6OiGEC5HSSiZz83CnZVhBbs6YQoFuPbj2ZhT9V/cgX5OG9HmkJMUCXPdSWEII55AZuZN4hlWkwIbV5P90DO3DCrB25zGee/0bBszawvELV60OTwhhYzIjdyalyNG2Ne2BhqvW4tX8Sb7/tS4tYh+jdYPyvFCvGHl9vayOUghhMzIjt0jeOjXwjd9Nu7IBLP/8BWb9uIm6w6P5eNleLl+/ZXV4QggbkURupYAA/KZOJtdvm/hmUAv6HF3NhonfUnd4NFNWH+D6LdllUQiROimtZAWFC1MCKNG2Bh17vszOnUvoc/hxPl9dgv81CKVN5QfwcJe/uUKIu5PskJU0bYrPrh2EPd2Wr65vxN/Pi4EzNtFozEoWbDkm2+YKIe5KZuRZjZcX9O5NAWDOuXNcK1uez6q3438n6/PJA3np80gJIksGopSsEhVCGDIjz8JU3rzkWLKYnhe38evMvgQe3M0zX8TR9tM1rIk/Y3V4QogsQmbkWV358qiffyb3/PlMrlqNRWv38s2qPTx2+Dw1i+Xj1UYlrI5QCGExmZHbgVLQogWe+YNp6XaW6V/1ZeG+Hzhx4HfafrqW0RuvsfXoBaujFEJYRBK53bRqhdqxg7L5vFk6ZzD9GoWy7/xtHh23mu7/F8euExetjlAI4WSSyO0oMBAmTMBt3Vp6RIay+NeJfHjfCdbuO0vjMat48ZtNxJ+SS88JkV1IIrez3LkBuFDjIVpPeZ9Nv3zE68XdiN51ikYfrqTXd5s5cOayxUEKITKbJHK7U4qztWrBtm14NmxAV9/zrOobQc/K/vy47QQNR6/g1Rm/ceisJHQhXJUkclfh7Q2vvgqdO5Nv5xZ692zBL3l28/RDBVmw5Rj1R62g3w+/ceSPK1ZHKoRwMEnkrqhaNVi6lFyL5/P6wM6sfqoMT9QozJzNx4gcGUP/mVskoQvhQhySyJVSjZVSu5VS8Uqp/o54TZFB5cvDzz/DuHEEFHuQN71/J7ZlAbpUL8ysX3+XhC6EC8lwIldKuQPjgSZAGaCzUqpMRl9XOIBSEBlpPh4+TEDThgxZPplV3Srx2EMPMmuTSeiv/bCFw2cloQthV46YkVcD4rXW+7XWN4DpQEsHvK5wpK5dYccOuH6doNd68VbLcqzoF0HUQw8ye/PvRI6Koe/3clJUCDtyxBL9B4Ajd9w/CjzkgNd1IG9gNUqNAhYBq4Eb1oZkITfA77vvWACsA/YmHR+ZdBNCpC46OtrqEP7ktL1WlFLdge4AQUFBxMTEOOutMUn7BUz1ZxiwFugNtAY2AoedGIv1EoGLmKQ9HtgNvAQcsDIoIWwmISHByXnsP2itM3QDagBL7rg/ABjwX19TpUoV7UyANkM1EhPN7amntPb317p0aa0/+ijlMbuJjo5O/xdfu6b1iBFa79mj9ZkzWv/xhz558aoetnCHLvX6Yh3Sf4F+cdpGvev4RYfFe68yNL4sTsZmX1aMD4jTd8mpjqiRbwBClVJFlFJeQCdgngNeN9MoZW5Tp8LJk/DFF1CunHmsTh1o0QI+/RQOHrQySifx9oY+fSA0FBYvhlKlCPxqCgMbhbL6tUier1eM6F2neGTMSp77Ko5tv8vmXEJkNRlO5FrrW0BPYAmwE5ihtd6e0dd1Fjc303Zdv765P2cOdOoEa9ZAv37m2MyZsHQpXL9uWZjO0aWLaVmcORNq1CBfDg9ea1yK2P71eblBKGv2naX5x6t55osNbDp8zupohRBJHFIj11ovwpxFtD1/f3jsMXNLdvIkjBxpmj5atICvvoIbN8zFfFxOhQrmr9beveDuDqNHk6dZM3o/XJKudYrwf2sO8vnqA7T5ZA21iuejZ2Qo1YveL1csEsJCsrIzDXr0gLVrYd8+08WXfKxMGbMqftkySEy0NkaHUgpKlACddI3Q2rWhVy9yX02gZ/1QVr9Wn4FNS7HnZAKdJ6+j3YS1RO86lXyORAjhZJLI74G/P9SrZz6fNAm+/BLuuw9GjTLHFi6ECRPg0CHrYnQopaB3b9i+Ha5eNScTAF8PRfe6xVjVL5K3W5blxIVrPP3FBpp/vJpFW49zWy4SLYRTSSJPJzc3qFoVBg+GRYvMfV9fiI2F8HAzWz9/Hi5edIHaetL+5/TqBRs2mPLLTz/h4+nO4zVCiO4TwYh2Fbh64zY9pm3i4Q9X8H3cEW7edqV/U4TIuiSRO1BEhKmfnzwJ06ZBnjwwY4bJgy1bmlx4we5NH+Hh8O678OKL0Lw5nD6Nl4cb7cML8XPveox7LAwfD3f6/rCFesOj+SL2AFdv3LY6aiFcmiTyTODmBmFh5vOuXU1tvWNHM1u/fh3Wr0+prdtutq6U+au0bRs8+qj5a7VlC5w7h7ubonmFAix8uTZTn65KgTw5GDJ/B7U/WM746HguXL1pdfRCuCRJ5E6Q3Anz1Vdmdh4UZGrrgwaZ++vXw5UrNqute3vDc8+BpyfMnw+lSsH48XDrFkopIksG8sMLNZnxXA3KF7yPEUt2U+v95by3eCenLl6zOnohXIokcguEhJja+rp1ZrZesaJpbaxa1dTW+/SBPXusjvIeDBqU0n/epctfHqpW5H6+eLoaC1+uTWSpQCav3E/tD6IZMGurXIZOCAeRRG4xf3/w8TGl5xMnTGNIrlyQkADHj5sqxsSJcDirbwdToYKpFY0fbzpcunaF3bv/fLhsgfv4uHMY0X0iaBdekJmbjlJ/VAw9pm1ky9Hz1sUthAuQRJ6FJK8yffNNqFzZJPSOHWH1aqhSxWwboDWsXGkWJGU5SkG+fGYhUalSUKuW6XQ5l7IKtHA+X95tXZ7Vr0XyQr1irNp7hhbjYnls8jpW7jktvehCpIMk8izMzy+ltn7yJDz1lGln7NcPAgKgVStYssTqKO/Cy8vUh3bsMMX/3bvh2jW4devPpwTm8qFf41Ks6V+fAU1KEX8qgSem/EKzsauZu/l3bknrohBpJoncJtzcIEcOc5J03TqIj4f27c3kNzHRbPbVp4+pbmSZ2XpgoKkLVa9u/hpVrAg//fSXp+Ty8eS5esVY9VokH7Qtz/Vbt/nf9M1EjIxhauwBrty49S8vLoRIJoncpgICICoKGjY090eNMqWYgQPhhRfMsUWLslAnTNeupv+8Rw9T+L/9195ybw93OlZ9kJ971WPyE+EE5/Zh6Pwd1Hx/OTP33uD0Jbv1aQrhPJLIXcCdtfX162HyZHN83ryUVaYjRphjlpWgk/vPt283ydzdHebO/Uv9HMDNTfFwmSB+eKEmM1+oQbWQ+1mw7ya1PljOgFlb2Hc6waIBCJF1SSJ3QW5J39UJE1I6YcqWNceaNze1dcs6Yby94ZFHzF+UZcvMSdFPPvlL/TxZlcL3M+mJcN6rk4N2VQoyc9PvNBi1gq5fbmD9/rNyYlSIJJLIXZy7u5mtN21q7n/5JXToYDphunUzxxYtsqC2rhSMHZvSf/7BB//61GBfN95tXZ41Sfuibzx0jo6T1tFyfCzzfjsmJ0ZFtieJPJu5c5VpcsfLsWNmTU9AgDmBqvVdJ8iZI3n/8759zbL/Fi3+0n/+l9j9vOn9cAnW9G/AO63KcenaLV7+9lfqjYjhs1X7uXRNtgAQ2ZMkckHXrimrTLt1M5Plvn1NOaZPH1i+PJMTu1KmZTE0FOrWTek/v3jxrk/P4eVOl+qFWda7HpMer8IDeXPwzsKd1HxvOe8s2MHRc1cyMVghsh5J5OJP/v7QqJH5fORIc03TXLng7bdNi2N0dCbX1pOvH7pjhyn0KwXHj6Nu3333RDc3RaOywcx4rgbzetYislQgU9ccpO7waF6ctomNh+RydCJ7kEQu7iq5tv7mmyaBe3mZ/bFWrTKrTMuWNSWZK1cyYQfHwMCUfsoRIwjv2tXU0v9DhYJ5GNs5jFX9IulWtygr956m7adraDU+lvlSRxcuThK5SLPateHrr00nzNSpEBxsWhyT91ufOPEf3YQZN2oUB5591jTHJzfI/4cCeXIwoElp1g1owNAWZTl/5QYvffsrdYdHM2HFPi5ckTq6cD2SyMU9S56tu7lBp05mlWmHDma2/scf5pylw1aZKsWZ2rVN//nLL5uFRMOHp/oXw9fbgydrhrDs1QgmPxFO4Xy+vL94F9XfW8brc7YSf0r60YXrkEQuMix5lenXX0OxYmYbgVy5Ujphli83CT1DtXVvbyhd2uzZsn//f/af38k9aYHRt92rs+jlOjSrkJ8ZG47ScPQKnpzyCzG7T5Eo1xgVNieJXDhcoUKmtp7cCVOtmvmYXFtPPp+ZLr6+ZqXTzz/DDz+YaX8aFwaVKZCbke0rsmZAfXo/XIIdxy/y1NQNNPxwBf+39iCXr8u+LsKeJJGLTOXvb3ZxLF06pbaeKxecOWO6C9u0SWcnTPL+540awZQp5rJz/9J//o+Y/Lx5uUEosa/VZ0zHSuTy9mDw3O1Uf3cZb83fwUG54IWwGUnkwmnu7ISpWxc8PKBt25ROmOT9YNatS2NtXSlz69IF6tUzZ2N7905z07uXhxutwh5gzou1mPlCTSJLBfJ/aw8SOSqGZ77YwIo9p6XsImxBErmwTM6cKbX1EydMU8r16/DKKyn7rc+fn4YXSu4/377dFOk9PODXX9Oc0JVSVCmcl7Gdw4jtX5+X6oey5egFnpzyCw1Hr+CL2AOyalRkaZLIRZbg7m5KMN7ef91vXWtza9o0DZ0wgYHw4ovmC4YOvev+56kJyu1D74dLENs/kjEdK5E7hydD5u+g+rvLeGPONvaevJTxwQrhYJLIRZaU3AnTooW5P2SISfQDB6Zc33n58n+prSsFs2fDsGFmy9zkms098PZw/7PsMvfFWjQul5/v4o7w8Icr6TxpHT9uOy6LjESWIYlcZHlKmdr6kCFmv/VvvjHH585N6YQZNuwuX9SqlSm3PPMMHDlipvTpWLFUsVAeRnWoyNr+9enXuCSH/7jC819vos7waMYu28upS9cyOkQhMkQSubAdDw/z8aOPUjphSpc2x6Ki/rbfure3uSB0jhxw6ZLpP//003TtApbPz5seEcVZ2S+SSY9XoXigH6N/3kPN95bT85tNske6sIyH1QEIkRHJnTDVqpn7Y8aYsvjixaYrcd06iIkBrf2p/fFEvHr0gNdfN9k+b17w8bn390zarKtR2WAOnLnM1+sO8X3cERZsOU6JID+6VC9M67AHyOXj6cihCvGvMjQjV0q1V0ptV0olKqXCHRWUEOl15yrT9etNheXwYRgwwDzWekhFbs6az+3A/KZVsUUL2LMn3e9XxN+XN5qXYf3AhgxvWwFvD3cGz93OQ+8uY8CsrWw/dsGBoxPi7jJaWtkGtAFWOiAWITLFk0+apB4fb8rlnp7w1lsQFvMhCy7W5Ua1Wtx+Y0iG3iOHlzsdqhZiXk9zcrRZ+fzM/vUozcauptX4WL6PO8LVG3ffjleIjMpQItda79Rap205nRAWCwgwC0ABBg+GCVO9iYvoQ+vQ7VyrWodf1muWvDCHw/vTv1RfKUXFQnkY0b4i6wc0ZHDzMiRcv0XfH7bw0LtLGTJvu7QwCodTjjg5o5SKAfporeP+4zndge4AQUFBVaZPn57h902ryMhIAKKjo532ns6UkJCAn5+f1WFkGmeNb88md6q+MwDvC+d41/99mowOICDgOkqBp2f6f0+01uw5l8jywzeJO3mb2xpK5HUjopAnpXyvcf99rvm9k59Lx4uMjNyotf5HGTvVRK6UWgoE3+WhQVrruUnPiSGVRH6n8PBwHReXpqc6hFIKwGU7CmJiYoiIiLA6jEzj1PFpze3Zc7n5vz64jxvLEvemREVBZCQ0aWK2FPD3T//Ln024zg8bj/LNL4c5dPYKvp7QsVoROlcrRGhQLseNIwuQn0vHU0rdNZGn2rWitW6YOSEJkQUphXubVrg3awIeHjSfO4vjnVexsPJg5q7My0MPQUICjB9vEnvt2ubqSWmVz8+b5+oVo1udoqzdf5aPFsTx1bqDTIk9QHjhvHSu9iBNy+cnh5d75o1RuBzpIxfibry9TW9jrVrkTLxM+zdK8XXNT6hU/jZeXmY33eROmAULzDVNjxxJ+8u7uSlqFfenRyUf1g5owIAmpTh7+Qavfv8b1d5dyuC529hx7O4Xnxbi7zLafthaKXUUqAEsVEotcUxYQmQRQUEwaRIsWQK7doFSFPA8/ecq0/h4Mys/cgTCwqBcOejbF7ZuTftb+CfN0pe/Wo9vu1WnQalApm84QtOxq2gxbjXT1h+STbvEf8po18psrXVBrbW31jpIa/2IowITIkupVAnGjjXXt2vXzvSf791LQADkyQOFC8PJk/D552a2fvSo2dyrY0fzdyAts3WlFDWK5WNMpzB+GdiANx8tw/WbiQyavY1qw5bR5/vf2HDwD5c91yPST0orQtyrn36COnWgRg2TpZO4u8NDD5k9YZo0MZswPvoorFhhZutvvmmet2lT6vut58npxdO1ivDjK3WY82ItWoUVYPHW47SfsJYGo1cwYcU+2eNF/EmW6Atxr7y9Tf3kiSfMmc/z52HGDLPayMPjL0/r0sXcbt82T01MhP/9D7ZsMZ0wlSv781+ND0opKhXKQ6VCeXi9WRkWbj3OjA1HeH/xLkYs2U1kyUA6hBckslQgnu4yL8uuJJELkV5BQeZ26BBMnw4ffwwffggN/9no5e5uLkoN5opIp0+bsvuePSb5du4MBQv+dyeMr7cHHcIL0SG8EPtOJ/B93FFmbjrK0p0n8ffzonXYA7QPL0QJF2tjFKmTP+FCZFThwuaKF2+/DS+9ZAriqdSxAwLMTL1+/VMA9OqV0gmTvAd7bOy/19aLBfjRv0kp1vSvz2dPhFOlcF6mxh6k0YcraTluNV+tO8SFK3KCNLuQGbkQjpC8/3mLFuaEaJ8+po4yeLA5G5qK5B0chwyBm0n5d+5cs4NjcLA5vzpkyD+/ztPdjYZlgmhYJogzCdeZu/kY38cd4Y0523h7wQ4alQmiXZWC1AkNwN1NOXDAIiuRGbkQjuSW9CvVty9cvgwlS8Jnn93TS3gm7X47fHhKJ0zx4ubYiy9C69Z374Tx9/Pm2dpFWPy/Oix4qTaPVXuQ2PgzPDV1AzXeW8Z7i3ayR/Z5cUmSyIXIDEFB5uoWP/1kCuRgrlZ0j5I7YZIvbzdkiNkmYMUKaNTITPrXroXo6JROGKUU5R64jyEtyrJ+YEMmdKlMhYJ5+Hz1ARp9uJJHP17NF7EHOJtw3TFjFZaT0ooQmaliRXNLSDCll9KlYeRIKFEiXS+XXFvv0sWU4ZWCgwfNOdY9e0wnzFdfQc6c5p8DLw83GpfLT+Ny+f8svczadJQh83fwzsKdRJQMpG3lB6hfOhBvD9kWwK4kkQvhDH5+sG2bWVRUs6a58GijRhl6yaS94Ojc2dxOnTIzdV9fGD3aXAKvSRNzq1MnpfTybO0i7Dpxkdmbfmf2r7+zdOdJcvt40LxiAdqEPUCVwnn/3GhO2IMkciGc5c7+81y5IDaWAnPnmn5Dj4z/KgYGQvv25vNXXjEvu2gR9O8Pc+bAH3/AmjUmsZcqlJsBTXPTr3EpYuPPMGvTUWZtOso36w9T6P4ctK70AC3DHqBYgOtuQ+tKpEYuhLMFBZnaR548BMTEmGWfS5c69C2Sa+tDh8Ivv0CBAqYbJnmVably5h+ExNuK6iEBjOkURtzrDzOqfUVC8vnycXQ8DUatoOW41UyNPcDpS1JPz8pkRi6EVcqW5bfRo4k4fx5GjIB69cxxz8y5aHNYGEybZlaZxsVBSIjpVW/VKnm/dQ9aty5I2yoFOXHhGvN/O8aczb8zNKmeXqu4P60qFaBR2WD8vCV1ZCXy3RDCSkqZfsLWrc39evUgPBzeeCNN/efpkTxbB4iIgL17zSrTxYuhTBlzfPQIH5o0Kcrs54ty6Nwl5mz+nbmbj9F7xm/4eG6lQekgWlYsQL2SAXKSNAuQRC5EVjJjhkniJUuaq1e0a5fpb3lnJwyYk6Y5csBrr5lOmIkTc9G3YykeK1eS4zfPMXfzMRZuPc7CLcfJ7eNB0/L5aVGxAA8VzSeLjiwiiVyIrCR5//MePeD6dXNbs8bUPpwkMNDU1ocONUldKThzBipWVOTPfz9NmtzPxI5lSPA9w7zfjjHvt2NM33CEwFzeNKuQn0crFiCsUB6nxSskkQuRNVWqZD7u2QPdupmax6hREBrq1DACA1M+P3kSNmwwJZiD+91o2zaQL94LpMdDifgWP8W6k0eZtv4wU2MPUjBvDirmuUVgiYuUzp9L2hkzmSRyIbKyEiXMitDk/vO4OLNJlwXc3aF6dXMD0wXz8MOweLEbSwYH8+STwcQNu8nkOWf59eIRFh84xcKxqyga4Evz8vlpXrGA7MyYSSSRC5HVJfefP/cc5M4NY8aYY926OaT/PL08Pf+63/r585DL25PFnwUTFxdM6bKnCW8KZ3Lt4+PoeMYuj6dEkB/NKxSgWYX80qPuQNJHLoRd5M5tPkZEwHffZUr/eXq5u0O+fKaevnSpqQjVjzhDWP4Avu1enTonH6HMkdrcOOzP6CV7aTBqBY3HrOTjZXvZdzrB6vBtTxK5EHZTqZLZJeuttyAmxhw7f97CgP4pMBAefvgkzzxj7vd83oOqxe/jSmxZAlc2ZnDzMtw8cR8fzDz8Z1Ifu2wv8ackqaeHlFaEsKM7+88vX4ayZaFTp0ztP8+I5P3Whw6F69fd8PYuwqGf4dfvNDnz3ORg6ClGHf+N0T/voUSQH03K5adp+fyUCPKTE6VpIDNyIezO19dc0fniRShVClautDqi/+TtbT4OHQqnTyl+mObFq20Lsn5gA8odq8X2L8szbPQ1Grz1Cw1GrWD4j7vYevQCOpWrLmVnMiMXwhUEBcHkybB5s7n45+7d5soTd7l+aFby104YH6aM9GHJEpg97z5WztIEvrmRsTPOMPrL8xQtd5WmFYNoXC6Yyg/mxU0WH/1JErkQriS5/3zXLtPlUq6c2f/cyf3n6RUYCI8/Do8/7p603/pDfJXnJkPfSWTjXHe2FDrDxEfiCAxw45FyQTQuG0z1ovnw8sjexYXsPXohXFXt2rBjB9SqBS1bplwI1EaSS+OPd/Ykfrs3Rw56MKa/P2OfKkeuQ6GMfi6ERx+/SMlnN9Lz680s2nqcy9dvWRu0RSSRC+GqvL2hXz/YssU0fXfoABMmwC17JrvAQOj6tAetKxdg2cTCLJmVkxZVgrgWW5Zl6xPoNm4XRdvspOPozUz/5XC22npXSitCuLrkRUMDB5orTowfD+PGpWyba0Pu7lCvtjv1aptFRbdu1+S7ny7y/vsw642czM55Df9mG6he1Y36JYNoXCGI4oGuuwBJErkQ2UVy//mcOabvPDERDh2CIkUsDizjPNzdiGqSh6gmcOuW5vvFNzlwMz+L113g5VcfxKfwWR6s8DttWila1fSnSuG8LrVTo5RWhMhOkvvPW7aEnTuhalV49dUst6AoIzw8FJ0f9WVgm2KsGl6ZzVtv8XhHd64f8eeLeRdo9/Ev5K97kHaD45mz8RiXrtnv/MHfSSIXIrsqW9Zc7y25//zgQasjyhQVQnMw8a0ADqzLx+4vKzG8dSWKBuVg0ecBtKkVQPF2O4n6bB1j5h/i0NnLVoebLhlK5EqpEUqpXUqpLUqp2UqpPA6KSwjhDMHBpv88JsbsqjhzJixbZnVUmSaXjyed6gSzbmYwl47kZtGqBJ7plIPjZ27yaqdgQkvd5sHIw/T4aD9r9p3h5u1Eq0NOk4zOyH8GymmtKwB7gAEZD0kI4XSlSpmyi4+P2VWxVStzDTgX5u6maFw1L+91CWV5/zrEH7jNa+8kkMfPjZk/J/DY5PUEVz9KRNdDTFh4jFOXrpkvnDYNQkKoV7++ufDptGmWjgMymMi11j9prZN7mdYBBTMekhDCMs2amf7zmjVhyBBzLNEes9KMKhKYk7efK8CW+QXZP7MME6KqUKMGbN3gzYvt/SlW9yQjnhjMtqc/4Oah31Fam5PF3btbnswd2bXyDPCdA19PCGEFHx/Tfw5w5QpUrgy9esGzz1q6/7kz+Xp70Lh8MI3HgtaarUcvsmhDIk92m8SzNyfjzm3mkHTB7CtXuD1gAO5RUZbFq1LbiEYptRQIvstDg7TWc5OeMwgIB9rof3lBpVR3oDtAUFBQlenTp2ck7nsSmXS9w+joaKe9pzMlJCTg5+e6PbKuPD47jM0vPp7i48bhcekSu157jYQSJdL0dXYY272qV78+Smtu4oEnKQurElFEDltIuXzulA9wp3geNzwyob0xMjJyo9Y6/O/HU03kqVFKPQU8BzTQWl9Jy9eEh4fruLi4DL3vvUjeBtNVd0+LiYkhIiLC6jAyjSuPzzZj0xpmzza96DlzwqVLqe7fYpux3YuQEFNO+ZuLgQXoOvQHNh0+x61Eja+XOzWK+VOvhD91QgMI8fd1yNsrpe6ayDP0f5JSqjHQD6iX1iQuhLAhpaBNG/P5okXwxBPw1FPw+utZcv/zTDNsmKmJX7kj3eXMSe7Rw5kRVYNL126yZt9ZVu45zcq9p1m68yQAhe7PQZ3QAOqG+lOruD+5fDwdGlZGC17jAG/g56RZ7zqt9fMZjkoIkXU1bWr6z994A6pUMQuLvLysjso5kuvggwahDx9GPfigSe5Jx3P5ePJI2WAeKRuM1ppDZ6+wcu9pVu09w7zNx/hm/WE+eyKchmWCHBpWhhK51rq4owIRQthIcv/5mTMmiQ8cCA0amJuri4qCqChWpFI6UkoR4u9LiL8vT9QI4ebtRDYfOU/ZArkdHpKs7BRCpJ+/v/kYHm76z1u2dPn+8/TydHejasj95PRyfOePJHIhRMa1aZPSf75kCQBuV69aHFT2IYlcCOEYPj7w2mvQsyfs2EH1qCiYOBFu37Y6MpcniVwI4XhlyrDlgw/gm28gLAwOH7Y6IpeWPZZpCSGcLiE01GzGtWgRFCgAa9ZAQIBtrh9qJzIjF0JkHqXM/i0eHuaC0DVqQN++cOGC1ZG5FEnkQgjneOYZ039+/jx07mx1NC5FErkQwnmS+89nz4arV83iIhfe/9xZJJELIZzP29t0uXTtmrL/+f79VkdlW5LIhRDWSN6/Jbn//OxZSEhwqeuHOoskciGEtZL3P69aFX780VytaOJEuHUr9a8VgCRyIURW0q6daVf85huT2G/csDoiW5BELoTIWipXNv3nU6aYDbmmTIH4eKujytIkkQshsh6lzIpQgD/+gOrVpf/8P0giF0JkbX36mP7zc+dg5EhzzEWv9pVeksiFEFlfcDB89hm89ZZZIVq5MixfbnVUWYYkciGEfSgFJUuaS8x17Wr6z0+csDoqy0kiF0LYi1LQtq3pP69XD/z8zMnQbFw/l0QuhLAnHx/o1csk8lmzzEw9m+5/LolcCGF//frB4sWm/7xNG6ujcTpJ5EII1xAWZvrPP/3UzMp79co21w+VRC6EcB1KmYtY3L4NQUHZZv9zSeRCCNfj5QX9+6f0n8fGmr1bXLR+LolcCOG6kvvPmzaF7783/efR0VZH5XByzU4hRPbQqZOZqT/7LFSsCNOnm33RXYDMyIUQ2cOd/edRUSaJL1vmEvVzSeRCiOzFx8dslwswb57Z/3zSJFvXzyWRCyGyr48+goULYdo0GDTI6mjSTRK5ECJ7S97/fPBgOHQIOnSw3f7nksiFEEIpyJnT9J5Xrpyy//nFi1ZHliaSyIUQIpmPT0r/+eXLcPWqubBFFq+fZyiRK6XeVkptUUptVkr9pJQq4KjAhBDCMsHB8MknZoY+cmSW3/88ozPyEVrrClrrSsACYHDGQxJCiCxk2DB44w3Tf/7MM1ZHc1cZSuRa6zsLSL6AXH9JCOFalDLtijt3Qo8e5jJzX35pdVR/oXQGr32nlBoGPAFcACK11qf/5Xndge4AQUFBVaZPn56h970XkZGRAES74NJcgISEBPz8/KwOI9O48vhkbPbjduMGeePiOFShgtPHFxkZuVFrHf7346kmcqXUUiD4Lg8N0lrPveN5AwAfrfWbqQUTHh6u4+LiUo/aQZRSAGT0j1ZWFRMTQ0REhNVhZBpXHp+Mzb6sGJ9S6q6JPNW9VrTWDdP4HtOARUCqiVwIIYTjZLRrJfSOuy2BXRkLRwghxL3K6O6H7yulSgKJwCHg+YyHJIQQ4l5kKJFrrds6KhAhhBDpIys7hRDC5iSRCyGEzUkiF0IIm5NELoQQNieJXAghbE4SuRBC2JwkciGEsDlJ5EIIYXOSyIUQwuYkkQshhM1JIhdCCJuTRC6EEDYniVwIIWxOErkQQthcRvcjtwWtNTExMVaHIYQQmUJm5EIIYXOSyIUQwuYkkQshhM1JIhdCCJuTRC6EEDYniVwIIWxOErkQQticJHIhhLA5SeRCCGFzSmvt/DdV6jRwyMlv6w+ccfJ7Oosrjw1ce3wyNvuyYnyFtdYBfz9oSSK3glIqTmsdbnUcmcGVxwauPT4Zm31lpfFJaUUIIWxOErkQQthcdkrkk6wOIBO58tjAtccnY7OvLDO+bFMjF0IIV5WdZuRCCOGSJJELIYTNuVwiV0o1VkrtVkrFK6X63+Vxb6XUd0mPr1dKhVgQZrqkYWy9lVI7lFJblFLLlFKFrYgzPVIb2x3Pa6uU0kqpLNH2lVZpGZ9SqkPS92+7UuobZ8eYXmn4uXxQKRWtlPo16WezqRVxpodSaopS6pRSatu/PK6UUmOTxr5FKVXZ2TEC5jJornID3IF9QFHAC/gNKPO35/QAJiR93gn4zuq4HTi2SCBn0ucvuNLYkp6XC1gJrAPCrY7bwd+7UOBXIG/S/UCr43bg2CYBLyR9XgY4aHXc9zC+ukBlYNu/PN4UWAwooDqw3oo4XW1GXg2I11rv11rfAKYDLf/2nJbAl0mf/wA0UEopJ8aYXqmOTWsdrbW+knR3HVDQyTGmV1q+bwBvAx8A15wZnAOkZXzdgPFa63MAWutTTo4xvdIyNg3kTvr8PuCYE+PLEK31SuCP/3hKS+D/tLEOyKOUyu+c6FK4WiJ/ADhyx/2jScfu+hyt9S3gApDPKdFlTFrGdqdnMTMFO0h1bEn/shbSWi90ZmAOkpbvXQmghFIqVim1TinV2GnRZUxaxjYE6KKUOgosAl5yTmhOca+/l5nCw9lvKDKfUqoLEA7UszoWR1BKuQGjgacsDiUzeWDKKxGY/6RWKqXKa63PWxmUg3QGvtBaj1JK1QC+UkqV01onWh2Yq3C1GfnvQKE77hdMOnbX5yilPDD/6p11SnQZk5axoZRqCAwCWmitrzsptoxKbWy5gHJAjFLqIKYWOc9GJzzT8r07CszTWt/UWh8A9mASe1aXlrE9C8wA0FqvBXwwG065gjT9XmY2V0vkG4BQpVQRpZQX5mTmvL89Zx7wZNLn7YDlOumsRRaX6tiUUmHAREwSt0uNFVIZm9b6gtbaX2sdorUOwdT/W2it46wJ956l5edyDmY2jlLKH1Nq2e/EGNMrLWM7DDQAUEqVxiTy006NMvPMA55I6l6pDlzQWh93ehRWnxXOhLPMTTGzmX3AoKRjb2F+8cH8EH0PxAO/AEWtjtmBY1sKnAQ2J93mWR2zo8b2t+fGYKOulTR+7xSmfLQD2Ap0sjpmB46tDBCL6WjZDDSyOuZ7GNu3wHHgJua/pmeB54Hn7/i+jU8a+1arfi5lib4QQticq5VWhBAi25FELoQQNieJXAghbE4SuRBC2JwkciGEsDlJ5EIIYXOSyIUQwub+HxgQTRAKQRc7AAAAAElFTkSuQmCC", "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": 7, "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": 8, "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": 9, "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}}}+ \\\\\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": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApAElEQVR4nO3dd3hUZd7/8feXEpAi0ha5BAEFV1xdSkISmkS6goCNFSLig8oCuoq7KFVEFB4VC/4QFRZcLCBFpcoCIRAQJSGERRBFKRoW7BQ1AqHdvz8y+EQMApnJnCmf13XNlVPue+ZzezDfzNznzDHnHCIiEr2KeR1ARES8pUIgIhLlVAhERKKcCoGISJRTIRARiXIlvA5QGFWqVHG1a9cuVN+ff/6ZsmXLBjZQiNOYo4PGHPn8HW9WVtb3zrmqp24Py0JQu3Zt1q9fX6i+aWlpJCUlBTZQiNOYo4PGHPn8Ha+ZZRe0XR8NiYhEORUCEZEoF5BCYGYdzexTM9tuZkMK2H+1mW0ws2NmdvMp+3qb2Tbfo3cg8oiIyNnzuxCYWXFgInAtcAXQw8yuOKXZLuAOYMYpfSsBjwAJQDzwiJlV9DeTiIicvUC8I4gHtjvndjrnjgAzga75GzjnvnDObQJOnNK3A5DinNvnnNsPpAAdA5BJRETOUiDOGroI+G++9d3k/YVf2L4XFdTQzPoCfQGqVatGWlraOQcFyMnJKXTfcKUxRweNOfIV1XjD5vRR59xkYDJAXFycK+wpVNF2uhlozNFCY458RTXeQHw0tAeomW+9hm9bUfc9Z3PnziUlJaWonl5EJCwF4h1BJlDPzOqQ90v8VqDnWfZdCozNN0HcHhgagEy/4ZxjypQpLF68mJiYGEaOHImZFcVLiYiEFb/fETjnjgH3kvdL/RNgtnNui5mNNrMuAGbWxMx2A7cAk8xsi6/vPuAx8opJJjDaty3gzIy5c+fSoUMHRo0aRe/evcnNzS2KlxIRCSsBmSNwzi0GFp+ybWS+5UzyPvYpqO8rwCuByHEmMTExDB48mBYtWvDwww+TnZ3N3LlzqVSpUjBeXkQkJEXdlcVmxogRI5gxYwbp6ek0bdqU7du3ex1LRMQzUVcITurRowepqans3buXxMRE1qxZ43UkERFPRG0hAGjRogXp6elUrlyZNm3aMGPGjDN3EhGJMFFdCADq1q3L2rVrSUxMJDk5mdGjR+Oc8zqWiEjQRH0hAKhUqRLLli3j9ttv55FHHtEZRSISVcLmyuKiVqpUKaZNm0a9evV0RpGIRBW9I8gn/xlFGRkZJCYmsm3bNq9jiYgUKRWCApw8o2jfvn0kJiayevVqryOJiBQZFYLTaN68ORkZGVStWpW2bdvy+uuvex1JRKRIqBD8jksvvZS1a9fSvHnzXyaSdUaRiEQaFYIzqFixIkuXLuV//ud/GD16NMnJyRw+fNjrWCIiAaOzhs5CTEwMU6dO5bLLLmPo0KFkZ2czb948qlat6nU0ERG/6R3BWTIzhgwZwuzZs9mwYQOJiYls3brV61giIn5TIThHt9xyC2lpafz8888kJiaSmprqdSQREb+oEBRCQkICGRkZ1KxZk44dOzJlyhSvI0mI2LRpE/369ePKK6+kfPnyFC9enPLly/PnP/+ZBx98kO+++87riCK/oUJQSLVq1WLNmjW0adOGu+++m8GDB3PixAmvY4mHXn31VWJjY5k0aRJbtmzBOUfFihU5evQomzdv5umnn+ajjz7yOqbIb6gQ+KFChQosWrSI/v3789RTT3HzzTfz888/ex1LPPDDDz9wzz33cOzYMR566CG+/vprcnJy+P777zl06BA7duxgzJgxNGzY0OuoIr+hQuCnEiVKMHHiRMaPH8+8efNo1aoVX375pdexJMjWr1/Pzz//TLFixRg2bBjVqlX7ZZ+ZcckllzBs2DAqVqz4O88i4g0VggAwM+6//34WLFjA1q1bSUhI4MMPP/Q6lgTRpZdeSokSJThx4gTx8fGMHz+eTZs26QJECQsqBAHUuXPnX+501rx5cxYtWuRxIgmW2rVrM3XqVCDvupNnn32WBg0aUKVKFQYMGMDXX3/tcUKR01MhCLCGDRuybt066tevT9euXRk/frz+KowCGRkZjB49mgkTJrB582Z27drFl19+yY033shLL71EgwYN+PTTT72OKVIgFYIiUL16dVatWkW3bt144IEHGDBgAEePHvU6lhSRjIwM2rRpQ+/evbn33nt/2V69enX++c9/0qVLF7799lvuu+8+D1OKnJ4KQREpU6YMc+bMYfDgwbz88st06tSJAwcOeB1LAuzAgQPcfPPNVKpUicGDBxfY5s477wRgxYoVHDp0KJjxRM6KCkERKlasGE888QRTp05l5cqVNGvWjJ07d3odSwJo3Lhx7N69m969exMTE1Ngm5PfSXXs2DH9MSAhSYUgCPr06UNKSgpff/01CQkJvP/++15HkgCZNm0aAG3atDltm++//x7IO7tMp49KKFIhCJKkpCQyMjKoWLEirVu35o033vA6kvjp888//+Wakauuuuq07T744AMg78yi0qVLByWbyLlQIQiievXqkZ6eTtOmTenVqxcPP/ywvpYijOX/mO/3/tKfP38+AG3bti3yTCKFoUIQZJUqVWLZsmX06dOHxx9/nB49emgCMUwdP378l+W9e/cW2GbhwoV88sknQN5HhCKhKCCFwMw6mtmnZrbdzIYUsL+Umc3y7c8ws9q+7bXN7JCZbfQ9Xg5EnlAXExPDlClTeOqpp5gzZw5JSUm64CgM1alT55flk3/15/f999/Tv39/IO9iw8TExKBlEzkXfhcCMysOTASuBa4AepjZFac0uxPY75yrCzwHPJlv3w7nXEPfo5+/ecKFmfHggw/yzjvv8NFHHxEfH8+mTZu8jiXnoF69ejRo0ACAwYMHs2zZMgCccyxfvpzmzZuzZ88eatWqxeTJk72MKvK7AvGOIB7Y7pzb6Zw7AswEup7Spivwqm/5LaCNmVkAXjvsdevWjffee4/jx4/raynC0OTJkylTpgz79u2jQ4cOlC9fnrJly9KuXTs+++wzGjZsyIoVK6hevbrXUUVOy/z9+gMzuxno6Jy7y7feC0hwzt2br81Hvja7fes7gASgHLAF+Az4ERjhnHvvNK/TF+gLUK1atdiZM2cWKm9OTg7lypUrVN+i9N133zF8+HC2b99O//79ufnmmwlUrQzVMRelYI55165dvP766/znP//hhx9+oEyZMtStW5c2bdrQoUMHihcvHpQcOs6Rz9/xXnPNNVnOubjf7HDO+fUAbgam5FvvBbxwSpuPgBr51ncAVYBSQGXftljgv8D5Z3rN2NhYV1grV64sdN+ilpOT42688UYHuL59+7ojR44E5HlDecxFRWOODtE2Zn/HC6x3BfxODcRHQ3uAmvnWa/i2FdjGzEoAFYC9zrlc59xeX0HK8hWIywKQKSyVLVuWOXPmMHToUCZPnsy1117L/v37vY4lIhEuEIUgE6hnZnXMLAa4FVhwSpsFQG/f8s3ACuecM7OqvslmzOwSoB4Q1d/BUKxYMcaOHcu0adNYvXo1iYmJbNu2zetYIhLB/C4EzrljwL3AUuATYLZzbouZjTazLr5mU4HKZrYd+Dtw8hTTq4FNZraRvEnkfs65ff5migS9e/cmNTWVvXv3kpCQQFpamteRRCRClQjEkzjnFgOLT9k2Mt/yYeCWAvq9DbwdiAyRqGXLlmRkZHD99dfTrl07Xn755V++yVJEJFB0ZXGIu/TSS1m7di2tW7fmrrvuYtCgQb+6olVExF8qBGGgQoUKvPvuu9xzzz0888wz3HDDDfz0009exxKRCKFCECZKlCjBCy+8wAsvvMDixYtp0aIF2dnZXscSkQigQhBm7rnnHhYvXkx2djbx8fGkp6d7HUlEwpwKQRhq3749a9eupXz58iQlJTFjxgyvI4lIGFMhCFP169cnIyODhIQEkpOTGTlypO5tICKFokIQxipXrkxKSgp9+vThscce49Zbb+XgwYNexxKRMKNCEOZO3tvg6aef5q233qJVq1a/3D5RRORsqBBEADPjH//4B/Pnz2fr1q00adKErKwsr2OJSJhQIYgg119/Pe+//z4lSpTIuyr5/vuhdm1atW4NtWvD9OleRxSREBSQr5iQ0PHnP/+ZdevW8WKLFlz5//4fAAaQnQ19++Y1Sk72LJ+IhB4VgghUrVo1HsnNpRgwkOfYSMO8HQeBO0vDPz0MFyQHDjTkggu8ThFc0Tbm7dvhvPNi0Zfz+k+FIEIV27274B25h4MbRKSI5ORAbm6M1zEiggpBpLr4YsjOZjwP/GrzkerViUmL/LOK0tI2kpSU5HWMoIq2MSclwYEDh8i70aH4Q5PFkWrMGChT5lebDprx1++/55133vEolIiEIhWCSJWcDJMnQ61aODOoVYvcCRP4pHFjbrrpJsaOHXvy/tEiEuVUCCJZcjJ88QWrVqyAL76g4j33kJaWRs+ePRk+fDi9evXi8GHNGYhEOxWCKFO6dGneeOMNxowZw/Tp00lKSuLrr7/2OpaIeEiFIAqZGcOGDePtt99m8+bNxMfHs3HjRq9jiYhHVAii2I033siaNWtwztG8eXNNIotEKRWCKNeoUSMyMzO56qqruOmmmxgzZowmkUWijAqBcOGFF5KWlkZycjIjRowgOTmZQ4cOeR1LRIJEhUCAvEnk119/nbFjx/Lmm2+SlJTEV1995XUsEQkCFQL5hZkxdOhQ5s6dy5YtW2jSpAkbNmzwOpaIFDEVAvmNbt268f7771O8eHFatGjBnDlzvI4kIkVIhUAK1KBBAzIzM2ncuDHdu3fn0Ucf1T2RRSKUCoGc1h/+8AdSU1O54447GDVqlO6JLBKhAlIIzKyjmX1qZtvNbEgB+0uZ2Szf/gwzq51v31Df9k/NrEMg8kjglCpVildeeeWXeyK3bNmS3af7imsRCUt+FwIzKw5MBK4FrgB6mNkVpzS7E9jvnKsLPAc86et7BXAr8CegI/Ci7/kkhJy8J/KiRYvYtm0bcXFxpKenex1LRAIkEPcjiAe2O+d2ApjZTKAr8HG+Nl2BUb7lt4AXzMx822c653KBz81su+/51gYg128MHDiQtLQ0Loim2zgBBw4cCNiYL7/8cjZv3kyzZs344x//SLVq1QLyvIEWyDGHi2gb88aN4zl27BhJSYO8jhI0VapUKZJ7TgSiEFwE/Dff+m4g4XRtnHPHzOwHoLJve/opfS8q6EXMrC/QF/JuxZiWlnbOQXfv3s3x48c5cODAOfcNZ4Ee86WXXsoXX3zB1q1b2bdvHxdeeCF5dT106DhHvmPHjuGci6oxV6hQoVC/+84kbO5Q5pybDEwGiIuLc4WpiklJSaSlpUXVXZyAIhnz0aNH+dvf/sakSZNISEhg+vTplC9fPqCv4Q8d58iXd4eyA1H1hYlFdYwDMVm8B6iZb72Gb1uBbcysBFAB2HuWfSUElSxZkpdffpmJEyeyePFimjZtys6dO72OJSKFEIhCkAnUM7M6ZhZD3uTvglPaLAB6+5ZvBla4vG82WwDc6jurqA5QD1gXgEwSJAMGDGDp0qV8+eWXxMfHF8nbVhEpWn4XAufcMeBeYCnwCTDbObfFzEabWRdfs6lAZd9k8N+BIb6+W4DZ5E0sLwHucc4d9zeTBFebNm1Yt24dVatWpV27dkyaNMnrSCJyDgIyR+CcWwwsPmXbyHzLh4FbTtN3DDAmEDnEO3Xr1iU9PZ0ePXrQr18/Nm/ezHPPPUfJkiW9jiYiZ6AriyVgKlSowMKFCxk0aBATJ06kY8eO7N271+tYInIGKgQSUMWLF2fcuHFMmzaNNWvWkJCQwMcff3zmjiLiGRUCKRK9e/dm1apV5OTkkJiYyKJFi7yOJCKnoUIgRSYxMZH169dTr149unTpwpNPPqnbYIqEIBUCKVI1atTgvffeo3v37gwZMoRevXrpNpgiIUaFQIpcmTJlePPNN3n88ceZPn06rVq14ssvv/Q6loj4qBBIUJgZw4cPZ968eXz88cfExcWxbp2uHRQJBSoEElRdu3Zl7dq1lCpViquvvpo33njD60giUU+FQILuqquuIjMzk8TERHr16sVDDz3E8eO6oFzEKyoE4okqVaqQkpJC//79GTduHF26dOGHH37wOpZIVFIhEM+ULFmSF198kZdeeolly5aRmJjIZ5995nUskaijQiCe69evH8uXL+e7774jISGBZcuWeR1JJKqoEEhIaNWqFZmZmdSsWZNrr72W5557ThefiQSJCoGEjDp16vDBBx/QrVs3/v73v9OnTx9yc3O9jiUS8VQIJKSUK1eOOXPm8MgjjzBt2jSSkpL46quvvI4lEtFUCCTkFCtWjFGjRvHWW2+xadMm4uLiyMzM9DqWSMRSIZCQddNNN/HBBx9QsmRJWrZsyfTp072OJBKRVAgkpDVo0OCXi89uu+02XXwmUgRUCCTkVa1alZSUFAYMGMC4cePo3LkzBw4c8DqWSMRQIZCwULJkSSZOnMikSZNYvnw5CQkJbN261etYIhFBhUDCSt++fVmxYgX79+8nISGBd9991+tIImFPhUDCTsuWLVm/fj2XXnop119/PU888YQuPhPxgwqBhKWLL76YNWvW0L17d4YOHUrPnj05ePCg17FEwpIKgYStk3c++9///V9mzZpFixYt2LVrl9exRMKOCoGENTNjyJAhLFy4kB07dhAXF8d7773ndSyRsKJCIBGhU6dOZGRkULFiRVq3bs2CBQu8jiQSNlQIJGJcfvnlZGRk0K5dO5577jn69evHkSNHvI4lEvL8KgRmVsnMUsxsm+9nxdO06+1rs83MeufbnmZmn5rZRt/jD/7kEbngggtYuHAhPXr0YNKkSbRp04Zvv/3W61giIc3fdwRDgFTnXD0g1bf+K2ZWCXgESADigUdOKRjJzrmGvof+jxW/FS9enL59+zJjxgyysrKIi4tjw4YNXscSCVn+FoKuwKu+5VeBbgW06QCkOOf2Oef2AylARz9fV+SMevTowZo1awBo3rw5M2bM8DiRSGgyfy7EMbMDzrkLfMsG7D+5nq/NIKC0c+5x3/rDwCHn3NNmlgZUBo4DbwOPu9MEMrO+QF+AatWqxc6cObNQmXNycihXrlyh+oaraB/z/v37GTVqFJs2beIvf/kLd999N8WLF/c4YeBF23EeOLAhx48fZ8KEzV5HCRp/j/E111yT5ZyL+80O59zvPoDlwEcFPLoCB05pu7+A/oOAEfnWHwYG+ZYv8v0sDywDbj9THuccsbGxrrBWrlxZ6L7hSmN2Ljc31w0YMMABrkOHDm7fvn3eBCtC0XacW7VyrkGD/V7HCCp/jzGw3hXwO/WMHw0559o6564s4DEf+MbMqgP4fhb0Gf8eoGa+9Rq+bTjnTv78CZhB3hyCSMDFxMQwceJEJk+ezIoVK4iPj2fLli1ex5IQkJ2dzfDhw2nUqBFVqlShdOnS1KxZk6SkJAYNGsShQ4e8jljk/J0jWACcPAuoNzC/gDZLgfZmVtE3SdweWGpmJcysCoCZlQQ6k/dOQ6TI3H333aSlpfHTTz+RmJjIvHnzvI4kHnrnnXeoX78+Y8eOZePGjRw+fJjSpUuzZ88eVq1axdSpUznvvPO8jlnk/C0ETwDtzGwb0Na3jpnFmdkUAOfcPuAxINP3GO3bVoq8grAJ2Ejeu4R/+plH5IyaNWtGVlYW9evX54YbbmDUqFGcOHHC61gSZHv27OG2227j8OHDDBs2jF27dpGTk8OBAwc4cuQIGzZsYNq0aV7HDIoS/nR2zu0F2hSwfT1wV771V4BXTmnzMxDrz+uLFNZFF13E6tWr6devH48++igffvghr732GuXLl/c6mgTJokWLOHToEF27dmXMmDG/2leiRAkaNWpEo0aNPEoXXLqyWKJW6dKl+de//sX48eNZuHAhiYmJbNu2zetYEiTHjh0DICMjg9WrV3ucxlsqBBLVzIz777+fZcuW8c0339CkSROWLFnidSwJgltvvZXLLruMr7/+mlatWlG6dGkuvPBC6tat63W0oFMhEAFat25NZmYmtWvX5rrrruPJJ5/UzW4iXOXKlZk9ezYNGjQAIDc3l2+++YZKlSp5nCz4VAhEfOrUqcP7779P9+7dGTJkCD169ODnn3/2OpYUAeccgwcPJjY2lquuuop169bx448/4pxj3bp1XscLOr8mi0UiTdmyZXnzzTdp3LgxQ4YMYevWrcydO5c6dep4HU0C6Omnn+app57i3nvvZcKECV7H8ZzeEYicwsx46KGHWLx4MdnZ2cTFxZGamup1LAmg559/HoB+/fp5nCQ0qBCInEbHjh3JzMykevXqtG/fnmeffVbzBhHg8OHD7NmzB4DvvvvO4zShQYVA5HfUrVuXtWvX0q1bN/7xj3/Qq1cvDh486HUs8UPp0qWpXr06AH/9619JSUnh6NGjQN6E8datWxkzZgxz5szxMmZQqRCInEH58uWZM2cOjz32GDNmzKBFixZkZ2d7HUv88PjjjwPw2Wef0b59e0qXLk2lSpU477zzqF+/PiNGjKBChQoepwweFQKRs1CsWDFGjBjBggUL2LFjB3FxcaxcudLrWFJIffr0ISUlhRtvvJEaNWpQokQJcnNzqVWrFp06dWLcuHE0a9bM65hBo7OGRM5B586dWbduHd26daNdu3Y888wz3HfffeTdjkPCSdu2bWnbtq3XMUKC3hGInKM//vGPZGRk0LlzZwYOHMgdd9wRFV9VLJFLhUCkEM4//3zeeecdRo0axWuvvUbLli3ZtWuX17FECkWFQKSQihUrxiOPPML8+fP57LPPiIuLY9WqVV7HEjlnKgQifurSpQvr1q2jUqVKtG3blgkTJuh6AwkrKgQiAXD55ZeTkZHBddddx3333ad5AwkrKgQiAVKhQgXmzp2reQMJOyoEIgGUf95g27ZtxMbGkpaW5nUskd+lQiBSBE7OG1SpUoW2bdsyfvx4zRtIyFIhECkiJ683uP7663nggQe4/fbb9T1FEpJUCESK0Pnnn8/bb7/NY489xvTp02nRogVffPGF17FEfkWFQKSInfyeooULF7Jz507d30BCjgqBSJB06tSJzMxMqlWrRvv27XnmmWc0byAhQYVAJIjq1atHRkYGN954I4MGDaJnz566L7J4ToVAJMjKlSvH7NmzeeKJJ5g1axZNmzZlx44dXseSKKZCIOIBM2Pw4MEsWbKE3bt3ExcXx5IlS7yOJVFKhUDEQ+3bt2f9+vXUqlWL6667jrFjx2reQIJOhUDEY5dccgkffPABPXv2ZPjw4dx00038+OOPXseSKOJXITCzSmaWYmbbfD8rnqbdEjM7YGaLTtlex8wyzGy7mc0ysxh/8oiEqzJlyvD666/z3HPPsWDBAhISEti6davXsSRK+PuOYAiQ6pyrB6T61gsyDuhVwPYngeecc3WB/cCdfuYRCVtmxsCBA1m+fDl79+4lPj6eefPmeR1LooC/haAr8Kpv+VWgW0GNnHOpwE/5t1neTV5bA2+dqb9INElKSiIrK4vLL7+cG264geHDh3P8+HGvY0kEM38mpszsgHPuAt+yAftPrhfQNgkY5Jzr7FuvAqT73g1gZjWBfzvnrjxN/75AX4Bq1arFzpw5s1CZc3JyKFeuXKH6hiuNOTwdOXKE559/nsWLFxMfH8/w4cM5//zzT9s+EsZ8LgYObMjx48eZMGGz11GCxt9jfM0112Q55+J+s8M597sPYDnwUQGPrsCBU9ru/53nSQIW5VuvAmzPt14T+OhMeZxzxMbGusJauXJlofuGK405fJ04ccJNmjTJlSxZ0tWpU8dt3LjxtG0jZcxnq1Ur5xo02O91jKDy9xgD610Bv1PP+NGQc66tc+7KAh7zgW/MrDqA7+e351Cc9gIXmFkJ33oNYM859BeJeGZG3759Wb16Nbm5uTRt2pTp06d7HUsijL9zBAuA3r7l3sD8s+3oq04rgZsL018kmiQmJpKVlUVcXBy33XYbAwcO5OjRo17HkgjhbyF4AmhnZtuAtr51zCzOzKacbGRm7wFzgDZmttvMOvh2DQb+bmbbgcrAVD/ziESsCy+8kNTUVO6//36ef/552rZtyzfffON1LIkAJc7c5PScc3uBNgVsXw/clW+95Wn67wTi/ckgEk1KlizJ+PHjadKkCXfffTeNGzfmrbfeomnTpl5HkzCmK4tFwlBycjJr166lVKlStGrVipdffllfTSGFpkIgEqYaNGjA+vXradu2Lf379+epp57i8OHDXseSMKRCIBLGKlWqxKJFixg5ciRLliyhRYsWZGdnex1LwowKgUiYK1asGI8++ihjxoxh+/btxMbGsnz5cq9jSRhRIRCJEM2aNSMzM5MLL7yQDh068MQTT2jeQM6KCoFIBKlXrx7p6enccsstDB06VF9pLWdFhUAkwpQrV44333yTZ599lgULFhAfH8/HH3/sdSwJYSoEIhHIzHjggQdITU1l//79JCQkMGfOHK9jSYhSIRCJYK1atSIrK4srr7yS7t278+CDD3Ls2DGvY0mIUSEQiXA1atRg1apVDBgwgKeffpp27drx7bfn8v2QEulUCESiQExMDBMnTuTVV18lPT2dxo0bk56e7nUsCREqBCJR5Pbbb2ft2rXExMRw9dVX89JLL+kUU1EhEIk2DRs2/OWrKQYMGMAdd9zBwYMHvY4lHlIhEIlC+b+a4rXXXqNZs2bs3LnT61jiERUCkSh18qspFi1aRHZ2NrGxsbz77rtexxIPqBCIRLlOnTqRlZVF7dq16dy5M6NGjeLEiRNex5IgUiEQES655BLef/99evfuzaOPPkrnzp3Zt2+f17EkSFQIRASAMmXK8K9//YuXXnqJ5cuXExsby4YNG7yOJUGgQiAivzAz+vXrx3vvvcexY8do1qwZr7zyitexpIipEIjIbyQkJLBhwwZatGjBnXfeyd133627n0UwFQIRKVDVqlVZunQpQ4cOZcqUKbr7WQRTIRCR0ypevDhjx45l/vz5bNu2jcaNG7N06VKvY0mAqRCIyBl16dKFrKwsLrroIq699loee+wxnWIaQVQIROSs1K1bl/T0dJKTkxk5ciRdunRh//79XseSAFAhEJGzVqZMGV577TVefPFFli1bplNMI4QKgYicEzOjf//+rF69mqNHj9KsWTOmTp3qdSzxgwqBiBRKYmLiL6eY3nXXXdx11106xTRM+VUIzKySmaWY2Tbfz4qnabfEzA6Y2aJTtk8zs8/NbKPv0dCfPCISXCdPMR02bBhTp06lefPmfP75517HknPk7zuCIUCqc64ekOpbL8g4oNdp9j3onGvoe2z0M4+IBFnx4sUZM2YMCxYsYMeOHcTGxrJ48WKvY8k58LcQdAVe9S2/CnQrqJFzLhX4yc/XEpEQdv3115OVlUWtWrXo1KkTI0eO5Pjx417HkrNg/tymzswOOOcu8C0bsP/kegFtk4BBzrnO+bZNA5oCufjeUTjnck/Tvy/QF6BatWqxM2fOLFTmnJwcypUrV6i+4Upjjg6hMubc3Fyef/55/v3vfxMXF8eIESOoUKFCwF9n4MCGHD9+nAkTNgf8uUOVv8f4mmuuyXLOxf1mh3Pudx/AcuCjAh5dgQOntN3/O8+TBCw6ZVt1wIBS5L2jGHmmPM45YmNjXWGtXLmy0H3DlcYcHUJtzFOmTHGlSpVyNWrUcGvXrg3487dq5VyDBvsD/ryhzN9jDKx3BfxOPeNHQ865ts65Kwt4zAe+MbPqAL6f355LdXLOfeXLlwv8C4g/l/4iErruvPNOPvjgA0qWLMnVV1/NCy+8cPIPQAkx/s4RLAB6+5Z7A/PPpXO+ImLkzS985GceEQkhjRs3Jisri/bt2/O3v/2N5ORkcnJyvI4lp/C3EDwBtDOzbUBb3zpmFmdmU042MrP3gDlAGzPbbWYdfLumm9lmYDNQBXjczzwiEmIqVqzIggULePzxx5k1axYJCQls3brV61iSTwl/Ojvn9gJtCti+Hrgr33rL0/Rv7c/ri0h4KFasGMOHDycxMZEePXrQpEkTpk6dSvfu3b2OJujKYhEJojZt2rBhwwauuuoq/vKXvzBw4ECOHDnidayop0IgIkFVo0YN0tLSuP/++3n++edJSkpi9+7dXseKaioEIhJ0MTExjB8/ntmzZ7N582YaNWpESkqK17GilgqBiHjmlltuITMzk2rVqtGhQwfd8MYjKgQi4qnLL7+cjIwMevbsyciRI+ncuTN79+71OlZUUSEQEc+VLVuW119/nRdffJHU1FQaN27MunXrvI4VNVQIRCQknLzhzZo1awBo0aIFEydO1NXIQaBCICIhpUmTJmzYsIF27dpx77336mrkIFAhEJGQU7lyZRYuXMiYMWOYNWsW8fHxfPLJJ17HilgqBCISkooVK8awYcNISUlh7969NGnShBkzZngdKyKpEIhISGvdujUbNmygUaNGJCcnM2DAAI6++iqkp3PBhxuhdm2YPt3rmGFNhUBEQt5FF13EihUrGDRoEAdeeoljffpA7uG8ndnZ0LevioEfVAhEJCyULFmScePGMaVqVc479aKzgwdh+HBvgkUAv759VEQk2Mp8/z0ADdn46x27dgU/TIRQIRCR8HLxxZCdzXge+O12KRR9NCQi4WXMGChT5tfbypTJ2y6FokIgIuElORkmT4ZatXBmUKtW3npystfJwpYKgYiEn+Rk+OILVq1YAV98oSLgJxUCEZEop0IgIhLlVAhERKKcCoGISJRTIRARiXIWjjd9MLPvgOxCdq8CfB/AOOFAY44OGnPk83e8tZxzVU/dGJaFwB9mtt45F+d1jmDSmKODxhz5imq8+mhIRCTKqRCIiES5aCwEk70O4AGNOTpozJGvSMYbdXMEIiLya9H4jkBERPJRIRARiXIRVQjMrKOZfWpm281sSAH7/25mH5vZJjNLNbNa+fb1NrNtvkfv4CYvHD/He9zMNvoeC4KbvPDOYsz9zGyzb1xrzOyKfPuG+vp9amYdgpu88Ao7ZjOrbWaH8h3nl4OfvnDONOZ87W4yM2dmcfm2ReRxztfuV2MOyHF2zkXEAygO7AAuAWKAD4ErTmlzDVDGt9wfmOVbrgTs9P2s6Fuu6PWYimq8vvUcr8dQRGM+P99yF2CJb/kKX/tSQB3f8xT3ekxFPObawEdej6EoxuxrVx5YDaQDcZF+nH9nzH4f50h6RxAPbHfO7XTOHQFmAl3zN3DOrXTOHfStpgM1fMsdgBTn3D7n3H4gBegYpNyF5c94w9XZjPnHfKtlgZNnQ3QFZjrncp1znwPbfc8X6vwZc7g645h9HgOeBA7n2xaxx9mnoDH7LZIKwUXAf/Ot7/ZtO507gX8Xsm8o8Ge8AKXNbL2ZpZtZtyLIVxTOasxmdo+Z7QCeAu47l74hyJ8xA9Qxs/+Y2Soza1m0UQPmjGM2s8ZATefcu+faN0T5M2bw8zhH5c3rzew2IA5o5XWWYDjNeGs55/aY2SXACjPb7Jzb4U3CwHLOTQQmmllPYAQQFnM+/jjNmL8CLnbO7TWzWGCemf3plHcQYcfMigHPAnd4HCVozjBmv49zJL0j2APUzLdew7ftV8ysLTAc6OKcyz2XviHGn/HinNvj+7kTSAMaFWXYADnX4zQT6FbIvqGi0GP2fTyy17ecRd5n0JcVTcyAOtOYywNXAmlm9gWQCCzwTZ5G6nE+7ZgDcpy9niQJ4GRLCfImeevwf5MtfzqlTSPff6R6p2yvBHxO3kRxRd9yJa/HVITjrQiU8i1XAbZRwMRUqD3Ocsz18i1fD6z3Lf+JX08i7iQ8JhH9GXPVk2MkbxJyT6j/uz7bMZ/SPo3/mziN2OP8O2P2+zhHzEdDzrljZnYvsJS8GfhXnHNbzGw0ef9jLADGAeWAOWYGsMs518U5t8/MHgMyfU832jm3z4NhnDV/xgvUByaZ2Qny3hU+4Zz72JOBnIOzHPO9vndBR4H9+D4W8rWbDXwMHAPucc4d92Qg58CfMQNXA6PN7ChwAugX6v+u4azHfLq+kXycT8fv46yvmBARiXKRNEcgIiKFoEIgIhLlVAhERKKcCoGISJRTIRARiXIqBCIiUU6FQEQkyv1/ctXoXp2E4NwAAAAASUVORK5CYII=", "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 \\\\\n", "g\\, := \\,( {x,y} )\\mapsto 10\\,x-2\\,y+1 \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": 20, "metadata": { "scrolled": true }, "outputs": [ { "data": { "application/javascript": "/* Put everything inside the global mpl namespace */\n/* global mpl */\nwindow.mpl = {};\n\nmpl.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(\n '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\nmpl.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 = document.createElement('div');\n this.root.setAttribute('style', 'display: inline-block');\n this._root_extra_style(this.root);\n\n parent_element.appendChild(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 (fig.ratio !== 1) {\n fig.send_message('set_device_pixel_ratio', {\n device_pixel_ratio: fig.ratio,\n });\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 fig.ws.close();\n };\n\n this.ws.onmessage = this._make_on_message_function(this);\n\n this.ondownload = ondownload;\n};\n\nmpl.figure.prototype._init_header = function () {\n var titlebar = document.createElement('div');\n titlebar.classList =\n 'ui-dialog-titlebar ui-widget-header ui-corner-all ui-helper-clearfix';\n var titletext = document.createElement('div');\n titletext.classList = 'ui-dialog-title';\n titletext.setAttribute(\n 'style',\n 'width: 100%; text-align: center; padding: 3px;'\n );\n titlebar.appendChild(titletext);\n this.root.appendChild(titlebar);\n this.header = titletext;\n};\n\nmpl.figure.prototype._canvas_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._root_extra_style = function (_canvas_div) {};\n\nmpl.figure.prototype._init_canvas = function () {\n var fig = this;\n\n var canvas_div = (this.canvas_div = document.createElement('div'));\n canvas_div.setAttribute(\n 'style',\n 'border: 1px solid #ddd;' +\n 'box-sizing: content-box;' +\n 'clear: both;' +\n 'min-height: 1px;' +\n 'min-width: 1px;' +\n 'outline: 0;' +\n 'overflow: hidden;' +\n 'position: relative;' +\n 'resize: both;'\n );\n\n function on_keyboard_event_closure(name) {\n return function (event) {\n return fig.key_event(event, name);\n };\n }\n\n canvas_div.addEventListener(\n 'keydown',\n on_keyboard_event_closure('key_press')\n );\n canvas_div.addEventListener(\n 'keyup',\n on_keyboard_event_closure('key_release')\n );\n\n this._canvas_extra_style(canvas_div);\n this.root.appendChild(canvas_div);\n\n var canvas = (this.canvas = document.createElement('canvas'));\n canvas.classList.add('mpl-canvas');\n canvas.setAttribute('style', 'box-sizing: content-box;');\n\n this.context = canvas.getContext('2d');\n\n var backingStore =\n this.context.backingStorePixelRatio ||\n this.context.webkitBackingStorePixelRatio ||\n this.context.mozBackingStorePixelRatio ||\n this.context.msBackingStorePixelRatio ||\n this.context.oBackingStorePixelRatio ||\n this.context.backingStorePixelRatio ||\n 1;\n\n this.ratio = (window.devicePixelRatio || 1) / backingStore;\n\n var rubberband_canvas = (this.rubberband_canvas = document.createElement(\n 'canvas'\n ));\n rubberband_canvas.setAttribute(\n 'style',\n 'box-sizing: content-box; position: absolute; left: 0; top: 0; z-index: 1;'\n );\n\n // Apply a ponyfill if ResizeObserver is not implemented by browser.\n if (this.ResizeObserver === undefined) {\n if (window.ResizeObserver !== undefined) {\n this.ResizeObserver = window.ResizeObserver;\n } else {\n var obs = _JSXTOOLS_RESIZE_OBSERVER({});\n this.ResizeObserver = obs.ResizeObserver;\n }\n }\n\n this.resizeObserverInstance = new this.ResizeObserver(function (entries) {\n var nentries = entries.length;\n for (var i = 0; i < nentries; i++) {\n var entry = entries[i];\n var width, height;\n if (entry.contentBoxSize) {\n if (entry.contentBoxSize instanceof Array) {\n // Chrome 84 implements new version of spec.\n width = entry.contentBoxSize[0].inlineSize;\n height = entry.contentBoxSize[0].blockSize;\n } else {\n // Firefox implements old version of spec.\n width = entry.contentBoxSize.inlineSize;\n height = entry.contentBoxSize.blockSize;\n }\n } else {\n // Chrome <84 implements even older version of spec.\n width = entry.contentRect.width;\n height = entry.contentRect.height;\n }\n\n // Keep the size of the canvas and rubber band canvas in sync with\n // the canvas container.\n if (entry.devicePixelContentBoxSize) {\n // Chrome 84 implements new version of spec.\n canvas.setAttribute(\n 'width',\n entry.devicePixelContentBoxSize[0].inlineSize\n );\n canvas.setAttribute(\n 'height',\n entry.devicePixelContentBoxSize[0].blockSize\n );\n } else {\n canvas.setAttribute('width', width * fig.ratio);\n canvas.setAttribute('height', height * fig.ratio);\n }\n canvas.setAttribute(\n 'style',\n 'width: ' + width + 'px; height: ' + height + 'px;'\n );\n\n rubberband_canvas.setAttribute('width', width);\n rubberband_canvas.setAttribute('height', height);\n\n // And update the size in Python. We ignore the initial 0/0 size\n // that occurs as the element is placed into the DOM, which should\n // otherwise not happen due to the minimum size styling.\n if (fig.ws.readyState == 1 && width != 0 && height != 0) {\n fig.request_resize(width, height);\n }\n }\n });\n this.resizeObserverInstance.observe(canvas_div);\n\n function on_mouse_event_closure(name) {\n return function (event) {\n return fig.mouse_event(event, name);\n };\n }\n\n rubberband_canvas.addEventListener(\n 'mousedown',\n on_mouse_event_closure('button_press')\n );\n rubberband_canvas.addEventListener(\n 'mouseup',\n on_mouse_event_closure('button_release')\n );\n rubberband_canvas.addEventListener(\n 'dblclick',\n on_mouse_event_closure('dblclick')\n );\n // Throttle sequential mouse events to 1 every 20ms.\n rubberband_canvas.addEventListener(\n 'mousemove',\n on_mouse_event_closure('motion_notify')\n );\n\n rubberband_canvas.addEventListener(\n 'mouseenter',\n on_mouse_event_closure('figure_enter')\n );\n rubberband_canvas.addEventListener(\n 'mouseleave',\n on_mouse_event_closure('figure_leave')\n );\n\n canvas_div.addEventListener('wheel', function (event) {\n if (event.deltaY < 0) {\n event.step = 1;\n } else {\n event.step = -1;\n }\n on_mouse_event_closure('scroll')(event);\n });\n\n canvas_div.appendChild(canvas);\n canvas_div.appendChild(rubberband_canvas);\n\n this.rubberband_context = rubberband_canvas.getContext('2d');\n this.rubberband_context.strokeStyle = '#000000';\n\n this._resize_canvas = function (width, height, forward) {\n if (forward) {\n canvas_div.style.width = width + 'px';\n canvas_div.style.height = height + 'px';\n }\n };\n\n // Disable right mouse context menu.\n this.rubberband_canvas.addEventListener('contextmenu', function (_e) {\n event.preventDefault();\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\nmpl.figure.prototype._init_toolbar = function () {\n var fig = this;\n\n var toolbar = document.createElement('div');\n toolbar.classList = 'mpl-toolbar';\n this.root.appendChild(toolbar);\n\n function on_click_closure(name) {\n return function (_event) {\n return fig.toolbar_button_onclick(name);\n };\n }\n\n function on_mouseover_closure(tooltip) {\n return function (event) {\n if (!event.currentTarget.disabled) {\n return fig.toolbar_button_onmouseover(tooltip);\n }\n };\n }\n\n fig.buttons = {};\n var buttonGroup = document.createElement('div');\n buttonGroup.classList = 'mpl-button-group';\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 /* Instead of a spacer, we start a new button group. */\n if (buttonGroup.hasChildNodes()) {\n toolbar.appendChild(buttonGroup);\n }\n buttonGroup = document.createElement('div');\n buttonGroup.classList = 'mpl-button-group';\n continue;\n }\n\n var button = (fig.buttons[name] = document.createElement('button'));\n button.classList = 'mpl-widget';\n button.setAttribute('role', 'button');\n button.setAttribute('aria-disabled', 'false');\n button.addEventListener('click', on_click_closure(method_name));\n button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n\n var icon_img = document.createElement('img');\n icon_img.src = '_images/' + image + '.png';\n icon_img.srcset = '_images/' + image + '_large.png 2x';\n icon_img.alt = tooltip;\n button.appendChild(icon_img);\n\n buttonGroup.appendChild(button);\n }\n\n if (buttonGroup.hasChildNodes()) {\n toolbar.appendChild(buttonGroup);\n }\n\n var fmt_picker = document.createElement('select');\n fmt_picker.classList = 'mpl-widget';\n toolbar.appendChild(fmt_picker);\n this.format_dropdown = fmt_picker;\n\n for (var ind in mpl.extensions) {\n var fmt = mpl.extensions[ind];\n var option = document.createElement('option');\n option.selected = fmt === mpl.default_extension;\n option.innerHTML = fmt;\n fmt_picker.appendChild(option);\n }\n\n var status_bar = document.createElement('span');\n status_bar.classList = 'mpl-message';\n toolbar.appendChild(status_bar);\n this.message = status_bar;\n};\n\nmpl.figure.prototype.request_resize = function (x_pixels, y_pixels) {\n // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n // which will in turn request a refresh of the image.\n this.send_message('resize', { width: x_pixels, height: y_pixels });\n};\n\nmpl.figure.prototype.send_message = function (type, properties) {\n properties['type'] = type;\n properties['figure_id'] = this.id;\n this.ws.send(JSON.stringify(properties));\n};\n\nmpl.figure.prototype.send_draw_message = function () {\n if (!this.waiting) {\n this.waiting = true;\n this.ws.send(JSON.stringify({ type: 'draw', figure_id: this.id }));\n }\n};\n\nmpl.figure.prototype.handle_save = function (fig, _msg) {\n var format_dropdown = fig.format_dropdown;\n var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n fig.ondownload(fig, format);\n};\n\nmpl.figure.prototype.handle_resize = function (fig, msg) {\n var size = msg['size'];\n if (size[0] !== fig.canvas.width || size[1] !== fig.canvas.height) {\n fig._resize_canvas(size[0], size[1], msg['forward']);\n fig.send_message('refresh', {});\n }\n};\n\nmpl.figure.prototype.handle_rubberband = function (fig, msg) {\n var x0 = msg['x0'] / fig.ratio;\n var y0 = (fig.canvas.height - msg['y0']) / fig.ratio;\n var x1 = msg['x1'] / fig.ratio;\n var y1 = (fig.canvas.height - msg['y1']) / fig.ratio;\n x0 = Math.floor(x0) + 0.5;\n y0 = Math.floor(y0) + 0.5;\n x1 = Math.floor(x1) + 0.5;\n y1 = Math.floor(y1) + 0.5;\n var min_x = Math.min(x0, x1);\n var min_y = Math.min(y0, y1);\n var width = Math.abs(x1 - x0);\n var height = Math.abs(y1 - y0);\n\n fig.rubberband_context.clearRect(\n 0,\n 0,\n fig.canvas.width / fig.ratio,\n fig.canvas.height / fig.ratio\n );\n\n fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n};\n\nmpl.figure.prototype.handle_figure_label = function (fig, msg) {\n // Updates the figure title.\n fig.header.textContent = msg['label'];\n};\n\nmpl.figure.prototype.handle_cursor = function (fig, msg) {\n fig.rubberband_canvas.style.cursor = msg['cursor'];\n};\n\nmpl.figure.prototype.handle_message = function (fig, msg) {\n fig.message.textContent = msg['message'];\n};\n\nmpl.figure.prototype.handle_draw = function (fig, _msg) {\n // Request the server to send over a new figure.\n fig.send_draw_message();\n};\n\nmpl.figure.prototype.handle_image_mode = function (fig, msg) {\n fig.image_mode = msg['mode'];\n};\n\nmpl.figure.prototype.handle_history_buttons = function (fig, msg) {\n for (var key in msg) {\n if (!(key in fig.buttons)) {\n continue;\n }\n fig.buttons[key].disabled = !msg[key];\n fig.buttons[key].setAttribute('aria-disabled', !msg[key]);\n }\n};\n\nmpl.figure.prototype.handle_navigate_mode = function (fig, msg) {\n if (msg['mode'] === 'PAN') {\n fig.buttons['Pan'].classList.add('active');\n fig.buttons['Zoom'].classList.remove('active');\n } else if (msg['mode'] === 'ZOOM') {\n fig.buttons['Pan'].classList.remove('active');\n fig.buttons['Zoom'].classList.add('active');\n } else {\n fig.buttons['Pan'].classList.remove('active');\n fig.buttons['Zoom'].classList.remove('active');\n }\n};\n\nmpl.figure.prototype.updated_canvas_event = function () {\n // Called whenever the canvas gets updated.\n this.send_message('ack', {});\n};\n\n// A function to construct a web socket function for onmessage handling.\n// Called in the figure constructor.\nmpl.figure.prototype._make_on_message_function = function (fig) {\n return function socket_on_message(evt) {\n if (evt.data instanceof Blob) {\n var img = evt.data;\n if (img.type !== 'image/png') {\n /* FIXME: We get \"Resource interpreted as Image but\n * transferred with MIME type text/plain:\" errors on\n * Chrome. But how to set the MIME type? It doesn't seem\n * to be part of the websocket stream */\n img.type = 'image/png';\n }\n\n /* Free the memory for the previous frames */\n if (fig.imageObj.src) {\n (window.URL || window.webkitURL).revokeObjectURL(\n fig.imageObj.src\n );\n }\n\n fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n img\n );\n fig.updated_canvas_event();\n fig.waiting = false;\n return;\n } else if (\n typeof evt.data === 'string' &&\n evt.data.slice(0, 21) === 'data:image/png;base64'\n ) {\n fig.imageObj.src = evt.data;\n fig.updated_canvas_event();\n fig.waiting = false;\n return;\n }\n\n var msg = JSON.parse(evt.data);\n var msg_type = msg['type'];\n\n // Call the \"handle_{type}\" callback, which takes\n // the figure and JSON message as its only arguments.\n try {\n var callback = fig['handle_' + msg_type];\n } catch (e) {\n console.log(\n \"No handler for the '\" + msg_type + \"' message type: \",\n msg\n );\n return;\n }\n\n if (callback) {\n try {\n // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n callback(fig, msg);\n } catch (e) {\n console.log(\n \"Exception inside the 'handler_\" + msg_type + \"' callback:\",\n e,\n e.stack,\n msg\n );\n }\n }\n };\n};\n\n// from https://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\nmpl.findpos = function (e) {\n //this section is from http://www.quirksmode.org/js/events_properties.html\n var targ;\n if (!e) {\n e = window.event;\n }\n if (e.target) {\n targ = e.target;\n } else if (e.srcElement) {\n targ = e.srcElement;\n }\n if (targ.nodeType === 3) {\n // defeat Safari bug\n targ = targ.parentNode;\n }\n\n // pageX,Y are the mouse positions relative to the document\n var boundingRect = targ.getBoundingClientRect();\n var x = e.pageX - (boundingRect.left + document.body.scrollLeft);\n var y = e.pageY - (boundingRect.top + document.body.scrollTop);\n\n return { x: x, y: y };\n};\n\n/*\n * return a copy of an object with only non-object keys\n * we need this to avoid circular references\n * https://stackoverflow.com/a/24161582/3208463\n */\nfunction simpleKeys(original) {\n return Object.keys(original).reduce(function (obj, key) {\n if (typeof original[key] !== 'object') {\n obj[key] = original[key];\n }\n return obj;\n }, {});\n}\n\nmpl.figure.prototype.mouse_event = function (event, name) {\n var canvas_pos = mpl.findpos(event);\n\n if (name === 'button_press') {\n this.canvas.focus();\n this.canvas_div.focus();\n }\n\n var x = canvas_pos.x * this.ratio;\n var y = canvas_pos.y * this.ratio;\n\n this.send_message(name, {\n x: x,\n y: y,\n button: event.button,\n step: event.step,\n guiEvent: simpleKeys(event),\n });\n\n /* This prevents the web browser from automatically changing to\n * the text insertion cursor when the button is pressed. We want\n * to control all of the cursor setting manually through the\n * 'cursor' event from matplotlib */\n event.preventDefault();\n return false;\n};\n\nmpl.figure.prototype._key_event_extra = function (_event, _name) {\n // Handle any extra behaviour associated with a key event\n};\n\nmpl.figure.prototype.key_event = function (event, name) {\n // Prevent repeat events\n if (name === 'key_press') {\n if (event.key === this._key) {\n return;\n } else {\n this._key = event.key;\n }\n }\n if (name === 'key_release') {\n this._key = null;\n }\n\n var value = '';\n if (event.ctrlKey && event.key !== 'Control') {\n value += 'ctrl+';\n }\n else if (event.altKey && event.key !== 'Alt') {\n value += 'alt+';\n }\n else if (event.shiftKey && event.key !== 'Shift') {\n value += 'shift+';\n }\n\n value += 'k' + event.key;\n\n this._key_event_extra(event, name);\n\n this.send_message(name, { key: value, guiEvent: simpleKeys(event) });\n return false;\n};\n\nmpl.figure.prototype.toolbar_button_onclick = function (name) {\n if (name === 'download') {\n this.handle_save(this, null);\n } else {\n this.send_message('toolbar_button', { name: name });\n }\n};\n\nmpl.figure.prototype.toolbar_button_onmouseover = function (tooltip) {\n this.message.textContent = tooltip;\n};\n\n///////////////// REMAINING CONTENT GENERATED BY embed_js.py /////////////////\n// prettier-ignore\nvar _JSXTOOLS_RESIZE_OBSERVER=function(A){var t,i=new WeakMap,n=new WeakMap,a=new WeakMap,r=new WeakMap,o=new Set;function s(e){if(!(this instanceof s))throw new TypeError(\"Constructor requires 'new' operator\");i.set(this,e)}function h(){throw new TypeError(\"Function is not a constructor\")}function c(e,t,i,n){e=0 in arguments?Number(arguments[0]):0,t=1 in arguments?Number(arguments[1]):0,i=2 in arguments?Number(arguments[2]):0,n=3 in arguments?Number(arguments[3]):0,this.right=(this.x=this.left=e)+(this.width=i),this.bottom=(this.y=this.top=t)+(this.height=n),Object.freeze(this)}function d(){t=requestAnimationFrame(d);var s=new WeakMap,p=new Set;o.forEach((function(t){r.get(t).forEach((function(i){var r=t instanceof window.SVGElement,o=a.get(t),d=r?0:parseFloat(o.paddingTop),f=r?0:parseFloat(o.paddingRight),l=r?0:parseFloat(o.paddingBottom),u=r?0:parseFloat(o.paddingLeft),g=r?0:parseFloat(o.borderTopWidth),m=r?0:parseFloat(o.borderRightWidth),w=r?0:parseFloat(o.borderBottomWidth),b=u+f,F=d+l,v=(r?0:parseFloat(o.borderLeftWidth))+m,W=g+w,y=r?0:t.offsetHeight-W-t.clientHeight,E=r?0:t.offsetWidth-v-t.clientWidth,R=b+v,z=F+W,M=r?t.width:parseFloat(o.width)-R-E,O=r?t.height:parseFloat(o.height)-z-y;if(n.has(t)){var k=n.get(t);if(k[0]===M&&k[1]===O)return}n.set(t,[M,O]);var S=Object.create(h.prototype);S.target=t,S.contentRect=new c(u,d,M,O),s.has(i)||(s.set(i,[]),p.add(i)),s.get(i).push(S)}))})),p.forEach((function(e){i.get(e).call(e,s.get(e),e)}))}return s.prototype.observe=function(i){if(i instanceof window.Element){r.has(i)||(r.set(i,new Set),o.add(i),a.set(i,window.getComputedStyle(i)));var n=r.get(i);n.has(this)||n.add(this),cancelAnimationFrame(t),t=requestAnimationFrame(d)}},s.prototype.unobserve=function(i){if(i instanceof window.Element&&r.has(i)){var n=r.get(i);n.has(this)&&(n.delete(this),n.size||(r.delete(i),o.delete(i))),n.size||r.delete(i),o.size||cancelAnimationFrame(t)}},A.DOMRectReadOnly=c,A.ResizeObserver=s,A.ResizeObserverEntry=h,A}; // eslint-disable-line\nmpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Left button pans, Right button zooms\\nx/y fixes axis, CTRL fixes aspect\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\\nx/y fixes axis\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n\nmpl.extensions = [\"eps\", \"jpeg\", \"pgf\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n\nmpl.default_extension = \"png\";/* global mpl */\n\nvar comm_websocket_adapter = function (comm) {\n // Create a \"websocket\"-like object which calls the given IPython comm\n // object with the appropriate methods. Currently this is a non binary\n // socket, so there is still some room for performance tuning.\n var ws = {};\n\n ws.binaryType = comm.kernel.ws.binaryType;\n ws.readyState = comm.kernel.ws.readyState;\n function updateReadyState(_event) {\n if (comm.kernel.ws) {\n ws.readyState = comm.kernel.ws.readyState;\n } else {\n ws.readyState = 3; // Closed state.\n }\n }\n comm.kernel.ws.addEventListener('open', updateReadyState);\n comm.kernel.ws.addEventListener('close', updateReadyState);\n comm.kernel.ws.addEventListener('error', updateReadyState);\n\n ws.close = function () {\n comm.close();\n };\n ws.send = function (m) {\n //console.log('sending', m);\n comm.send(m);\n };\n // Register the callback with on_msg.\n comm.on_msg(function (msg) {\n //console.log('receiving', msg['content']['data'], msg);\n var data = msg['content']['data'];\n if (data['blob'] !== undefined) {\n data = {\n data: new Blob(msg['buffers'], { type: data['blob'] }),\n };\n }\n // Pass the mpl event to the overridden (by mpl) onmessage function.\n ws.onmessage(data);\n });\n return ws;\n};\n\nmpl.mpl_figure_comm = function (comm, msg) {\n // This is the function which gets called when the mpl process\n // starts-up an IPython Comm through the \"matplotlib\" channel.\n\n var id = msg.content.data.id;\n // Get hold of the div created by the display call when the Comm\n // socket was opened in Python.\n var element = document.getElementById(id);\n var ws_proxy = comm_websocket_adapter(comm);\n\n function ondownload(figure, _format) {\n window.open(figure.canvas.toDataURL());\n }\n\n var fig = new mpl.figure(id, ws_proxy, ondownload, element);\n\n // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n // web socket which is closed, not our websocket->open comm proxy.\n ws_proxy.onopen();\n\n fig.parent_element = element;\n fig.cell_info = mpl.find_output_cell(\"
\");\n if (!fig.cell_info) {\n console.error('Failed to find cell for figure', id, fig);\n return;\n }\n fig.cell_info[0].output_area.element.on(\n 'cleared',\n { fig: fig },\n fig._remove_fig_handler\n );\n};\n\nmpl.figure.prototype.handle_close = function (fig, msg) {\n var width = fig.canvas.width / fig.ratio;\n fig.cell_info[0].output_area.element.off(\n 'cleared',\n fig._remove_fig_handler\n );\n fig.resizeObserverInstance.unobserve(fig.canvas_div);\n\n // Update the output cell to use the data from the current canvas.\n fig.push_to_output();\n var dataURL = fig.canvas.toDataURL();\n // Re-enable the keyboard manager in IPython - without this line, in FF,\n // the notebook keyboard shortcuts fail.\n IPython.keyboard_manager.enable();\n fig.parent_element.innerHTML =\n '';\n fig.close_ws(fig, msg);\n};\n\nmpl.figure.prototype.close_ws = function (fig, msg) {\n fig.send_message('closing', msg);\n // fig.ws.close()\n};\n\nmpl.figure.prototype.push_to_output = function (_remove_interactive) {\n // Turn the data on the canvas into data in the output cell.\n var width = this.canvas.width / this.ratio;\n var dataURL = this.canvas.toDataURL();\n this.cell_info[1]['text/html'] =\n '';\n};\n\nmpl.figure.prototype.updated_canvas_event = function () {\n // Tell IPython that the notebook contents must change.\n IPython.notebook.set_dirty(true);\n this.send_message('ack', {});\n var fig = this;\n // Wait a second, then push the new image to the DOM so\n // that it is saved nicely (might be nice to debounce this).\n setTimeout(function () {\n fig.push_to_output();\n }, 1000);\n};\n\nmpl.figure.prototype._init_toolbar = function () {\n var fig = this;\n\n var toolbar = document.createElement('div');\n toolbar.classList = 'btn-toolbar';\n this.root.appendChild(toolbar);\n\n function on_click_closure(name) {\n return function (_event) {\n return fig.toolbar_button_onclick(name);\n };\n }\n\n function on_mouseover_closure(tooltip) {\n return function (event) {\n if (!event.currentTarget.disabled) {\n return fig.toolbar_button_onmouseover(tooltip);\n }\n };\n }\n\n fig.buttons = {};\n var buttonGroup = document.createElement('div');\n buttonGroup.classList = 'btn-group';\n var button;\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 /* Instead of a spacer, we start a new button group. */\n if (buttonGroup.hasChildNodes()) {\n toolbar.appendChild(buttonGroup);\n }\n buttonGroup = document.createElement('div');\n buttonGroup.classList = 'btn-group';\n continue;\n }\n\n button = fig.buttons[name] = document.createElement('button');\n button.classList = 'btn btn-default';\n button.href = '#';\n button.title = name;\n button.innerHTML = '';\n button.addEventListener('click', on_click_closure(method_name));\n button.addEventListener('mouseover', on_mouseover_closure(tooltip));\n buttonGroup.appendChild(button);\n }\n\n if (buttonGroup.hasChildNodes()) {\n toolbar.appendChild(buttonGroup);\n }\n\n // Add the status bar.\n var status_bar = document.createElement('span');\n status_bar.classList = 'mpl-message pull-right';\n toolbar.appendChild(status_bar);\n this.message = status_bar;\n\n // Add the close button to the window.\n var buttongrp = document.createElement('div');\n buttongrp.classList = 'btn-group inline pull-right';\n button = document.createElement('button');\n button.classList = 'btn btn-mini btn-primary';\n button.href = '#';\n button.title = 'Stop Interaction';\n button.innerHTML = '';\n button.addEventListener('click', function (_evt) {\n fig.handle_close(fig, {});\n });\n button.addEventListener(\n 'mouseover',\n on_mouseover_closure('Stop Interaction')\n );\n buttongrp.appendChild(button);\n var titlebar = this.root.querySelector('.ui-dialog-titlebar');\n titlebar.insertBefore(buttongrp, titlebar.firstChild);\n};\n\nmpl.figure.prototype._remove_fig_handler = function (event) {\n var fig = event.data.fig;\n if (event.target !== this) {\n // Ignore bubbled events from children.\n return;\n }\n fig.close_ws(fig, {});\n};\n\nmpl.figure.prototype._root_extra_style = function (el) {\n el.style.boxSizing = 'content-box'; // override notebook setting of border-box.\n};\n\nmpl.figure.prototype._canvas_extra_style = function (el) {\n // this is important to make the div 'focusable\n el.setAttribute('tabindex', 0);\n // reach out to IPython and tell the keyboard manager to turn it's self\n // off when our div gets focus\n\n // location in version 3\n if (IPython.notebook.keyboard_manager) {\n IPython.notebook.keyboard_manager.register_events(el);\n } else {\n // location in version 2\n IPython.keyboard_manager.register_events(el);\n }\n};\n\nmpl.figure.prototype._key_event_extra = function (event, _name) {\n // Check for shift+enter\n if (event.shiftKey && event.which === 13) {\n this.canvas_div.blur();\n // select the cell after this one\n var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n IPython.notebook.select(index + 1);\n }\n};\n\nmpl.figure.prototype.handle_save = function (fig, _msg) {\n fig.ondownload(fig, null);\n};\n\nmpl.find_output_cell = function (html_output) {\n // Return the cell and output element which can be found *uniquely* in the notebook.\n // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n // IPython event is triggered only after the cells have been serialised, which for\n // our purposes (turning an active figure into a static one), is too late.\n var cells = IPython.notebook.get_cells();\n var ncells = cells.length;\n for (var i = 0; i < ncells; i++) {\n var cell = cells[i];\n if (cell.cell_type === 'code') {\n for (var j = 0; j < cell.output_area.outputs.length; j++) {\n var data = cell.output_area.outputs[j];\n if (data.data) {\n // IPython >= 3 moved mimebundle to data attribute of output\n data = data.data;\n }\n if (data['text/html'] === html_output) {\n return [cell, data, j];\n }\n }\n }\n }\n};\n\n// Register the function which deals with the matplotlib target/channel.\n// The kernel may be null if the page has been refreshed.\nif (IPython.notebook.kernel !== null) {\n IPython.notebook.kernel.comm_manager.register_target(\n 'matplotlib',\n mpl.mpl_figure_comm\n );\n}\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ ":20: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6. This is consistent with other Axes classes.\n", " plot3d = Axes3D(fig)\n" ] } ], "source": [ "%matplotlib notebook\n", "\n", "import ipympl\n", "from mpl_toolkits.mplot3d import Axes3D\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "def f(x,y):\n", " return 4*x+2*y-6*x*y\n", "def g(x,y):\n", " return 10*x-2*y+1\n", "\n", "x = np.arange(-3, 3, 0.25)\n", "y = np.arange(-3, 3, 0.25)\n", "X, Y = np.meshgrid(x, y)\n", "Z1 = f(X,Y)\n", "Z2 = g(X,Y)\n", "\n", "fig = plt.figure()\n", "plot3d = Axes3D(fig)\n", "plot3d.plot_surface(X,Y,Z1) \n", "plot3d.plot_surface(X,Y,Z2, color='r') \n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2020年度課題(中筋さんありがとう)\n", "\n", "1. 次に示した「例題:二分法とNewton法の収束性」および「解答例」をコピペして，pythonが動作することを確認せよ．\n", "\n", "1. 対象の関数を $f(x) = \\exp(-x)-2\\exp(-2x)$ として解答せよ．\n", "提出は2.だけでよい．\n", "\n", "ただし，func, dfuncは以下を使え．下の「exp関数に関する注意」参照" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def func(x):\n", " return np.exp(-x)-2*np.exp(-2*x)\n", "\n", "def df(x):\n", " return -np.exp(-x) + 4*np.exp(-2*x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 例題:二分法とNewton法の収束性\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "代数方程式に関する次の課題に答えよ．(2004年度期末試験)\n", "1. $\\exp(-x) = x^2$を二分法およびニュートン法で解け.\n", "1. $n$回目の値$x_n$と小数点以下10桁まで求めた値$x_f=0.7034674225$との差$\\Delta x_n$の絶対値(abs)のlogを$n$の関数としてプロットし，その収束性を比較せよ．また，その傾きの違いを両解法の原理から説明せよ.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 解答例 \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "+ funcで関数を定義．\n", "+ 関数をplotして概形を確認．\n", "+ 組み込みコマンドで正解を確認しておく．\n", "$$\n", "0.7034674224983916520498186018599021303429\n", "$$\n", "テキストからプログラムをコピーして走らせてみる．環境によっては，printf分の中の\"\\\"が文字化けしているので，その場合は修正して使用せよ．\n", "\n", "+ プロットのためにリストをlist_bisecで作成している．\n", "+ 同様にNewton法での結果をlist_newtonに入れる．\n", "+ list_bisec, list_newtonを片対数プロットして同時に表示．\n", "\n", "2分法で求めた解は，Newton法で求めた解よりもゆっくりと精密解へ収束している．これは，二分法が原理的に計算回数について一次収束なのに対して，Newton法は2次収束であるためである．解の差($\\delta$)だけでなく，関数値$f(x),\\epsilon$をとっても同様の振る舞いを示す．\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-2*x - exp(-x)\n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from sympy import *\n", "\n", "x = symbols('x')\n", "\n", "\n", "def func(x):\n", " return exp(-x)-x**2\n", "\n", "def df(x):\n", " return diff(func(x), x)\n", "\n", "print(df(x))\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApqElEQVR4nO3deXgUZdb38e9hMwxhUzAwLIIISDQSCYaAgOwgYgDZZfUBGQdQUYIwOggibigyOIDjDCBEwMgOggMjEORhSwg7CSAYFVEcHRcg+MiW8/6RZt4Mk0BId7rSXedzXX3R1XV33+cY/FGpqq4SVcUYY0zwK+J0AcYYY/zDAt8YY1zCAt8YY1zCAt8YY1zCAt8YY1yimNMF5KZChQpao0aNfL//7NmzlCpVyncFBQC39ey2fsF6dgtvet61a9e/VLViTusKbeDXqFGDlJSUfL9/06ZNtGjRwncFBQC39ey2fsF6dgtvehaRL3NbZ7t0jDHGJSzwjTHGJSzwjTHGJSzwjTHGJSzwjTHGJXwS+CIyR0S+E5GDuawXEXlLRI6JyH4RaeCLeY0xxuSdr7bw5wIdrrL+fqC25zEUeNtH8xpjjMkjnwS+qm4GfrzKkM5AvGbZAZQTkcq+mDuHWhg9ejQHDhzALv1sjDH/n/gqFEWkBrBaVe/MYd1q4FVV3eJZ3gCMUdWUK8YNJes3AMLCwqISEhKuu46vv/6a3/3ud5w9e5ZbbrmFTp060a5dO8qUKXP9TQWYjIwMQkNDnS7Db9zWL1jPbuFNzy1bttylqg1zXKmqPnkANYCDuaxbDTTNtrwBaHi1z4uKitL8ysjI0NGjR2t0dLQCGhISogMGDNBt27ZpZmZmvj+3sEtMTHS6BL9yW7+q1rNbeNMzkKK55Kq/ztL5GqiWbbmq57UCUapUKTp27EhSUhJ79uzhkUceYfny5TRp0oTIyEjefvttzpw5U1DTG2NMoeSvwF8FDPCcrRMDnFLVk/6YODIykpkzZ/LNN9/wzjvvULRoUYYNG8Zvf/tbHnvsMfbt2+ePMowxxnG+Oi3zfWA7UFdETojIYBF5TEQe8wz5CEgHjgF/A4b5Yt7rERoaytChQ9m1axdJSUl0796defPmERkZyb333suCBQs4d+6cv8syxhi/8dVZOn1UtbKqFlfVqqo6W1X/oqp/8axXVR2uqrVUNUKvOFjrTyJCdHQ07777Ll9//TVTpkzh+++/p1+/flSrVo1nn32W48ePO1WeMcYUGFd/0/bGG2/k6aef5vDhw6xbt47GjRvz2muvUbNmTbp27cqGDRvs1E5jTNBwdeBfVqRIEdq1a8fKlStJT09nzJgxbNmyhTZt2hAeHs6MGTPsIK8xJuBZ4F/hlltu4eWXX+arr74iPj6e0qVLM2LECKpUqcITTzzBp59+6nSJxhiTLxb4uQgJCaF///4kJyezY8cOYmNj+ctf/kLdunXp2LEja9euJTMz0+kyjTEmzyzw86BRo0bMnz+f48ePM2HCBHbv3s39999PeHg4M2fOJCMjw+kSjTHmmizwr0OlSpUYP348x48fZ/78+ZQuXZrhw4dTtWpVRo8ebWf3GGMKNQv8fChRogR9+/YlOTmZbdu20b59e6ZOncqtt95Kjx492LZtm53dY4wpdCzwvSAiNG7cmA8++IDPP/+cuLg41q9fz7333ktMTAwffPABFy9edLpMY4wBLPB9plq1arz66qucOHGCGTNm8NNPP9G7d29q1arFlClTOHXqlNMlGmNczgLfx0qVKsWwYcM4fPgwK1eu5NZbbyUuLo5q1aoxatQo289vjHGMBX4BKVKkCLGxsSQmJpKSksKDDz7ItGnTuPXWW3n44YfZvXu30yUaY1zGAt8PoqKiWLBgAenp6YwcOZLVq1cTFRVF69at+fvf/24HeI0xfmGB70fVq1fnjTfe4KuvvuL111/nyJEjdOzYkfr16xMfH8/58+edLtEYE8Qs8B1QtmxZ4uLiSE9PZ968eWRmZjJw4EBq1arFm2++adftMcYUCAt8B5UoUYIBAwZw4MAB1qxZQ61atRg1ahTVq1fnj3/8I999953TJRpjgoivboDSQUSOiMgxERmbw/rqIpIoIntEZL+IdPTFvMFCROjYsSObNm1ix44dtGzZkpdffplbbrmFESNG8MUXXzhdojEmCHgd+CJSFJgB3A+EA31EJPyKYX8EFqnq3UBvYKa38warRo0asWzZMg4dOkTfvn3561//ym233Ua/fv04ePCg0+UZYwKYL7bwo4FjqpququeBBKDzFWMUKON5Xhb4xgfzBrW6desya9YsPv/8c5588klWrFhBREQEXbp0ISkpyenyjDEBSLw9JVBEugMdVHWIZ7k/0EhVR2QbUxn4B1AeKAW0UdVdOXzWUGAoQFhYWFRCQkK+68rIyCA0NDTf7y9sTp06xYoVK1i6dClnzpyhQYMG9OvXj8jISEQECL6er8Vt/YL17Bbe9NyyZctdqtowx5Wq6tUD6A7MyrbcH5h+xZingVGe542BNKDI1T43KipKvZGYmOjV+wurM2fO6BtvvKGVKlVSQGNiYnT16tWamZkZtD3nxm39qlrPbuFNz0CK5pKrvtil8zVQLdtyVc9r2Q0GFnn+gdkOhAAVfDC364SGhjJq1Cg+//xzZs6cycmTJ+nUqRMNGjTgk08+sZuyGGNy5YvA3wnUFpGaIlKCrIOyq64YcxxoDSAi9cgK/O99MLdrhYSE8Pvf/56jR48yd+5cfvnlFyZMmEBERAQLFy7k0qVLTpdojClkvA58Vb0IjADWAYfIOhsnVUQmikisZ9go4FER2Qe8Dwzy/OphvFS8eHEGDhxIWloa48aNQ0To27cv9erVY968eXZ5ZmPMv/nkPHxV/UhV66hqLVV9yfPa86q6yvM8TVXvVdX6qhqpqv/wxbzm/ytatCitWrVi//79LF26lFKlSjFo0CDq1q3LnDlzuHDhgtMlGmMcZt+0DTJFihThoYceYvfu3axcuZLy5cszePBg6taty+zZsy34jXExC/wgJSLExsayc+dOVq9eTYUKFRgyZAh16tRh1qxZFvzGuJAFfpATER544AGSkpJYs2YNFStW5NFHH7VdPca4kAW+S1y+Xk9SUhKrV6/mxhtvZPDgwdSrV4/4+Hg7uGuMC1jgu8zlLf6dO3eycuVKSpcuzcCBA7njjjvsdE5jgpwFvktd3se/e/duli1bxg033EDfvn256667WLJkiX2By5ggZIHvciJC165d2bt3Lx988AGZmZn06NGDqKgo1qxZY7dfNCaIWOAbIOt0zp49e3Lw4EHi4+M5ffo0nTp1okmTJmzcuNHp8owxPmCBb/5D0aJF6d+/P4cPH+add97hxIkTtG7dmrZt25KcnOx0ecYYL1jgmxwVL16coUOHcvToUaZOncrevXtp1KgRXbp0ITU11enyjDH5YIFvriokJISRI0eSnp7OxIkTSUxMJCIiggEDBvD55587XZ4x5jpY4Js8KV26NOPGjSM9PZ24uDgWL15M3bp1eeKJJ+xm68YECAt8c11uuukmJk+ezLFjx3jkkUeYOXMmt956K+PHj+fMmTNOl2eMuQoLfJMvVapU4Z133iEtLY2OHTsyceJEatWqxVtvvcW5c+ecLs8YkwMLfOOVOnXqsGjRIpKTk7nzzjt58sknqVevHu+//759ecuYQsYC3/jEPffcw4YNG1i7di1lypTh4YcfJjo62s7hN6YQ8Ungi0gHETkiIsdEZGwuY3qKSJqIpIrIQl/MawoXEaF9+/bs3r2b+Ph4vv/+e1q3bk3Hjh05cOCA0+UZ43peB76IFAVmAPcD4UAfEQm/Ykxt4A/Avap6BzDS23lN4VWkSBH69+/PkSNHeP3119m+fTv169fnf/7nf/j66yvvb2+M8RdfbOFHA8dUNV1VzwMJQOcrxjwKzFDVnwBU1c7jc4GQkBDi4uL47LPPeOqpp1iwYAG1a9dm3LhxdkaPMQ4Qby+OJSLdgQ6qOsSz3B9opKojso1ZAXwK3AsUBSao6tocPmsoMBQgLCwsKiEhId91ZWRkEBoamu/3B6LC3vPJkyeZNWsWGzdupHz58gwaNIgHHniAokWL5uvzCnu/BcF6dgdvem7ZsuUuVW2Y40pV9eoBdAdmZVvuD0y/YsxqYDlQHKgJfAWUu9rnRkVFqTcSExO9en8gCpSek5KStGnTpgpoeHi4fvTRR5qZmXndnxMo/fqS9ewO3vQMpGguueqLXTpfA9WyLVf1vJbdCWCVql5Q1c/J2tqv7YO5TQCKjo5m8+bNLFu2jPPnz9OxY0c6dOjAwYMHnS7NmKDmi8DfCdQWkZoiUgLoDay6YswKoAWAiFQA6gDpPpjbBKjL1+FPTU3lT3/6Ezt37qR+/fo89thjdqkGYwqI14GvqheBEcA64BCwSFVTRWSiiMR6hq0DfhCRNCARGK2qP3g7twl8JUqU4Mknn+TYsWOMGDGC2bNnU7t2bSZPnmzf2DXGx3xyHr6qfqSqdVS1lqq+5HnteVVd5Xmuqvq0qoaraoSq5v9orAlKN954I9OmTePgwYM0b96cMWPGEB4ezrJly+yuW8b4iH3T1hQqdevW5cMPP2TdunWULFmSbt260apVK/bt2+d0acYEPAt8Uyi1a9eOvXv3MmPGDA4cOECDBg147LHH+Ne//uV0acYELAt8U2gVK1aMYcOGcfToUUaMGMGsWbOoXbs206ZN48KFC06XZ0zAscA3hV758uWZNm0a+/fvJzo6mpEjRxIZGcmuXbucLs2YgGKBbwJGeHg4a9euZeXKlfz666/ExcXRtWtXu9WiMXlkgW8CiogQGxtLamoqQ4YM4eOPP6ZevXo8//zz/PLLL06XZ0yhZoFvAlJISAh9+/bl8OHDPPTQQ7z44ovUq1ePpUuX2mmcxuTCAt8EtKpVq7Jw4UI++eQTypUrR/fu3Wnfvj1HjhxxujRjCh0LfBMUmjdvzq5du3jrrbdITk4mIiKCMWPGkJGR4XRpxhQaFvgmaBQrVozHH3+cTz/9lH79+jF58mTq1avH4sWLbTePMVjgmyB08803M2fOHLZt20bFihXp2bMn7dq1s908xvUs8E3Qaty4MTt37mT69Ons3LmTiIgInn32Wc6ePet0acY4wgLfBLWiRYsyfPhwjhw5Qp8+fXjllVe44447+PDDD50uzRi/s8A3rhAWFsa8efP45JNPCA0NJTY2ls6dO/Pll186XZoxfmOBb1ylefPm7Nmzh8mTJ7N+/XrCw8OZPHmyXZvHuIJPAl9EOojIERE5JiJjrzKum4ioiOR8g11j/KB48eKMHj2aQ4cO0a5dO8aMGUODBg3YunWr06UZU6C8DnwRKQrMAO4HwoE+IhKew7jSwJNAkrdzGuML1atXZ/ny5axcuZLTp0/TtGlTHn30UX788UenSzOmQPhiCz8aOKaq6ap6HkgAOucw7kXgNeBXH8xpjM/ExsaSlpbG6NGjeffdd7n99tuZP3++nbtvgo54+5daRLoDHVR1iGe5P9BIVUdkG9MAeE5Vu4nIJiBOVVNy+KyhwFCAsLCwqISE/N8JMSMjg9DQ0Hy/PxC5reeC6PfYsWO8+eabHDp0iKioKJ566imqVKni0zm84bafMVjP16tly5a7VDXn3eaq6tUD6A7MyrbcH5iebbkIsAmo4VneBDS81udGRUWpNxITE716fyByW88F1e/Fixd1xowZWqZMGQ0JCdGXXnpJz58/XyBzXS+3/YxVrefrBaRoLrnqi106XwPVsi1X9bx2WWngTmCTiHwBxACr7MCtKayKFi3KsGHDOHToEA888ADPPfccDRo0YMeOHU6XZoxXfBH4O4HaIlJTREoAvYFVl1eq6ilVraCqNVS1BrADiNUcdukYU5j89re/ZcmSJaxcuZKff/6ZJk2a8Pjjj3P69GmnSzMmX7wOfFW9CIwA1gGHgEWqmioiE0Uk1tvPN8Zplw/qjhgxghkzZnDHHXewevVqp8sy5rr55Dx8Vf1IVeuoai1Vfcnz2vOquiqHsS1s694EmtKlS/PWW2+xbds2ypYty4MPPkivXr345z//6XRpxuSZfdPWmOsQExPD7t27mThxIitWrKBevXrMmzfPTuE0AcEC35jrVKJECcaNG8fevXsJDw9n0KBBdOjQgS+++MLp0oy5Kgt8Y/KpXr16bN68menTp7Nt2zbuvPNO/vznP5OZmel0acbkyALfGC8UKVKE4cOHk5qaSrNmzXjiiSdo1qwZhw8fdro0Y/6LBb4xPlC9enU++ugj5s2bx6FDh4iMjOS1117j4sWLTpdmzL9Z4BvjIyLCgAEDSEtLo2PHjowdO5aYmBgOHDjgdGnGABb4xvhcpUqVWLp0KYsWLeL48eNERUUxceJEu+a+cZwFvjEFQETo0aMHaWlpdOvWjfHjx3PPPfewd+9ep0szLmaBb0wBqlChAu+//z7Lly/n22+/5Z577uGFF16wrX3jCAt8Y/ygS5cupKam0rNnTyZMmEB0dDT79u1zuizjMhb4xvjJTTfdxIIFC1i+fDknT57knnvuYdKkSba1b/zGAt8YP7u8td+tWzfGjRtH48aNSU1Ndbos4wIW+MY44KabbuL9999n8eLFfPnllzRo0IDXX3+dS5cuOV2aCWIW+MY4qHv37qSmptKxY0eeeeYZmjdvzrFjx5wuywQpC3xjHHbzzTezbNky3nvvPVJTU6lfvz4zZ860K3Aan/NJ4ItIBxE5IiLHRGRsDuufFpE0EdkvIhtE5BZfzGtMsBAR+vXr9+9r8gwfPpz27dtz4sQJp0szQcTrwBeRosAM4H4gHOgjIuFXDNtD1o3L7wKWAJO9ndeYYFSlShX+/ve/8/bbb7N161YiIiJYuHChbe0bn/DFFn40cExV01X1PJAAdM4+QFUTVfUXz+IOsm50bozJgYjw2GOPsX//fsLDw+nbty+9e/fmhx9+cLo0E+DE2y0HEekOdFDVIZ7l/kAjVR2Ry/jpwLeqOimHdUOBoQBhYWFRCQkJ+a4rIyOD0NDQfL8/ELmtZzf0e+nSJRISEpg7dy5ly5bliSeeoHnz5k6X5Vdu+DlfyZueW7ZsuUtVG+a4UlW9egDdgVnZlvsD03MZ24+sLfwbrvW5UVFR6o3ExESv3h+I3Nazm/rds2eP3nHHHQro8OHD9ezZs06X5Ddu+jlf5k3PQIrmkqu+2KXzNVAt23JVz2v/QUTaAM8Bsap6zgfzGuMakZGRpKSk0L17d2bMmEGDBg3YtWuX02WZAOOLwN8J1BaRmiJSAugNrMo+QETuBt4hK+y/88GcxrhOSEgIw4cPZ/369WRkZBATE8Mrr7xiX9YyeeZ14KvqRWAEsA44BCxS1VQRmSgisZ5hrwOhwGIR2Ssiq3L5OGPMNbRu3ZoDBw7w0EMP8eyzz9KiRQu7gbrJE5+ch6+qH6lqHVWtpaoveV57XlVXeZ63UdUwVY30PGKv/onGmKspX748CQkJvPfee+zfv5/69euzYMECp8syhZx909aYAHX5y1r79u3jrrvuol+/fjz88MP8/PPPTpdmCikLfGMCXI0aNdi0aROTJk1i0aJF1K9fn//93/91uixTCFngGxMEihYtynPPPce2bdsoXrw4LVq0YNy4cXatffMfLPCNCSLR0dHs2bOHAQMGMGnSJJo3b056errTZZlCwgLfmCBTunRp3n33XRISEjh06BCRkZEsXLjQ6bJMIWCBb0yQ6tWrF3v37uWuu+6ib9++DBw4kDNnzjhdlnGQBb4xQezyAd3x48czf/58+4auy1ngGxPkihUrxoQJE9i0aRPnzp2jcePGvPnmm2RmZjpdmvEzC3xjXKJZs2bs3buXTp06MWrUKDp16sR339mVTtzEAt8YF7nxxhtZunQpM2bMYOPGjdSvX5+NGzc6XZbxEwt8Y1xGRBg2bBjJycmUK1eONm3aMG7cOC5evOh0aaaAWeAb41J33XUXKSkpDBo0iEmTJtGqVSu7h26Qs8A3xsVKlSrFnDlzmD9/Pnv27CEyMpI1a9Y4XZYpIBb4xhj69u3Lrl27qFq1Kp06dSIuLo7z5887XZbxMQt8YwwAderUYceOHQwbNowpU6bQvHlzvvzyS6fLMj7kk8AXkQ4ickREjonI2BzW3yAiH3jWJ4lIDV/Ma4zxrZCQEGbMmMGiRYtIS0vj7rvv5sMPP3S6LOMjXge+iBQFZgD3A+FAHxEJv2LYYOAnVb0NmAq85u28xpiC06NHD3bv3k2NGjWIjY0lLi7OrrwZBCTrJudefIBIY2CCqrb3LP8BQFVfyTZmnWfMdhEpBnwLVNSrTN6wYUNNSUnJV00vfJjKtrTjlCtXLl/vD1Q///yzq3p2W7/g/54zMzP57LPP+OabbyhTpgzh4eHccMMNfpsf3PlzLpN5mr/9vn2+3isiu1S1YU7rinlVVZYqwFfZlk8AjXIbo6oXReQUcBPwrysKHQoMBQgLC2PTpk35KujEiXNcunTJdXf+cVvPbusXnOm5YsWKFC9enK+++oqdO3dSvXp1ypQp47f53fhzLlnyUr7z76pU1asH0B2YlW25PzD9ijEHgarZlj8DKlztc6OiotQbiYmJXr0/ELmtZ7f1q+psz59++qnWr19fAX322Wf1woULfpnXfs7XB0jRXHLVFwdtvwaqZVuu6nktxzGeXTplgR98MLcxxk9q167N9u3bGTx4MC+//DJt27bl22+/dboscx18Efg7gdoiUlNESgC9gVVXjFkFDPQ87w5s9PxLZIwJICVLlmTWrFm8++67JCUl0aBBAzZv3ux0WSaPvA58Vb0IjADWAYeARaqaKiITRSTWM2w2cJOIHAOeBv7r1E1jTOAYNGgQSUlJhIaG0qpVK15//XVsG67w88VBW1T1I+CjK157PtvzX4EevpjLGFM4REREkJKSwuDBg3nmmWfYtm0bc+fOpWzZsk6XZnJh37Q1xuRbmTJlWLRoEVOnTmX16tU0bNiQ/fv3O12WyYUFvjHGKyLCyJEjSUxM5OzZs8TExBAfH+90WSYHFvjGGJ9o2rQpe/bsoVGjRgwcOJBhw4Zx7tw5p8sy2VjgG2N8JiwsjI8//pjRo0fz9ttv07x5c7766qtrv9H4hQW+McanihUrxuTJk1myZAmHDh2iQYMGJCYmOl2WwQLfGFNAunXrRnJyMhUrVqRNmza88cYbduqmwyzwjTEF5vbbbycpKYmuXbsyevRoevXqRUZGhtNluZYFvjGmQJUuXZrFixfz2muvsXTpUho1asTRo0edLsuVLPCNMQVORHjmmWdYt24d//znP2nYsKHdWMUBFvjGGL9p06YNKSkp3HbbbcTGxvLCCy+QmZnpdFmuYYFvjPGrGjVqsGXLFgYMGMCECRPo2rUrp0+fdrosV7DAN8b4XcmSJZk7dy7Tpk1jzZo1REdHc+TIEafLCnoW+MYYR4gITzzxBBs2bODHH38kOjqa1atXO11WULPAN8Y46r777vuP/fqTJk2y/foFxALfGOO46tWrs2XLFvr27cu4cePo2bOnna9fALwKfBG5UUQ+FpGjnj/L5zAmUkS2i0iqiOwXkV7ezGmMCU4lS5YkPj6eKVOmsHz5cpo0aUJ6errTZQUVb7fwxwIbVLU2sIGc72T1CzBAVe8AOgB/EpFyXs5rjAlCIsLTTz/N2rVrOXHiBPfccw+7d+92uqyg4W3gdwbmeZ7PA7pcOUBVP1XVo57n3wDfARW9nNcYE8Tatm1LcnIylSpVYvTo0UyfPt2uw+MD3gZ+mKqe9Dz/Fgi72mARiQZKAJ95Oa8xJsjddtttbN++nZiYGB5//HGGDh3K+fPnnS4roMm1/tUUkfVApRxWPQfMU9Vy2cb+pKr/tR/fs64ysAkYqKo7chkzFBgKEBYWFpWQkJCHFnKWkZFBaGhovt8fiNzWs9v6BXf2fPr0aRYvXsz8+fOJiIjghRdeoHz5HGMmaHjzc27ZsuUuVW2Y40pVzfcDOAJU9jyvDBzJZVwZYDfQPa+fHRUVpd5ITEz06v2ByG09u61fVXf3nJCQoCEhIVq9enXdu3evs0UVMG9+zkCK5pKr3u7SWQUM9DwfCKy8coCIlACWA/GqusTL+YwxLtWrVy+2bNnCpUuXaNKkCcuXL3e6pIDjbeC/CrQVkaNAG88yItJQRGZ5xvQEmgODRGSv5xHp5bzGGBeKiopi586dRERE8NBDD/HSSy/ZwdzrUMybN6vqD0DrHF5PAYZ4ns8H5nszjzHGXFa5cmU2bdrEkCFD+OMf/8jBgweZM2cOJUuWdLq0Qs++aWuMCTghISG89957vPLKK3zwwQfcd999nDx58tpvdDkLfGNMQBIRxo4dy/Lly0lLSyM6Opo9e/Y4XVahZoFvjAlonTt3ZsuWLYgITZs2tYO5V2GBb4wJeJGRkSQnJ//7YO5rr71mB3NzYIFvjAkKlSpVIjExkT59+jB27FgeeeQRzp0753RZhYpXZ+kYY0xhUrJkSRYsWMDtt9/O+PHjSU9PZ9myZVSoUMHp0goF28I3xgQVEeH5558nISGB5ORkYmJi7PaJHhb4xpig1KtXLxITEzl9+jQxMTFs3LjR6ZIcZ4FvjAlajRs3Jjk5mSpVqtC+fXvmzJnjdEmOssA3xgS1GjVqsHXrVlq2bMngwYP5wx/+4Np75lrgG2OCXtmyZVmzZg2/+93vePXVV+nVqxf/93//53RZfmdn6RhjXKF48eK8/fbb1K5dm9GjR3PixAlWrlzJzTff7HRpfmNb+MYY1xARRo0axZIlS9i3bx8xMTEcPnzY6bL8xgLfGOM6Dz30EJs2beLs2bM0btyYTz75xOmS/MIC3xjjStHR0SQlJVG5cmXatm3LggULnC6pwFngG2Nc6/IZPPfeey/9+vVj0qRJQX0NHq8CX0RuFJGPReSo589c7ywsImVE5ISITPdmTmOM8aXy5cuzbt06+vfvz7hx43j00Ue5cOGC02UVCG+38McCG1S1NrDBs5ybF4HNXs5njDE+V6JECebNm8e4ceOYPXs2Dz74IGfOnHG6LJ/zNvA7A/M8z+cBXXIaJCJRQBjwDy/nM8aYAiEiTJw4kdmzZ7N+/XqaN2/ON99843RZPiXe7K8SkZ9VtZznuQA/XV7ONqYIsBHoR9aNzhuq6ohcPm8oMBQgLCwsKiEhId+1ZWRkEBoamu/3ByK39ey2fsF69pfk5GQmTJhA6dKlefXVV6lZs6Zf5/em55YtW+5S1YY5rlTVqz6A9cDBHB6dgZ+vGPtTDu8fATzjeT4ImH6tOVWVqKgo9UZiYqJX7w9EbuvZbf2qWs/+tHv3bq1cubKWLVvW7zV4Mx+Qornk6jW/aauqbXJbJyL/FJHKqnpSRCoD3+UwrDHQTESGAaFACRHJUNWr7e83xhhH3X333Wzfvp3777+fdu3aER8fT+/evZ0uyyve7sNfBQz0PB8IrLxygKr2VdXqqloDiAPiLeyNMYHglltuYevWrTRu3Jg+ffowZcqUgD5t09vAfxVoKyJHydo//yqAiDQUkVneFmeMMU67fNpmjx49iIuL46mnngrYq216dfE0Vf0BaJ3D6ynAkBxenwvM9WZOY4zxt5CQEBISEqhSpQp/+tOf+Oabb4iPjyckJMTp0q6LXS3TGGPyoEiRIkydOpWqVasSFxfHd999x4oVKyhXrpzTpeWZXVrBGGOuw6hRo1i4cCHbtm2jadOmnDhxwumS8swC3xhjrlOfPn1Yu3Ytx48fp0mTJqSlpTldUp5Y4BtjTD60atWKzZs3c+HCBZo2bcrWrVudLumaLPCNMSafIiMj2bZtGxUqVKBNmzasWrXK6ZKuygLfGGO8ULNmTbZu3UpERARdu3Zl9uzZTpeUKwt8Y4zxUsWKFdm4cSNt27ZlyJAhvPLKK4XyC1oW+MYY4wOhoaGsWrWKvn378uyzzzJy5MhC9wUtOw/fGGN8pESJEsTHx3PzzTczdepUvv/+e+bOnUuJEiWcLg2wwDfGGJ8qUqQIU6ZMISwsjLFjx/Ljjz+yZMmSQnFZa9ulY4wxPiYijBkzhlmzZvHxxx/Tpk0bfvjhB6fLssA3xpiCMnjwYJYuXcrevXtp1qyZ49/KtcA3xpgC1KVLF9auXcuJEye49957OXLkiGO1WOAbY0wBa9GiBZ988gm//vorzZo1Y/fu3Y7UYYFvjDF+cPfdd7NlyxZ+85vf0KJFCzZt2uT3GrwKfBG5UUQ+FpGjnj/L5zKuuoj8Q0QOiUiaiNTwZl5jjAlEtWvXZuvWrVSrVo0OHTqwcuV/3SSwQHm7hT8W2KCqtYENnuWcxAOvq2o9IJqc731rjDFBr0qVKmzevJn69evTrVs33nvvPb/N7W3gdwbmeZ7PA7pcOUBEwoFiqvoxgKpmqOovXs5rjDEB66abbmL9+vXcd999DBgwgD//+c9+mVe8ud6DiPysquU8zwX46fJytjFdyLrd4XmgJrAeGKuql3L4vKHAUICwsLCohISEfNeWkZFRKL7o4E9u69lt/YL1HGzOnz/Piy++yJYtWxg0aBADBgxARLzquWXLlrtUtWGOK1X1qg+yAvpgDo/OwM9XjP0ph/d3B04Bt5L1zd6lwOBrzRsVFaXeSExM9Or9gchtPbutX1XrORhduHBBBw4cqIA+/fTTmpmZ6VXPQIrmkqvXvLSCqrbJbZ2I/FNEKqvqSRGpTM775k8Ae1U13fOeFUAMUHivIWqMMX5SrFgx5syZQ5kyZXjzzTc5deoUffr0KZi5vHz/KmAg8Krnz5wOOe8EyolIRVX9HmgFpHg5rzHGBI0iRYowbdo0ypUrx4svvsixY8do0aIFRYsW9ek83gb+q8AiERkMfAn0BBCRhsBjqjpEVS+JSBywwbOffxfwNy/nNcaYoCIiTJw4kbJly7Jv3z6fhz14Gfiq+gPQOofXU8g6UHt5+WPgLm/mMsYYNxg1alSBfSnLvmlrjDEuYYFvjDEuYYFvjDEuYYFvjDEuYYFvjDEuYYFvjDEuYYFvjDEuYYFvjDEu4dXVMguSiHxP1rd386sC8C8flRMo3Naz2/oF69ktvOn5FlWtmNOKQhv43hKRFM3tEqFBym09u61fsJ7doqB6tl06xhjjEhb4xhjjEsEc+H91ugAHuK1nt/UL1rNbFEjPQbsP3xhjzH8K5i18Y4wx2VjgG2OMSwR04ItIBxE5IiLHRGRsDutvEJEPPOuTRKSGA2X6VB56flpE0kRkv4hsEJFbnKjTl67Vc7Zx3UREPXdcC2h56VlEenp+1qkistDfNfpaHv5uVxeRRBHZ4/n73dGJOn1FROaIyHcicjCX9SIib3n+e+wXkQZeT5rb3c0L+wMoCnwG3AqUAPYB4VeMGQb8xfO8N/CB03X7oeeWwG88z3/vhp4940oDm4EdQEOn6/bDz7k2sAco71m+2em6/dDzX4Hfe56HA184XbeXPTcHGgAHc1nfEfg7IEAMkOTtnIG8hR8NHFPVdFU9DyQAna8Y0xmY53m+BGjtua9uoLpmz6qaqKq/eBZ3AFX9XKOv5eXnDPAi8Brwqz+LKyB56flRYIaq/gSgqt/5uUZfy0vPCpTxPC8LfOPH+nxOVTcDP15lSGcgXrPsAMqJSGVv5gzkwK8CfJVt+YTntRzHqOpF4BRwk1+qKxh56Tm7wWRtIQSya/bs+VW3mqqu8WdhBSgvP+c6QB0R2SoiO0Skg9+qKxh56XkC0E9ETgAfAY/7pzTHXO//79fk1U3MTeElIv2AhsB9TtdSkESkCPAmMMjhUvytGFm7dVqQ9VvcZhGJUNWfnSyqgPUB5qrqFBFpDLwnIneqaqbThQWKQN7C/xqolm25que1HMeISDGyfg38wS/VFYy89IyItAGeA2JV9Zyfaiso1+q5NHAnsElEviBrX+eqAD9wm5ef8wlglapeUNXPgU/J+gcgUOWl58HAIgBV3Q6EkHWRsWCVp//fr0cgB/5OoLaI1BSREmQdlF11xZhVwEDP8+7ARvUcDQlQ1+xZRO4G3iEr7AN9vy5co2dVPaWqFVS1hqrWIOu4RayqpjhTrk/k5e/2CrK27hGRCmTt4kn3Y42+lpeejwOtAUSkHlmB/71fq/SvVcAAz9k6McApVT3pzQcG7C4dVb0oIiOAdWQd4Z+jqqkiMhFIUdVVwGyyfu07RtbBkd7OVey9PPb8OhAKLPYcnz6uqrGOFe2lPPYcVPLY8zqgnYikAZeA0aoasL+95rHnUcDfROQpsg7gDgrkDTgReZ+sf7QreI5LjAeKA6jqX8g6TtEROAb8Ajzi9ZwB/N/LGGPMdQjkXTrGGGOugwW+Mca4hAW+Mca4hAW+Mca4hAW+Mca4hAW+Mca4hAW+Mca4xP8DQoHrp2VsaNwAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def func(x):\n", " return np.exp(-x)-x**2\n", "def df(x):\n", " return -2*x - np.exp(-x)\n", "\n", "x0=0.0\n", "x1=1.0\n", "x = np.linspace(x0, x1, 100)\n", "y = func(x)\n", "plt.plot(x, y, color = 'k')\n", "plt.plot([x0,x1],[0,0])\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7034674224983918" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.optimize import fsolve\n", "x0 = fsolve(func, 0.0)[0]\n", "x0" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " x1 x2 f1 f2\n", " +0.0000000000 +1.0000000000 +1.0000000000 -0.6321205588\n", " +0.5000000000 +1.0000000000 +0.3565306597 -0.6321205588\n", " +0.5000000000 +0.7500000000 +0.3565306597 -0.0901334473\n", " +0.6250000000 +0.7500000000 +0.1446364285 -0.0901334473\n", " +0.6875000000 +0.7500000000 +0.0301753280 -0.0901334473\n", " +0.6875000000 +0.7187500000 +0.0301753280 -0.0292404858\n", " +0.7031250000 +0.7187500000 +0.0006511313 -0.0292404858\n", " +0.7031250000 +0.7109375000 +0.0006511313 -0.0142486319\n", " +0.7031250000 +0.7070312500 +0.0006511313 -0.0067872536\n", " +0.7031250000 +0.7050781250 +0.0006511313 -0.0030651888\n", " +0.7031250000 +0.7041015625 +0.0006511313 -0.0012063109\n", " +0.7031250000 +0.7036132812 +0.0006511313 -0.0002774104\n", " +0.7033691406 +0.7036132812 +0.0001869053 -0.0002774104\n", " +0.7033691406 +0.7034912109 +0.0001869053 -0.0000452413\n", " +0.7034301758 +0.7034912109 +0.0000708348 -0.0000452413\n", " +0.7034606934 +0.7034912109 +0.0000127975 -0.0000452413\n", " +0.7034606934 +0.7034759521 +0.0000127975 -0.0000162218\n", " +0.7034606934 +0.7034683228 +0.0000127975 -0.0000017121\n", " +0.7034645081 +0.7034683228 +0.0000055427 -0.0000017121\n", " +0.7034664154 +0.7034683228 +0.0000019153 -0.0000017121\n", " +0.7034673691 +0.7034683228 +0.0000001016 -0.0000017121\n", "\n" ] } ], "source": [ "x1, x2 = 0.0, 1.0\n", "f1, f2 = func(x1), func(x2)\n", "print('%+15s %+15s %+15s %+15s' % ('x1','x2','f1','f2'))\n", "print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))\n", "\n", "list_bisec = [[0],[abs(x1-x0)]]\n", "for i in range(0, 20):\n", " x = (x1 + x2)/2\n", " f = func(x)\n", " if (f*f1>=0.0):\n", " x1, f1 = x, f\n", " list_bisec[0].append(i)\n", " list_bisec[1].append(abs(x1-x0))\n", " else:\n", " x2, f2 = x, f\n", " list_bisec[0].append(i)\n", " list_bisec[1].append(abs(x2-x0))\n", "\n", " print('%+15.10f %+15.10f %+15.10f %+15.10f' % (x1,x2,f1,f2))\n", "\n", "list_bisec\n", "print()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.7182818284590451" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df(-1.0)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0000000000 -0.6321205588285576659757226\n", "0.7330436052 -0.0569084480040253914978621\n", "0.7038077863 -0.0006473915387465445370196\n", "0.7034674683 -0.0000000871660306156485376\n", "0.7034674225 -0.0000000000000014988010832\n", "\n" ] } ], "source": [ "x1 = 1.0\n", "f1 = func(x1)\n", "list_newton = [[0],[x1]]\n", "print('%-15.10f %+24.25f' % (x1,f1))\n", "for i in range(0, 4):\n", " x1 = x1 - f1 / df(x1)\n", " f1 =func(x1)\n", " print('%-15.10f %+24.25f' % (x1,f1))\n", " list_newton[0].append(i)\n", " list_newton[1].append(abs(x1-x0))\n", "\n", "list_newton\n", "print()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqfklEQVR4nO3deXxU5dn/8c+VjUDClkxYEwiBhEXcILKoLFahCCi2WhW3Kij6VKpd7E99bJ+2T9sHq9W2IpUqIkqVpWgtClZUKCgCsgiChCWELWzZ2JJA1uv3x0xojElMmOVMMtf79cprMvec5eIwme+cc+5zH1FVjDHGhKYwpwswxhjjHAsBY4wJYRYCxhgTwiwEjDEmhFkIGGNMCItwuoDGcLlcmpyc7HQZxhjTpGzcuDFPVRNqe61JhUBycjIbNmxwugxjjGlSRGR/Xa/Z4SBjjAlhFgLGGBPCLASMMSaEWQgYY0wIsxAwxpgQ5mjvIBGJAf4ClAL/VtXXnazHGGNCjc/3BERktojkiMi2Gu1jRGSniGSKyGOe5u8Ci1T1PuB6X9dijDGmfv44HDQHGFO9QUTCgRnAtUA/YKKI9AMSgYOeySr8UItbYS4s/x0c2+63VRhjTFPk8xBQ1VVAQY3mQUCmqmapaikwH5gAZOMOgjprEZEpIrJBRDbk5uaeX1HFebDqKcjbeX7zG2NMMxWoE8Nd+c83fnB/+HcF3gJuFJEXgHdqm1FVX1TVdFVNT0io9arnb1RQVArAoRNnz2t+Y4xprhw9MayqRcA9/l7PqbNlxAFHT56lq79XZowxTUig9gQOAUnVnid62owxxjgoUCGwHkgVkR4iEgXcCiwO0LqNMcbUwR9dROcBa4DeIpItIpNVtRyYCrwPZAALVfVLX6/bGGNM4/j8nICqTqyjfSmw1NfrM8YYc/5s2AhjjAlhIRUCpeX+ux7NGGOaopAIgTbRkQC8se4Aj7+1lSMnzzhckTHGBIeQCIG4mCgArkx1sWjjQUY8/W9+8+528gtLHK7MGGOcFRIhUOWWy5JY/tORTLi4C6+s3suwp1bwzLKdnDxT5nRpxhjjiJAKAYCkuFY8/b2L+eAnI/hWnw5MX57JsN8vZ8aKTIpLy50uzxhjAirkQqBKz4RYnr9tAEseupLLkuN4+v2dDH/q37yyei8lPj6BrKqcPluGqvp0ucYY4y1Hxw4KBhd0acvLd1/Gxv3Hefr9Hfz6ne3M+ngvD13dixsHJBIR3ricLC2vZE9uIdsPnyLjyCkyjp4i48hpCopKGZbq4olxfenTqY2f/jXGGNM4IR8CVQZ2b8+8+4awOjOfp5ft5NE3tzJzZRY/HpXG+As7ExYmX5vneFEpGUdOsd3zk3HkNJk5pymrcH/jbxERRu9OrRnVtyPxsVG8vu4AY//8MbdclsSPR6XRoXV0oP+ZxhjzFRYC1YgIV6a6uKJXPB9m5PDMsp08NO9z/rIikx9+KxVF3d/uj5xm++FTHD31n6GpO7RuQd/ObRiRlkC/Lm3o17k1yfExX9mTmDI8henLM3n1030s3nyYH1zVi8lX9iA6MtyJf64xxiBN6Th1enq6btiwofEz5mTAX4bA9+bABd9p8GyVlco7Xxzmjx/sYl9+MQARYUKvDrH07dyGvp1bex7b4Ipt0eDl7s0r4sn3Mnj/y2N0aRvNo9f24bqLutS6t2GMMd4SkY2qml7ba7YnUI+wMGHCJV0Ze2FnPt2TT3xMFKkdY2kR4d039x6uGP56Zzprs/L57ZLtPDx/M7NX7+MX4/qSnhzno+qNMeabhVbvoPPc64kMD2NEWgL9u7b1OgCqG5ISz+IHr+SZ713M0ZNnuGnmGh58fRMHPHsdxhjjb6ERArEdIaIlbHvT6Uq+JixMuHFgIiseGcmPr0lj+Y4crnl2Jf+3NMMuYjPG+F1ohECrOLjqcdjxLmwPznvZtIqK4OFrUvn3z0Yy4ZIuvPRxFiOfXsFra/ZRVlHpdHnGmGYqNE4MA1SUw0sjoTAXpn4G0W19WpuvbTt0kt8tyWBNVj49E2J4/Nq+XJTUloiwMMLDhIgw+cqjiJ1UNsbUrr4Tw46GgIjcAIwD2gAvq+qy+qb3KgQADm2CWVfDwLth/B/PfzkBoqp8mJHDtKUZZOUV1TtteI1QcD+GnXveMiqcX4zvx4i0hABVb4wJFn4JARGZDYwHclS1f7X2McCfgXBglqo+2YBltQf+oKqT65vO6xAA+Nd/w9oZcM+/oPtQ75YVIGUVlSz78hjHi0upqFTKK5WKykr3Y0XV8xrtnueVnseN+49TUFTKkoeuJLF9K6f/ScaYAPJXCAwHCoHXqkJARMKBXcAoIBv3DeYn4g6EaTUWMUlVczzzPQO8rqqb6lunT0KgpBD+MhQiW8IDH0NEw/v3N2X78ooYP/0TUjvGsvD+oUQ2cjgMY0zTVV8InPcngaquAgpqNA8CMlU1S1VLgfnABFXdqqrja/zkiNvvgffqCgARmSIiG0RkQ25u7vmW+x8tYmHcM5C3Ez75k/fLayKSXTE8eeOFfH7gBH94f6fT5RhjgoSvvw52BQ5We57taavLD4FrgJtE5IHaJlDVF1U1XVXTExJ8dDw7bTT0vxE+/gPk7vLNMpuA8Rd14fbB3fjrqiyW7zgWsPUWFJUyec56frpwC29tyuZYteE2jDHOcvSKYVV9DnjOkZWPeRIyP4J3Hoa7l0BYaBwe+cX4fmzcf5yfLtzC0oeH0bltS7+ur7CknHte+YyMo6dpFRXOm5uyAeiZEMMVvVxc3tPF0JR42raK9Gsdxpja+ToEDgFJ1Z4netqCT2wHGP1bWDwVPn/N3WMoBERHhjPj9gFcN/0THpr3OfPuG9Lo4bIbqqS8ggfmbmTb4VPMvGMgV/fpwPYjp/h0Tx6rM/P5+4ZsXluzHxHo36Utl/eK54qeLi5LjqNllA2qZ0wgeNVFVESSgXernRiOwH1i+GrcH/7rgdtU9UvvS/XRieHqVOHV6+DIF+5rB1p38t2yg9zbnx/iRws28+BVPfnZt/v4fPkVlcpD8z5nydYj/OF7F3PTwMSvTVNaXsnmgydYnZnHp3vy+PzACcorlajwMC7t1o4rerlHdL0osZ2dyDbGC/7qHTQPGAm4gGPAL1X1ZREZC/wJd4+g2ar6u/NaQS18HgIAeZnwwuXQ+1q4+VXfLjvIPbroCxZuPMir9wxiuA+vH1BVnnh7G2+sO8ATY/ty3/CUBs1XVFLO+n0FfLonn9WZeWw/cgpViIkKZ3BKPMNSXYxIS6CHK8YujjOmEYL2YrHG8ksIAKx6Gpb/FibOd4dBiDhTWsENM1aTV1jC0oeH0bGNb25y88yynUxfnskDI3ry2LXnv5dxvKiUNVnuQFidmXduOO/E9i0ZkZbAiLQELu/lIraFDYZrTH0sBL5JeSm8OALOnoIH10KL1r5fR5DKzDnNddNXc3FSW16/dwjhXt7T4JXVe/n1O9u5JT2JJ2+80Kff2A/kF7Nydy4rd+ayZk8eRaUVRIQJ6cntGe4JhX6d29hegjE1WAg0xMHP4OXRMPh+uPb3/llHkFq0MZtH/r6Fh65O5Sej0s57OVXnGUb368hfbh/gtxPO4D6fsHH/cVbuymXlrlwyjpwCIKF1C4anJjCidwLDerloHxPltxqMaSosBBpqySOwfhbc+xEkDvTfeoLQTxdu4a3Ps3l98mAu7+Vq9PwrduZw36sbSE9uz5x7BgX8lpk5p86yanceK3fl8vHuXE4UlyECFyW2O3fo6JKkdl7v6RjTFFkINNTZUzBjMLRsD/evhPDQ6bteXFrOddM/4dTZcpY+NIyE1g0fTmPj/gJun7WOXh1imXffEFpHO7vdKiqVL7JPsGpXHit35bD54AkqFVyxLRjTvyNj+3dmUI84v+6pGBNMLAQaY8cSmH8bXP1LGPYT/64ryOw4eooJz6/msuQ4Xp00qEHfmnccPcXNM9cQH9uCvz8wtFH3Wg6UE8WlrNqdx/vbjrJ8Rw5nyiqIj4li9AWdGHthJ4akxFsXVNOsWQg01oI7YPcH8F+fQnxP/68viMz/7ACPvbWVn45K44dXp9Y77cGCYm584VNEYNEDl5MUF/yjk54prWDlrhyWbD3K8oxjFJVW0K5VJKP7dWTshZ25vKeLqAgLBNO8WAg01qkjMGMQdLkU7vonhFBvE1XlRws2886Ww7xx3xCGpMTXOl3u6RK+N/NTjheXsfD+ofTu1PR6VJ0tq2DVrlze23aUD7cf43RJOW2iIxjVz72HcGWqy6f3lDbGKRYC52P9y7DkJ3DDC3DJbYFZZ5AoLCnn+umfUFTqPj8QX+MQz6mzZUx8cS1ZuUX87d7BDOze3qFKfaekvIJPduexdOtRlm0/yumz5bRuEcE1/Tpybf9ODE9LCPjJbmN8xULgfFRWwivXuoecnroBYhrfY6Yp2374FDf8ZTVDU+J55e7LCPOcHzhbVsH3Z3/Gxv3HmfX9dEb27uBwpb5XWl7J6j15vLf1CMu2H+NEcRkxUeF8Lz2JB6/q1aiT5sYEAwuB85WzA2ZeCRd8B258KXDrDRJ/W7ufn7+9jUfH9OG/RvakvKKSH7y+iWXbj/HnWy9hwiX1jRLePJRVVLJmTz5vf36If245TFR4GHdfkcz9w1No18quQTBNg4WAN1b8H6z8PdzxJvS6JrDrdpiqMnXe5/xr21EWTBnCwg0HWbghm19d14+7r+jhdHkBl5VbyJ8+3M07XxwmNiqC+4anMOnKHjZshQl6FgLeKC+BF66AihL4wVqIigns+h12+mwZ46d/wtGTZykpr/T6quLmYMfRUzyzbBcfbD9GXEwU/zWiJ3cO7W7nDEzQ8svtJUNGRAu47k9w4gBsfsPpagKudXQkz08cQJgIdw3tzo+vqb/baCjo06kNL92Vzj8fvIILurThd0szGP7UCuau2UdpeaXT5RnTKLYn0FDT06FNF/j+YmfW77AzpRV2o5c6rMvK5w/LdrJ+33ES27fk4atT+c6lXe2KZBM0bE/AF/peB/s+geICpytxhAVA3QanxLPw/qHMuecy2reK4meLvmD0n1bx7heHqaxsOl+yTGiyEGiovuNBK2DX+05XYoKQiDCydwcWT72CmXcMJCJMmPrG54yb/gkfZRyjKe1xm9DieAiISIyIbBCR8U7XUq8uA6BNV8h4x+lKTBATEcb078R7Dw/nT7dcQnFpOZNf3cB3X/iUl1ZlsezLo+w8epozpRVOl2oM4MWN5kVkNjAeyKm6x7CnfQzwZ9y3l5ylqk9+w6IeBRaebx0BIwJ9xsOmV6G0KOR6CZnGCQ8Tbri0K+Mu6sybG7OZ8e9Mfrc04yvTdGoTTbf4ViTHt6J7fAzJ8TF0j29F9/hWjo/EakKHN/cYHg4UAq9Vu9F8OO4bzY8CsnHfaH4i7kCYVmMRk4CLgXggGshT1XfrW6ejJ4YB9q5y35j+5rnQ73rn6jBN0sniMvYXFLEvv5gD+e7H/Z7H3NMlX5k2PibKExDuYEiOj+HKVFdQjtJqgl99J4bPe09AVVeJSHKN5kFApqpmeVY8H5igqtNw7zXULGwkEAP0A86IyFJVrawxzRRgCkC3bt3Ot1zf6HY5tIxzHxKyEDCN1LZVJBe1asdFie2+9lpRSTkHCv4TCvs9AfHZ3gLe3nwIVYiODOOuoclMGZ5iYWB8xteXOnYFDlZ7ng0MrmtiVX0CQETuxr0n8LVO1qr6IvAiuPcEfFlso4VHQO+x7hAoL4UIGzbA+EZMiwj6dm5D385tvvba2bIKMnMKmf3JXmZ9nMXcNfu5c2h3CwPjE46fGAZQ1TnfdCgoaPQdDyUnYd/HTldiQkR0ZDj9u7bl2Vsu4YOfjGBM/07M+jiLYb9fwbSlGeQVlnzzQoypg69D4BCQVO15oqet+Ui5CiJjrJeQcUTPhFj+WC0MXrIwMF7ydQisB1JFpIeIRAG3As3rEtvIaEgdBTuXuoebNsYBVWGw7Mcj+PYFHf8TBu9lkG9hYBrBm95B84CRgAs4BvxSVV8WkbHAn3D3CJqtqr/zTalB0DuoytZF8OZkmLQMutV5ysOYgMnMKeT55bv555bDREeEc9fl3ZkyLOVrNwRqiPKKSg4UFLMnt4jMnEIycwrZn1/EwOT2/PBbqTZqahNko4j62tmT8FRPGHw/fNtnGWeM1zJzCpm+fDeLtxymZWQ4dw1N5r5hPWoNg+LScvbkFLEn1/1BX/W4L7+Isor/fC50aN2Czu1asuXgCTq1iebn4/sy7sLOSAjddrWpsxDwh7/dBPm74aHNIXUPYtM0ZOacZvryzHNhcOfQ7iS1b3Xugz4rt4hDJ86cmz48TOge14qUhFh6dYilZ0IMvTrEkpIQS9uW7gvXNu4/zi/e3sb2I6e4speLX0+4gJ4JsU79E00jWAj4w8Y58M7D8MBq6NT/Gyc3xgmZOad57qNM3vniMKrQKiqcngn/+ZDv6fnQ7xbfihYR3zxIYEWl8re1+/nDsp2cLavgvmEpTP1WL1pF2SGiYGYh4A+FOfCHNBjxKFz1uNPVGFOvwyfOoEDnNtHn7hftjdzTJUx7L4O3Nh2ia7uW/M91/Rjdr6MdIgpSNpS0P8R2gG5DYUfTuLzBhLYu7VrStV1LnwQAQELrFjx78yUsvH8osS0iuH/uRu6Zs579+UU+Wb4JHAsBb/QdD8e2QUGW05UY44hBPeJ496Er+fm4vqzfW8CoP67i2Q92cbbMRkltKiwEvNHHMxxShu0NmNAVGR7GvcNSWP7ISMZc0InnPtrNqD+u5KOMY06XZhrAzgl4a+YwiGwJk5c5XYkxQeHTzDx+8c9t7Mkt4pq+Hfnldf1IimvV4PnPlFaQX1RCfmEpBUWl5BWWcKK4jKE94+nfta0fK2++/DKKqPHoex2s+D84fRRad3K6GmMcd3kvF+89PJzZq/fy5w93c82zK5l6VS9G9E4gv6jU8+Hu/pDPL3J/0OcXlpx77Uwdh5Iiw4VfX9+f2wY7PJpwM2N7At46th1eGArjnoXLJjtdjTFB5fCJM/x2yXaWbj36tdeiIsKIj4kiLiaK+NgWxMdEuZ/HRuGKaUFctd8jI4RH39zKql253HpZEr+ecEGDurQaN+si6k+qMH0AtE+GO//hdDXGBKWN+49TUFRKXEwUrlj3B39si4hGdSmtqFSe/WAnM1bs4ZKkdrxwxwA6t23px6qbD+si6k8i7kNCe1fBmRNOV2NMUBrYvT2j+nVkYPf2dI+PoXV0ZKOvKQgPE3727T7MvGMAu4+d5rrpn7AuK99PFYcOCwFf6HMdVJbDrvedrsSYZm9M/868/eAVtImO5PZZ63hl9V6a0hGNYGMh4AtdB0JsJ9hh9xgwJhBSO7bm7alXMLJ3B379znZ+unCL365NqKhUVu7KZf2+gmYZNtY7yBfCwtwXjm1+A0qLIarh3eGMMeenTXQkL945kOdXZPLHD3ex89hpZt4xsFHdUeuTV1jCgvUHeWPdgXOD7XWPb8WNAxL57oCuJLZvHn/ndmLYV/asgLk3wC2vuwPBGBMwy3cc4+H5m4kIE6ZPHMCVqa7zWo6qsnH/ceau3c/SrUcoq1Cu6BXPHYO7c6asgkUbs/l0j/s8xOU947lpYCJj+ncK+gH0rHdQIFSUwdO9oPe18J2ZTldjTMjZm1fE/XM3kJlTyKNj+jBleEqDTz4XlZTz9uZDzF2znx1HT9O6RQQ3DkzkjiHd6dXhq8NlHywo5h+fH2LRxmwOFBQT2yKCcRd25qb0RNK7tw/KQfSCNgREJAz4DdAG2KCqr9Y3fVCHAMA/HoCd78HPMiE80ulqjAk5RSXl/GzRFpZuPcq4izrz9E0X1fstPTPnNH9be4A3N2ZzuqScvp3bcNfQ7ky4pMs3frtXVdbvO87fNxxkydYjFJdW0D2+FTcNSOS7AxPp2i54uq/6JQREZDYwHshR1f7V2scAf8Z9e8lZqvpkPcv4DnADkA8sUdWP6ltn0IdAxruw4Ha4823oeZXT1RgTklSVv67K4ql/7SC1Q2v+eudAkl0x514vq6jkg+3HmLtmP2uy8okKD2PshZ24c2h3BnQ7v2/yRSXl/GvbURZtzGZNVj4i1Q4XXdCZllHOXtjmrxAYDhQCr1WFgIiEA7uAUUA27hvPT8QdCNNqLGKS5+e4qv5VRBap6k31rTPoQ6C0GJ5KgUtvh3HPOF2NMSHt4925/HDe51RWKn++9VL6dWnDvM8OMO+zAxw7VULXdi25fUg3bk5PwnUe92Kuy8GCYt7adIhFmw5ysODMucNF9w3vQa8OrX22nsbw2+EgEUkG3q0WAkOBX6nqtz3PHwdQ1ZoBUDX/HUCpqi4UkQWqekst00wBpgB069Zt4P79+8+73oBYcAdkb4Afb3f3GjLGOOZgQTH3z91IxtFThItQXqmMSEvgrqHdGdm7A+E+ur9CbSorlfX7Cli0MZslW48QGR7GvPuG0K9LG7+tsy6BDIGbgDGqeq/n+Z3AYFWdWsf8rYDpQDGwQ1Vn1Le+oN8TANiyAP4xBSZ/CEmXOV2NMSHvTGkFzy3fTWWlctvgbnSPj/nmmXzsYEExN/91DaXllcyfMoTUjoHdIwjaYSNUtVhVJ6vqD78pAJqMtNEQFmEXjhkTJFpGhfPomD48PravIwEAkBTXitfvHUxYmHD7rHXsywueO7D5OgQOAUnVnid62kJHy/bQYzhkvOMeXM4YY4CUhFhev3cw5ZXKbS+t5WBBsdMlAb4PgfVAqoj0EJEo4FZgsY/XEfz6jHffcjInw+lKjDFBJK1ja+ZOHkRhSTm3z1rH0ZNnnS7p/ENAROYBa4DeIpItIpNVtRyYCrwPZAALVfVL35TahPQZB4jdhN4Y8zUXdGnLa5MHU1BUym2z1pJ7usTRes47BFR1oqp2VtVIVU1U1Zc97UtVNU1Ve6rq73xXahPSuhMkDYKM0NsJMsZ8s0uS2vHKPZdx5MRZ7nx5HceLSh2rxfow+kuf8XB0Kxzf53QlxpggdFlyHLO+n05WXhF3zl7HyTNljtRhIeAvVYPI7VjibB3GmKB1RS8Xf71jIDuPnubuVz6jsKQ84DVYCPhLXAp07O/uJWSMMXW4qk8Hpk8cwBfZJ5k0Zz1nSv1zX4S6WAj4U5/xcGAtFOY4XYkxJoiN6d+JZ2++mPX7Cpgyd4PfbpBTGwsBf+o7HlDYudTpSowxQW7CJV156saL+Hh3Hg++vonS8sqArNdCwJ869of2yXZIyBjTIN9LT+I3N/Tnox05PDz/c8or/B8EFgL+JOI+JJS1Es6edLoaY0wTcOeQ7vx8XF/e23aUR/6+hYpK/448YCHgb32vg8oy2P2B05UYY5qIe4el8MjoNN7efJj/fmsrlX4MguC+MWZzkDgIYju6Lxy7sN7bJRhjzDlTv5XK2bJKnl+RSXRkGL+6/gK/3LrSQsDfwsKg91j4YiGUnYHI4LnlnDEmuP10dBpnyyqY9cleWkSG8/i1fXweBHY4KBD6joeyIti7yulKjDFNiIjwxLi+3DmkO21bRtqeQJPVbaj78cgXkPZtZ2sxxjQpIsL/TvDPoSCwPYHAiIqBtkmQt8vpSowxTZC/AgAsBALHlWohYIwJOhYCgeJKg7zdUBmYqwCNMaYhLAQCxZXqPjl8+rDTlRhjzDmOhoCIdBORt0Vktog85mQtfudKcz/aISFjTBDx5vaSs0UkR0S21WgfIyI7RSSzAR/sFwKLVHUScOn51tIknAuB3c7WYYwx1XjTRXQO8DzwWlWDiIQDM4BRQDawXkQWA+HAtBrzTwLWAotEZBIw14tagl9sR2jR1vYEjDFB5bxDQFVXiUhyjeZBQKaqZgGIyHxggqpOA8bXXIaIPAL80rOsRcArtUwzBZgC0K1bt/Mt13ki1kPIGBN0fH1OoCtwsNrzbE9bXf4FPCQiM4F9tU2gqi+qarqqpickJPisUEe40iDXQsAYEzwcvWJYVbcBoTOqmisVtrzhHlY6uq3T1RhjjM/3BA4BSdWeJ3raDFQ7OZzpbB3GGOPh6xBYD6SKSA8RiQJuBRb7eB1Nl3UTNcYEGW+6iM4D1gC9RSRbRCarajkwFXgfyAAWquqXvim1GYjrAWERFgLGmKDhTe+giXW0LwXszuq1CY+EuBQLAWNM0LBhIwLNlWYhYIwJGhYCgeZKhYIsqChzuhJjjLEQCDhXGlSWw/F9TldijDEWAgHn6u1+tENCxpggYCEQaK5e7kcLAWNMELAQCLTothDbyUYTNcYEBQsBJ7hSIXen01UYY4yFgCOqbjWp6nQlxpgQZyHgBFcalJyEwhynKzHGhDgLASck2BhCxpjgYCHgBBtIzhgTJCwEnNC6C0TGWA8hY4zjLAScEBbmvl4gz3oIGWOcZSHglKoeQsYY4yALAae40uDkQSgtcroSY0wIsxBwStXJ4Xy71aQxxjkBCwERSRGRl0VkUbW2GBF5VUReEpHbA1VLUDjXQ8gOCRljnNOgEBCR2SKSIyLbarSPEZGdIpIpIo/VtwxVzVLVyTWavwssUtX7gOsbVXlTF5cCEmbdRI0xjmro7SXnAM8Dr1U1iEg4MAMYBWQD60VkMRAOTKsx/yRVre3y2ERgq+f3ioaX3QxERkO77jaGkDHGUQ0KAVVdJSLJNZoHAZmqmgUgIvOBCao6DRjfwPVn4w6CzdSxVyIiU4ApAN26dWvgYpsI6yFkjHGYN+cEugIHqz3P9rTVSkTiRWQmcKmIPO5pfgu4UUReAN6pbT5VfVFV01U1PSEhwYtyg1BCmvvEcGVo7QQZY4JHQw8HeU1V84EHarQVAfcEqoag40qDihI4cQDiejhdjTEmBHmzJ3AISKr2PNHTZhrKeggZYxzmTQisB1JFpIeIRAG3Aot9U1aIsIHkjDEOa2gX0XnAGqC3iGSLyGRVLQemAu8DGcBCVf3Sf6U2Q63ioFW8jSFkjHFMQ3sHTayjfSmw1KcVhRrrIWSMcZANG+E0V5odDjLGOMZCwGmuNCjOh6J8pysxxoQgCwGnnRtIzg4JGWMCz0LAaa5U96MdEjLGOMBCwGntukF4CxtDyBjjCAsBp4WFQ3wv6yFkjHGEhUAwSLAeQsYYZ1gIBANXGpzYD2Vnna7EGBNiLASCgSsNtBIKspyuxBgTYiwEgoH1EDLGOMRCIBjE93I/2slhY0yAWQgEg6gYaNvNBpIzxgSchUCwcKXa4SBjTMBZCASLqtFEKyudrsQYE0IsBIKFKxXKiuH0YacrMcaEEAuBYGF3GTPGOCCgISAiKSLysogsqtZ2g4i8JCILRGR0IOsJKna/YWOMAxocAiIyW0RyRGRbjfYxIrJTRDJF5LH6lqGqWao6uUbb26p6H/AAcEtjim9WYjtAdFsbSM4YE1ANur2kxxzgeeC1qgYRCQdmAKOAbGC9iCwGwoFpNeafpKo59Sz/555lhSYRu8uYMSbgGhwCqrpKRJJrNA8CMlU1C0BE5gMTVHUaML4hyxURAZ4E3lPVTQ2tp1lypUHmR05XYYwJId6eE+gKHKz2PNvTVisRiReRmcClIvK4p/mHwDXATSLyQC3zTBGRDSKyITc318tyg5wrFQqPwtmTTldijAkRjTkc5DVVzcd97L9623PAc/XM8yLwIkB6err6tUCnnTs5nAmJA52txRgTErzdEzgEJFV7nuhpM+fD1dv9aOcFjDEB4m0IrAdSRaSHiEQBtwKLvS8rRLXvDmGRNoaQMSZgGtNFdB6wBugtItkiMllVy4GpwPtABrBQVb/0T6khIDwS4lLsWgFjTMA0pnfQxDralwJLfVZRqLOB5IwxAWTDRgQbV5r7DmMVZU5XYowJARYCwcaVBpXlcHyf05UYY0KAhUCwSbCB5IwxgWMhEGziPfcbtjGEjDEBYCEQbKLbQOvO1kPIGBMQFgLByHoIGWMCxEIgGFXdalKb9ygZxhjnWQgEI1calJyEwvpG3jbGGO9ZCAQju9WkMSZALASC0bkQsB5Cxhj/shAIRm26QGSM9RAyxvidhUAwErEeQsaYgLAQCFZVPYSMMcaPLASCVUIanDwIpUVOV2KMacYsBIJV1cnh/Exn6zDGNGsWAsGqKgRy7byAMcZ/AhYCIpIiIi+LyKIa7TEiskFExgeqliYhLgUkzE4OG2P8qkEhICKzRSRHRLbVaB8jIjtFJFNEHqtvGaqapaqTa3npUWBhw0sOEREtoH2yhYAxxq8aenvJOcDzwGtVDSISDswARgHZwHoRWQyEA9NqzD9JVb82BoKIjAK2A9GNrjwUWA8hY4yfNSgEVHWViCTXaB4EZKpqFoCIzAcmqOo0oKGHdkYCMUA/4IyILFXVygbO2/y50mDPCqisgLBwp6sxxjRD3pwT6AocrPY829NWKxGJF5GZwKUi8jiAqj6hqj8C3gBeqi0ARGSK55zBhtzcXC/KbYJcaVBRAicOOF2JMaaZaujhIK+paj7wQB2vzalnvheBFwHS09NDa2zl6gPJxfVwthZjTLPkzZ7AISCp2vNET5vxFZfnVpN2ctgY4yfehMB6IFVEeohIFHArsNg3ZRkAWsVBK5eFgDHGbxraRXQesAboLSLZIjJZVcuBqcD7QAawUFW/9F+pIcp6CBlj/KihvYMm1tG+FFjq04rMVyWkQcY7TldhjGmmbNiIYOdKg+J8KMp3uhJjTDNkIRDs7FaTxhg/shAIdtZDyBjjRxYCwa5tEkREWwgYY/zCQiDYhYVDfKr1EDLG+IWFQFNg9xs2xviJhUBT4EqDE/uh7KzTlRhjmhkLgabAlQpaCQV7nK7EGNPMWAg0BdZN1BjjJxYCTUF8L0gbA9HtnK7EGNPMBGwoaeOFqFZw2wKnqzDGNEO2J2CMMSHMQsAYY0KYhYAxxoQwCwFjjAlhFgLGGBPCLASMMSaEWQgYY0wIsxAwxpgQJqrqdA0NJiK5wH4vFuEC8nxUjj9Yfd6x+rxj9XknmOvrrqoJtb3QpELAWyKyQVXTna6jLlafd6w+71h93gn2+upih4OMMSaEWQgYY0wIC7UQeNHpAr6B1ecdq887Vp93gr2+WoXUOQFjjDFfFWp7AsYYY6qxEDDGmBDW7EJARMaIyE4RyRSRx2p5vYWILPC8vk5EkgNYW5KIrBCR7SLypYg8XMs0I0XkpIhs9vz8T6Dqq1bDPhHZ6ln/hlpeFxF5zrMNvxCRAQGsrXe1bbNZRE6JyI9qTBPQbSgis0UkR0S2VWuLE5EPRGS357F9HfN+3zPNbhH5fgDre1pEdnj+//4hIu3qmLfe94If6/uViByq9n84to556/1792N9C6rVtk9ENtcxr9+3n9dUtdn8AOHAHiAFiAK2AP1qTPMDYKbn91uBBQGsrzMwwPN7a2BXLfWNBN51eDvuA1z1vD4WeA8QYAiwzsH/76O4L4RxbBsCw4EBwLZqbU8Bj3l+fwz4fS3zxQFZnsf2nt/bB6i+0UCE5/ff11ZfQ94LfqzvV8AjDfj/r/fv3V/11Xj9GeB/nNp+3v40tz2BQUCmqmapaikwH5hQY5oJwKue3xcBV4uIBKI4VT2iqps8v58GMoCugVi3j00AXlO3tUA7EensQB1XA3tU1ZuryL2mqquAghrN1d9nrwI31DLrt4EPVLVAVY8DHwBjAlGfqi5T1XLP07VAoq/X21B1bL+GaMjfu9fqq8/z2XEzMM/X6w2U5hYCXYGD1Z5n8/UP2XPTeP4ITgLxAamuGs9hqEuBdbW8PFREtojIeyJyQWArA0CBZSKyUUSm1PJ6Q7ZzINxK3X98Tm/Djqp6xPP7UaBjLdMEy3achHvPrjbf9F7wp6mew1Wz6zicFgzbbxhwTFV31/G6k9uvQZpbCDQJIhILvAn8SFVP1Xh5E+7DGxcD04G3A1wewJWqOgC4FnhQRIY7UEO9RCQKuB74ey0vB8M2PEfdxwWCsi+2iDwBlAOv1zGJU++FF4CewCXAEdyHXILRROrfCwj6v6XmFgKHgKRqzxM9bbVOIyIRQFsgPyDVudcZiTsAXlfVt2q+rqqnVLXQ8/tSIFJEXIGqz7PeQ57HHOAfuHe7q2vIdva3a4FNqnqs5gvBsA2BY1WHyDyPObVM4+h2FJG7gfHA7Z6g+poGvBf8QlWPqWqFqlYCL9WxXqe3XwTwXWBBXdM4tf0ao7mFwHogVUR6eL4p3gosrjHNYqCqF8ZNwPK6/gB8zXP88GUgQ1WfrWOaTlXnKERkEO7/o0CGVIyItK76HfcJxG01JlsM3OXpJTQEOFnt0Eeg1PkNzOlt6FH9ffZ94J+1TPM+MFpE2nsOd4z2tPmdiIwB/h9wvaoW1zFNQ94L/qqv+jmm79Sx3ob8vfvTNcAOVc2u7UUnt1+jOH1m2tc/uHuu7MLda+AJT9v/4n6zA0TjPoSQCXwGpASwtitxHxb4Atjs+RkLPAA84JlmKvAl7p4Oa4HLA7z9Ujzr3uKpo2obVq9RgBmebbwVSA9wjTG4P9TbVmtzbBviDqMjQBnu49KTcZ9n+gjYDXwIxHmmTQdmVZt3kue9mAncE8D6MnEfT696H1b1mOsCLK3vvRCg+uZ63ltf4P5g71yzPs/zr/29B6I+T/ucqvdctWkDvv28/bFhI4wxJoQ1t8NBxhhjGsFCwBhjQpiFgDHGhDALAWOMCWEWAsYYE8IsBIwxJoRZCBhjTAj7/7opiKrpNJqhAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "X = list_bisec[0]\n", "Y = list_bisec[1]\n", "plt.plot(X, Y)\n", "\n", "X = list_newton[0]\n", "Y = list_newton[1]\n", "plt.plot(X, Y)\n", "\n", "plt.yscale(\"log\") # y軸を対数目盛に\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## exp関数に関する注意\n", "exp関数のimport元で振る舞いが違うみたい．\n", "```python\n", "ValueError: sequence too large; cannot be greater than 32\n", "```\n", "がplot作成の前段階で出る．\n", "numpyでやるときには，np.exp(-x)などとしてる．\n", "\n", "でも，diffには通らない．そのあたり，覚悟して使う関数を決めないと．．．" ] }, { "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" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": { "height": "12px", "width": "252px" }, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": "block", "toc_window_display": true }, "vscode": { "interpreter": { "hash": "f3f87633aac09da3bda522f97956bee375b5501d1579e6458804e567301cb62a" } } }, "nbformat": 4, "nbformat_minor": 2 }