{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"Table of Contents
\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"
\n",
"\n",
"\n",
"数値計算プレ試験解答例(2015)\n",
"
\n",
"
\n",
"\n",
" 18/12/14 関西学院大学 西谷滋人 \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 数値解の収束性:25点\n",
"\n",
"次の関数\n",
"$$\n",
"f(x)=-4\\,\\exp(-x)+2\\,\\exp(-2\\, x)\n",
"$$\n",
"は$x=-1..0$に解$-\\ln(2)$を持つ．二分法によって数値解を求めよ．繰り返しは１０回程度でいい．また，収束の様子を片対数(logplot)でプロットせよ．"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FXXe/vH3NyEkIaGFcqRIUzorssEA8lvKonRE3FWExayKBiyAS6gLSMQWMPBQlgdFRUCRsEtRaYKuQcBCswBS3NBBiooIgaXm+/sjkcdCCTmTM6fcr+uai4TMzLk/Jt5M5syZY6y1iIhI8AhzO4CIiDhLxS4iEmRU7CIiQUbFLiISZFTsIiJBRsUuIhJkVOwiIkFGxS4iEmRU7CIiQaaQGw9aunRpW6VKlXxte/LkSWJiYpwN5Oc0c2jQzKHBm5k3bNjwnbW2zNXWc6XYq1Spwvr16/O17YoVK2jRooWzgfycZg4Nmjk0eDOzMWZPXtbTqRgRkSDjWLEbY8KNMZ8bYxY5tU8REbl2Th6x9wO2Org/ERHJB0eK3RhTEegAvOLE/kREJP+cOmIfDwwCsh3an4iI5JPx9o02jDEdgfbW2keNMS2AAdbajpdYLwlIAvB4PPHp6en5erysrCxiY2O9SBx4NHNo0MyhwZuZW7ZsucFa2/CqK1prvVqA54H9wG7gEHAKeONK28THx9v8ysjIyPe2gUozhwbNHBq8mRlYb/PQy16firHWDrXWVrTWVgHuBT6w1vbwdr+XkpGRwZtvvlkQuxYRCRoBdR37kiVLePXVV/nPf/7jdhQREb/laLFba1fYS5xfd8qAAQMoVKgQzz33XEE9hIhIwAuoI3aPx8Mdd9zB66+/zs6dO92OIyLilwKq2AG6du2qo3YRkSsIuGIvXbo0Dz/8MDNmzGD37t1uxxER8TsBV+wAgwcPJiwsjNTUVLejiIj4nYAs9ooVK/Lggw8ybdo09uzJ010sRURCRkAWO8DQoUMxxuhcu4jIrwRssVeqVImHH36YadOmsWvXLrfjiIj4jYAtdsg5ag8PD+fZZ591O4qIiN8I6GKvUKECvXr1Yvr06ezYscPtOCIifiGgix1gyJAhRERE8Mwzz7gdRUTELwR8sZcrV45HHnmEmTNn8vXXX7sdR0TEdQFf7JBz1B4dHc3IkSPdjiIi4rqgKPayZcvSr18/0tPT2bhxo9txRERcFRTFDjl3fixevDgjRoxwO4qIiKuCpthLlizJwIEDeeedd1izZo3bcUREXBM0xQ7Qr18/ypQpw/Dhw92OIiLimqAq9tjYWIYOHcr777/PBx984HYcERFXBFWxAzzyyCNcf/31DBky5Kc32xYRCSlBV+xRUVE89dRTrFu3jvnz57sdR0TE54Ku2AESExOpU6cOw4YN4/z5827HERHxqaAs9p9uDLZ9+3amT5/udhwREZ8KymIH6Ny5M40bNyYlJYX//ve/bscREfGZoC12YwyjR4/mwIEDTJgwwe04IiI+E7TFDtCsWTM6derE888/z3fffed2HBERnwjqYgdITU0lKytLt/UVkZDhdbEbY6KMMWuNMV8aY74yxjzlRDCn1KlTh549e/K///u/ejMOEQkJThyxnwH+aK2tD9wMtDXGNHZgv45JSUkhIiKCYcOGuR1FRKTAeV3sNkdW7qcRuYtfveSzfPnyJCcnM2fOHN0gTESCniPn2I0x4caYL4AjwHvWWr9rz4EDB+LxeOjfv79uNSAiQc04WXLGmBLAAqCPtXbzr76WBCQBeDye+PT09Hw9RlZWFrGxsfnadvHixaSlpZGSkkLz5s3ztQ83eDNzoNLMoUEzX5uWLVtusNY2vOqK1lpHF+BJYMCV1omPj7f5lZGRke9tz58/b3/3u9/ZqlWr2tOnT+d7P77mzcyBSjOHBs18bYD1Ng897MRVMWVyj9QxxkQDtwPbvN1vQQgPD2fs2LHs2rWLSZMmuR1HRKRAOHGOvRyQYYzZCKwj5xz7Igf2WyBuv/122rdvz9NPP823337rdhwREcc5cVXMRmttA2vtTdbaetbaUU4EK0hpaWmcPHmSJ5980u0oIiKOC/pXnl5K7dq1eeyxx5g6dSobN250O46IiKNCstgBRo4cSYkSJXjiiSd0+aOIBJWQLfa4uDiefvppMjIyWLBggdtxREQcE7LFDpCUlES9evVITk7m9OnTbscREXFESBd7oUKFGD9+PLt37yYtLc3tOCIijgjpYgdo1aoVf/rTn3juuefYu3ev23FERLwW8sUOMHbsWACSk5NdTiIi4j0VO1C5cmX+/ve/M3fuXP7973+7HUdExCsq9lwDBgygWrVq9OnTh3PnzrkdR0Qk31TsuaKiopgwYQJbt27Vm1+LSEBTsf9Mx44d6dSpEykpKezfv9/tOCIi+aJi/5UJEyaQnZ3N3/72N7ejiIjki4r9V6pWrcqwYcOYO3cuy5YtczuOiMg1U7FfwoABA6hRowaPP/64XpEqIgFHxX4JkZGRTJ48mczMTJ5//nm344iIXBMV+2XcdtttdO/endTUVLZt88s3hBIRuSQV+xWMGzeOIkWK0Lt3b93aV0QChor9CjweD2PGjOHDDz9kxowZbscREckTFftV9OzZk6ZNm5KcnKz3SBWRgKBiv4qwsDBeeuklTpw4Qf/+/d2OIyJyVSr2PKhbty5DhgzhjTfe4N1333U7jojIFanY82jYsGHUqlWL3r17k5WV5XYcEZHLUrHnUWRkJK+88gp79uxh+PDhbscREbksFfs1aNq0KY8++igTJ07k008/dTuOiMgleV3sxpjrjTEZxpgtxpivjDH9nAjmr55//nkqVKjAgw8+qNsNiIhfcuKI/TyQbK2tAzQGHjPG1HFgv36pWLFivPzyy2zdupWnn37a7TgiIr/hdbFbaw9aaz/L/fgEsBWo4O1+/Vnbtm25//77GT16NBs2bHA7jojILzh6jt0YUwVoAKxxcr/+aNy4cZQtW5YHHniAs2fPuh1HROQi49Q9UIwxscCHwLPW2vmX+HoSkATg8Xji09PT8/U4WVlZxMbGehPVMR999BHDhw8nMTGRBx54oMAex59m9hXNHBo087Vp2bLlBmttw6uuaK31egEigGVA/7ysHx8fb/MrIyMj39sWhB49etjw8HC7fv36AnsMf5vZFzRzaNDM1wZYb/PQsU5cFWOAV4Gt1tpx3u4v0EycOBGPx0NiYqKukhERv+DEOfamwH3AH40xX+Qu7R3Yb0AoWbIkr7zyClu2bGHkyJFuxxERceSqmNXWWmOtvclae3PussSJcIGiXbt2PPzww6SlpfHxxx+7HUdEQpxeeeqQsWPHUqlSJRITE3UvGRFxlYrdIUWLFmXGjBns3LmT5ORkt+OISAhTsTuoWbNmDBw4kKlTp7Jo0SK344hIiFKxO2zUqFHUr1+fnj176h2XRMQVKnaHRUZG8sYbb/Djjz/Ss2dPvQm2iPicir0A1KtXj9TUVBYuXMiLL77odhwRCTEq9gLSt29f2rRpQ//+/dmyZYvbcUQkhKjYC0hYWBjTp0+naNGidOvWTa9KFRGfUbEXoOuuu45p06axceNGhgwZ4nYcEQkRKvYC1rFjR/r06cOECRNYuHCh23FEJASo2H3ghRdeoEGDBtx///3s37/f7TgiEuRU7D4QGRlJeno6Z86coXv37pw/f97tSCISxFTsPlKjRg2mTJnCqlWrGDVqlNtxRCSIqdh96L777uOvf/0rzzzzDO+9957bcUQkSKnYfWzy5MnUqVOHv/zlL3zzzTduxxGRIKRi97GYmBj+9a9/cerUKe69916dbxcRx6nYXVC7dm1efPFFVq1axYgRI9yOIyJBRsXukh49epCUlERqaipvv/2223FEJIio2F00YcIE4uPjSUxMJDMz0+04IhIkVOwuioqKYt68eRQqVIg//elPnDp1yu1IIhIEVOwuq1y5Mm+++SabNm0iKSlJ928XEa+p2P1AmzZtGDVqFLNmzWLixIluxxGRAKdi9xN///vfufPOO0lOTiYjI8PtOCISwFTsfiIsLIwZM2ZQvXp17rnnHvbu3et2JBEJUCp2P1KsWDHeeustzp49S5cuXfRkqojkiyPFboyZZow5YozZ7MT+QlnNmjV58803+fzzz3nwwQf1ZKqIXDOnjtinA20d2lfI69ChA6mpqcyZM4fnnnvO7TgiEmAcKXZr7UrgqBP7khwDBw6kR48eDB8+nNWrV7sdR0QCiHHqV31jTBVgkbW23mW+ngQkAXg8nvj09PR8PU5WVhaxsbH5TBlYzp49S79+/di1axeTJk2ievXqbkfymVD6Pv9EM4cGb2Zu2bLlBmttw6ut57Ni/7mGDRva9evX5+txVqxYQYsWLfK1bSA6ePAgN998MxEREaxdu5by5cu7HcknQu37DJo5VHgzszEmT8Wuq2L8XLly5Xj22Wc5duwYnTt31pUyInJVKvYAcOONNzJ79mw2bNhAYmIi2dnZbkcSET/m1OWOs4FPgJrGmP3GmJ5O7Ff+T6dOnRg7dizz5s1j0KBBbscRET9WyImdWGu7ObEfubInnniCnTt3MnbsWKpWrcpjjz3mdiQR8UOOFLv4hjGG8ePHs2fPHvr27UulSpXo1KmT27FExM/oHHuACQ8PZ/bs2fz+97+na9eurFmzxu1IIuJnVOwBKCYmhkWLFlGuXDk6duzI119/7XYkEfEjKvYA5fF4ePfddwFo27Ythw4dcjmRiPgLFXsAq169OosXL+bw4cO0b9+e48ePux1JRPyAij3AJSQkMG/ePDZt2kTnzp05ffq025FExGUq9iDQtm1bZsyYwYoVK+jWrRvnz593O5KIuEjFHiS6d+/OxIkTeeutt0hKStKrU0VCmK5jDyJ9+vTh+++/56mnnqJYsWL8z//8D8YYt2OJiI+p2IPMyJEj+fHHHxk/fjzFihVj1KhRbkcSER9TsQcZYwzjxo3jxIkTPP300xQtWpSBAwe6HUtEfEjFHoSMMbz00kucPHmSQYMGERUVRZ8+fdyOJSI+omIPUuHh4cycOZPTp0/Tt29fChcuTK9evdyOJSI+oKtiglhERATp6el06NCB3r17M23aNLcjiYgPqNiDXGRkJHPnzqV169Y89NBDTJ8+3e1IIlLAVOwhICoqirfeeotWrVrx4IMP8tprr7kdSUQKkIo9RERHR/POO+9w22230bNnT52WEQliKvYQEh0dzdtvv83tt9/OQw89xNSpU92OJCIFQMUeYn4q93bt2tGrVy8mTZrkdiQRcZiKPQRFRUWxYMECunTpQt++fRkzZozbkUTEQSr2EFW4cGHmzJnDvffey+DBgxkxYgTWWrdjiYgD9AKlEBYREcEbb7xBTEwMzzzzzMV7zISF6d97kUCmYg9x4eHhvPzyyxQvXpxx48Zx/PhxXnnlFQoV0o+GSKDS/72CMYa0tDRKlCjBk08+yQ8//EB6ejrR0dFuRxORfHDkd25jTFtjzHZjTKYxZogT+xTfMsYwYsQIJk+ezMKFC2nTpg3Hjh1zO5aI5IPXxW6MCQcmA+2AOkA3Y0wdb/cr7nj00UdJT0/n008/pVmzZhw4cMDtSCJyjZw4Yk8AMq21O621Z4F0oLMD+xWX3HPPPSxdupRdu3bRpEkTvvrqK7cjicg1MN5e4maM+TPQ1lr7UO7n9wGNrLWPX26bhg0b2vXr11/zYz218Cs+3rKXEiVK5DtvIDp27JgrM2dlZbFp0yaysy9Qr149ihf3XQa3ZnaTZg4NxbKP8/IjbfK1rTFmg7W24dXW89mTp8aYJCAJwOPxsGLFimvex/79Z7hw4ULInft1c+YbbriBnTt38uWXG7n++uspWbKkTx5X3+fQEIozR0dfyFf/XRNrrVcL0ARY9rPPhwJDr7RNfHy8za+MjIx8bxuo3J75+++/t82aNbOAHTVqlM3Ozi7wx3R7Zjdo5tDgzczAepuHXnbiHPs6oLoxpqoxpjBwL/COA/sVPxEXF8fy5cvp0aMHTz75JPfffz9nzpxxO5aIXIbXp2KsteeNMY8Dy4BwYJq1Vs+2BZnIyEhmzpxJ9erVGTlyJLt27WLevHmUKVPG7Wgi8iuOXMdurV1ira1hrb3BWvusE/sU/2OM4cknn2T27NmsW7eOhIQENm/e7HYsEfkV3RRErtm9997Lhx9+yJkzZ2jSpAkLFy50O5KI/IyKXfIlISGBtWvXUrNmTTp37swzzzyju0OK+AkVu+RbxYoVWbVqFX/5y18YMWIEd999N1lZWW7HEgl5KnbxSnR0NDNnziQtLY0FCxbQqFEjtm/f7nYskZCmYhevGWNITk5m+fLlHDlyhFtuuYUFCxa4HUskZKnYxTGtWrViw4YN1KpVi7vuuovBgwdz/vx5t2OJhBwVuziqUqVKrFq1it69ezNmzBj++Mc/6g6RIj6mYhfHRUZGMmXKFGbNmsVnn31GgwYNWL58uduxREKGil0KTPfu3Vm3bh1ly5alTZs2DB06lHPnzrkdSyToqdilQNWuXZu1a9fy8MMPk5qaSvPmzdm9e7fbsUSCmopdClyRIkWYOnUq6enpfPXVV9SvX5/Zs2e7HUskaKnYxWe6du3KF198Qd26denevTuJiYkcP37c7VgiQUfFLj5VtWpVVq5cSUpKCrNmzaJ+/fqsXLnS7VgiQUXFLj5XqFAhRo4cyerVqwkPD6dFixYMHDiQ06dPux1NJCio2MU1TZo04YsvvqBXr16kpaXRsGFD8vNeuCLySyp2cVVsbCxTpkxh6dKlHDt2jMaNGzNs2DDOnj3rdjSRgKViF7/Qtm1bNm/ezH333cdzzz1Hr169WLNmjduxRAKSil38RokSJXjttddYtGgRWVlZ3HrrrSQnJ3Pq1Cm3o4kEFBW7+J0OHTowffp0kpKSGDduHHXr1uXdd991O5ZIwFCxi1+KiYlhypQpfPjhh0RGRtKuXTu6devGoUOH3I4m4vdU7OLXmjVrxpdffklKSgrz58+nVq1a/OMf/+DChQtuRxPxWyp28XuRkZGMHDmSjRs3csstt9CnTx8SEhL05KrIZajYJWDUrFmT5cuXk56ezsGDB2ncuDEPPPAAhw8fdjuaiF9RsUtAMcbQtWtXtm/fzsCBA5k1axY1atRg3LhxuvZdJJeKXQJS0aJFGTNmDJs2baJJkyYkJydTr149Fi5ciLXW7XgirvKq2I0xdxtjvjLGZBtjGjoVSiSvatasydKlS1m8eDHh4eHccccd3HbbbXz22WduRxNxjbdH7JuBuwDdnk9cY4yhffv2bNy4kUmTJvHll18SHx/Pfffdx969e92OJ+JzXhW7tXartXa7U2FEvBEREcHjjz/Ojh07GDJkCHPnzqVGjRokJyfz3XffuR1PxGeME+cjjTErgAHW2svems8YkwQkAXg8nvj09PR8PVZWVhaxsbH52jZQaeb8OXLkCNOnT2fZsmVERUVxzz338Oc//5mYmBiHUjpL3+fQ4M3MLVu23GCtvfppb2vtFRfgfXJOufx66fyzdVYADa+2r5+W+Ph4m18ZGRn53jZQaWbvbNmyxXbp0sUCNi4uzo4ePdpmZWU5tn+n6PscGryZGVhv89CxVz0VY629zVpb7xLL2/n6J0fEx2rXrs38+fNZt24djRo1YvDgwVSrVo20tDROnjzpdjwRx+lyRwkZDRs2ZMmSJXz00UfUr1+fgQMHUrVqVcaMGcOJEyfcjifiGG8vd+xijNkPNAEWG2OWORNLpODceuutLF++nI8++ogGDRowePBgKleuTEpKCkePHnU7nojXvL0qZoG1tqK1NtJa67HWtnEqmEhBu/XWW1m2bBlr1qyhefPmPPXUU1SqVIn+/fuzb98+t+OJ5JtOxUjIS0hIYMGCBWzatIkuXbowceJEqlWrRmJiIl9++aXb8USumYpdJFe9evV4/fXX2bFjB4899hjz58/n5ptvplWrVixevJjs7Gy3I4rkiYpd5FcqV67M+PHj2bdvH6NHj2b79u107NiRWrVqMWnSJD3RKn5PxS5yGSVLlmTQoEHs2rWLN998k1KlStG3b18qVKhAnz592LJli9sRRS5JxS5yFREREXTr1o1PPvmENWvWcOeddzJ16lTq1q1Ly5YtSU9P58yZM27HFLlIxS5yDRISEpg5cyb79+9n9OjR7Nmzh27dulGhQgUGDBjAtm3b3I4oomIXyY8yZcowaNAgMjMzWbZsGS1atGDChAnUrl2bpk2b8uqrr+pcvLhGxS7ihbCwMFq3bs3cuXPZt28fL7zwAkePHuWhhx7C4/HQo0cP3nvvPb35tviUil3EIddddx0DBgxgy5YtfPzxxyQmJrJ48WJat25NpUqVSE5O5rPPPtM7PEmBU7GLOMwYQ5MmTXjxxRc5ePAg//znP7nllluYNGkS8fHx1K5dm5SUFLZu3ep2VAlSKnaRAhQVFcXdd9/NW2+9xaFDh3jppZcoX748o0aNok6dOvzud79j1KhRunRSHKViF/GRuLg4kpKS+OCDDzhw4AATJ06kZMmSpKSkULduXRITExk6dChr167V6RrxiopdxAXlypWjT58+rFy5kgMHDjB58mTKlCnDCy+8QKNGjahYsSK9e/dmyZIlnD592u24EmBU7CIuK1euHI8++ihjx47lyJEjzJgxgyZNmjBr1iw6dOhAXFwcnTp1YsqUKezevdvtuBIACrkdQET+T1xcHImJiSQmJnLmzBkyMjJYsmQJixcvZtGiRQDUrFmTNm3a0Lp1a5o3bx5y7xkqV6cjdhE/FRkZSdu2bZk4cSKZmZls3bqV8ePHU7VqVaZOnUrHjh2Ji4ujefPmjBo1itWrV3P27Fm3Y4sfULGLBABjDLVq1aJfv34sXbqUH374gffff5/+/fuTlZVFSkoKf/jDHyhZsiStW7fm2WefZdWqVbqHTYjSqRiRABQVFUWrVq1o1aoVqampHD16lJUrV/LBBx+wYsUKhg8fDuQc9SckJNC0aVOaNm1K48aNKV26tMvppaCp2EWCQFxcHHfeeSd33nknAEePHmXVqlWsXr2a1atXk5aWRmpqKgDVq1enSZMmNGrUiISEBG666SYKFy7sZnxxmIpdJAjFxcXRuXNnOnfuDMCpU6dYv349n3zyCZ988gnvvvsuM2fOBHKO6uvXr098fDwNGzakQYMG1K1bV2UfwFTsIiGgSJEiNGvWjGbNmgFgrWXv3r2sWbOGtWvXsmHDBmbNmsWUKVOAnHvQ161bl5tvvpn69etTv3596tWrR5kyZdwcQ/JIxS4SgowxVK5cmcqVK3PPPfcAkJ2dTWZmJp9//vnFZenSpUyfPv3idh6Ph3r16lG3bl3q1KlDnTp1qFWrlgrfz6jYRQTIuQVxjRo1qFGjBl27dr3494cPH2bjxo1s2rSJzZs3s2nTJl599VVOnjx5cZ1SpUpRs2bNi9tXr16dG2+8kRtuuIGiRYu6MU5IU7GLyBV5PB5uv/12br/99ot/l52dzb59+9iyZQvbt29n27ZtbNu2jeXLl//iCB+gbNmyVKtWjWrVqlG1alWqVKlClSpVqFy5sq67LyBeFbsx5gWgE3AW2AE8YK095kQwEfFfYWFhF0/ltGvX7hdfO3HiBJmZmezYsYPMzEwyMzPZtWsXH3/8MXPmzPnNm46ULVuW66+/nooVK1KxYkUqVKhAhQoVKFeuHOXLl+e6664jLi4OY4wvRwxo3h6xvwcMtdaeN8aMBoYCg72PJSKBqmjRojRo0IAGDRr85mvnzp3jwIED7Nmzh927d7Ny5UrCw8PZt28fO3bs4MMPP+TYsd8eG0ZERODxePB4PJQtW5ayZctSpkwZSpcuTZkyZShVqhSlSpUiLi6OuLg4SpYsSWRkpC/G9UteFbu1dvnPPv0U+LN3cUQkmEVERFw8FdO8eXMqV65MixYtfrHOyZMnOXjwIN988w3ffPMNhw8f5tChQxw8eJAjR45w5MgRNm/ezLfffnvFO19GR0dTokQJSpQoQfHixSlevDjFihWjWLFiFC1alNjY2It/FilShJiYGGJiYihSpAjR0dFER0cTFRVFVFQUkZGRv1jCwvz7RftOnmN/EJjj4P5EJATFxMRw4403cuONN15xPWstp06d4ttvv+X777/n6NGjF//84Ycf+OGHHzh27BjHjh3jxx9/5OjRo+zZs4fjx49z/PhxTp48me/73oeFhVG4cGEiIiIoVKjQxSU8PPziEhYWRlhYGMaYi6eRjDH07t37N/+YOc1cbTBjzPvAdZf40jBr7du56wwDGgJ32cvs0BiTBCQBeDye+PT09HwFzsrKCrm72Wnm0KCZfctay5kzZ/jvf//L6dOnL/559uxZzpw5w+nTpzl37tzFz8+fP8+5c+c4d+4cFy5c4Ny5c5w/f54LFy78YsnOziY7OxtrLdZasrOzLz4ewF133cVNN92Ur8wtW7bcYK1tmKfhvFmA+4FPgCJ53SY+Pt7mV0ZGRr63DVSaOTRo5tDgzczAepuHjvX2qpi2wCCgubX2lDf7EhERZ3j7DMA/gKLAe8aYL4wxLzqQSUREvODtVTFXfnZDRER8zr+v2RERkWumYhcRCTIqdhGRIKNiFxEJMip2EZEgc9VXnhbIgxrzLbAnn5uXBr5zME4g0MyhQTOHBm9mrmytveq7mrhS7N4wxqy3eXlJbRDRzKFBM4cGX8ysUzEiIkFGxS4iEmQCsdinuh3ABZo5NGjm0FDgMwfcOXYREbmyQDxiFxGRK/D7YjfG3G2M+coYk22MuewzycaYtsaY7caYTGPMEF9mdJoxJs4Y854x5j+5f5a8zHp/y/1vs9kYM9sYE+XrrE65hplLGGPmGmO2GWO2GmOa+DqrU/I6c+664caYz40xi3yZ0Wl5mdkYc70xJsMYsyX357ufG1m9cbU+Mjkm5n59ozHm904+vt8XO7AZuAtYebkVjDHhwGSgHVAH6GaMqeObeAViCPBva2114N+5n/+CMaYC0BdoaK2tB4QD9/o0pbOuOnOuCcC71tpaQH1gq4/yFYS8zgzQj8Ce9Sd5mfk8kGytrQM0Bh4LpP+f89hH7YDquUvd6Q3IAAAC8ElEQVQSMMXJDH5f7Nbardba7VdZLQHItNbutNaeBdKBzgWfrsB0BmbkfjwDuPMy6xUCoo0xhYAiwDc+yFZQrjqzMaY40Ax4FcBae9Za+9u3tA8cefo+G2MqAh2AV3yUqyBddWZr7UFr7We5H58g5x+0Cj5L6L289FFnYGbuGyN9CpQwxpRzKoDfF3seVQD2/ezz/QTWD8Kveay1B3M/PgR4fr2CtfYAkAbsBQ4CP1prl/suouOuOjNQFfgWeC33tMQrxpgYnyV0Xl5mBhhPzjuVZfskVcHK68wAGGOqAA2ANQUby1F56aMC7Syv3mjDKXl5w+xgc6WZf/6JtdYaY35z6VLuucnO5JTdMeBfxpge1to3CiKvE7ydmZyf198Dfay1a4wxE8j5VX6E42Ed4sD3uSNwxFq7wRjTomBSOsuB7/NP+4kF5gFPWGuPO5syuPlFsVtrb/NyFweA63/2ecXcv/NbV5rZGHPYGFPOWnsw99ezI5dY7TZgl7X229xt5gO3An5b7A7MvB/Yb6396ehtLlc+L+06B2ZuCtxhjGkPRAHFjDFvWGt7FFBkrzkwM8aYCHJKfZa1dn4BRS0oeemjAu2sYDkVsw6oboypaowpTM6TiO+4nMkb7wB/zf34r8ClfmvZCzQ2xhQxxhigFYH95NpVZ7bWHgL2GWNq5v5VK2CLb+IViLzMPNRaW9FaW4Wcn+sP/LnU8+CqM+f+PL8KbLXWjvNhNqfkpY/eARJzr45pTM6p1IO/3lG+WWv9egG6kHOkdgY4DCzL/fvywJKfrdce+BrYQc4pHNezezFzKXKuGPgP8D4Qd5mZnwK2kXPl0OtApNvZfTDzzcB6YCPwFlDS7ewFPfPP1m8BLHI7d0HPDPw/wOZ+j7/IXdq7nf0a5/xNHwG9gd65HxtyrpzZAWwi5+o2xx5frzwVEQkywXIqRkREcqnYRUSCjIpdRCTIqNhFRIKMil1EJMio2EVEgoyKXUQkyKjYRUSCzP8Hl/ZXsrspAhQAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"def func(x):\n",
" return -4*np.exp(-x)+2*np.exp(-2*x)\n",
"\n",
"x1, x2 = -1.0, 0.0\n",
"x = np.linspace(x1,x2, 101)\n",
"y = func(x)\n",
"plt.plot(x, y, color = 'k')\n",
"plt.plot([x1,x2],[0,0])\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" x1 x2 f1 f2\n",
" -1.0000000000 +0.0000000000 +3.9049848840 -2.0000000000\n",
" -1.0000000000 -0.5000000000 +3.9049848840 -1.1583214259\n",
" -0.7500000000 -0.5000000000 +0.4953780742 -1.1583214259\n",
" -0.7500000000 -0.6250000000 +0.4953780742 -0.4922979148\n",
" -0.7500000000 -0.6875000000 +0.4953780742 -0.0447964325\n",
" -0.7187500000 -0.6875000000 +0.2128474183 -0.0447964325\n",
" -0.7031250000 -0.6875000000 +0.0810265592 -0.0447964325\n",
" -0.6953125000 -0.6875000000 +0.0173789137 -0.0447964325\n",
" -0.6953125000 -0.6914062500 +0.0173789137 -0.0138911236\n",
" -0.6933593750 -0.6914062500 +0.0016980959 -0.0138911236\n",
" -0.6933593750 -0.6923828125 +0.0016980959 -0.0061079375\n",
" -0.6933593750 -0.6928710938 +0.0016980959 -0.0022077800\n",
" -0.6933593750 -0.6931152344 +0.0016980959 -0.0002555572\n",
" -0.6932373047 -0.6931152344 +0.0007210905 -0.0002555572\n",
" -0.6931762695 -0.6931152344 +0.0002327219 -0.0002555572\n",
" -0.6931762695 -0.6931457520 +0.0002327219 -0.0000114288\n",
" -0.6931610107 -0.6931457520 +0.0001106438 -0.0000114288\n",
" -0.6931533813 -0.6931457520 +0.0000496068 -0.0000114288\n",
" -0.6931495667 -0.6931457520 +0.0000190888 -0.0000114288\n",
" -0.6931476593 -0.6931457520 +0.0000038299 -0.0000114288\n",
" -0.6931476593 -0.6931467056 +0.0000038299 -0.0000037995\n",
"\n"
]
}
],
"source": [
"x1, x2 = -1.0, 0.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",
"x0 = -np.log(2.0)\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": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGXa+PHvPSkEkhBKSAgTIEBCCYQakWJBQQUBAUWFtaAi2MXd9d3Vn7u+23x33XVdy9pAEFQEwYKoiK4VpBp6hxhaAiShJpT05/dHhjWGBCaZck4m9+e65iJzzpnz3DkZ5p6nnOcRYwxKKaXqH4fVASillLKGJgCllKqnNAEopVQ9pQlAKaXqKU0ASilVT2kCUEqpekoTgFJK1VOaAJRSqp7SBKCUUvVUsNUBVEVERgIjIyMjJ3Xs2NHqcJRSqk5Zs2bNYWNMiwsdJ3aeCiI1NdWkpaVZHYZSStUpIrLGGJN6oeO0CUgppeopWyYAERkpIlNPnDhhdShKKRWwbJkAjDEfG2MmR0VFWR2KUkoFLFsmAKWUUr5nywSgTUBKKeV7tkwA2gSklFK+Z8sEoJRSyvcCMgEs3nyIaUsyrA5DKaVsLSATwJfbspm5fI/VYSillK3ZMgFoJ7BSSvmeLROAdgIrpZTv2TIBKKWU8j1NAEopVU9pAlBKqXoqYBNAcWkZdp7qWimlrGbLBODpKKAUZxQ5+YX868tdXo5MKaUChy0TgKejgG7v35axfeJ54atdzF29z8vRKaVUYLDlkpCeEhH+en0KOfmFPLFgM7GNw7iic4zVYSmllK3YsgbgDSFBDl6+pTdd4iK5f/ZaNuw/bnVISillKwGbAAAiGgQz446LaB4Ryl0zf2DvkVNWh6SUUrYR0AkAICYyjFl39aXUGCbMWM2Rk4VWh6SUUrbgtwQgIu1FZLqIvOevMs/q0CKC6RNSOXiigImz0jhTVOrvEJRSynbcSgAiMkNEckRkc6XtQ0Vkh4iki8hj5zuHMSbDGDPRk2A90adtM54f14sNmcd5aM5aSkrLrApFKaVswd0awExgaMUNIhIEvAQMA5KB8SKSLCIpIvJJpYcthuAM7daSP17XlS+35fDkwi16o5hSql5zaxioMWaJiCRU2twXSDfGZACIyFxglDHmr8AIbwbpTbf3T+DA8QJe/e5HWkWF8eCVSVaHpJRSlvCkD8AJ7K/wPNO1rUoi0lxEXgV6icjj5zlusoikiUhabm6uB+FV7zfXdGJ0z1Y888VO3luT6ZMylFLK7vx2I5gx5ghwrxvHTRWRg8DI0NDQPr6IxeEQ/j62B7knC3ns/Y3ERDbgso4tfFGUUkrZlic1gCygdYXn8a5tHvPHgjChwQ5eubUPiTER3Pf2GjZn6epjSqn6xZME8AOQJCLtRCQUGAcs9EZQ/loSsnFYCLPu6ktUwxDunPkD+4+e9ml5SillJ+4OA50DrAA6iUimiEw0xpQADwKfA9uAecaYLd4Iyp9LQsY2Lr9RrLC4lAlvrObYqSKPzldYUsq2g3l617FSyvbEjkMhRWQkMDIxMXHSrl3+mdJ5VcYRbpu+mpT4KGbffTFhIUHnPb6wpJTdh0+xM/sku7Lz2ZV9kp05+ew9cprSMkNkWDDLH7uSyLAQv8SvlFJnicgaY0zqBY+zYwI4KzU11aSlpfmtvE83HuTBOWu5OjmWl2/pQ5BDKCopc33Q57MrO7/8Az8nnz2uD3oAh0BC83ASYyLoGBtJVMMQnlq0jf8dmcydA9v5LX6llAL3E4Atp4OuUAPwa7nDu8dxKC+ZP3+ylRteWU5+QfE5H/Rtm4eTFBPB0G4t6RgbSVJMJO1bhJ9TY1i85RAzl+/h9v4JBDnEr7+HUkq5w5YJwBjzMfBxamrqJH+XPfGSdpwsKOHjjQdo36L8gz4pJpKk2Ag6tIi4YNPQWXcOTODBd9bxzfYchiTH+jhqpZSqOVsmAKtNGZLElCGe3SF8TdeWxEWFMWPZbk0ASilbsuV00P4aBupLIUEObu+fwPIfj7D9UJ7V4Sil1DlsmQD8OQzUl8b3bU1YiIM3vt9jdShKKXUOWyaAQNGkUSjX945nwfosjnp4f4FSSnmbLRNAIDQBnXXngAQKS8qYs3qf1aEopdTP2DIBBEoTEEBSbCSXJkXz5oo9FOsiNEopG7FlAgg0dw1sR3ZeIYs2HbQ6FKWU+i9NAH5weccWtI8OZ8ayPVaHopRS/6UJwA8cDuGOgQls2H+ctfuOWR2OUkoBNk0AgdQJfNYNveOJDAtmxve7rQ5FKaUAmyaAQOoEPiu8QTDjLmrNZ5sPcfDEGZ+W9fGGAzz4zlodeqqUOi9bJoBAdXv/BIwxvLlir8/KSM/J53/e28AnGw9ywyvL2XNY1yVQSlVNE4AftW7WiKuTWzJn9T7OFJV6/fyFJaU8NGc9jUKDefXW3hw/XcT1ryxnzV7td1BKnUsTgJ/dOTCB46eL+XCdV5ZP/pm/L97BtoN5/GNsd4Z2i+PD+wfSOCyY8dNW6hBUpdQ5NAH4Wd92zejaqjEzl+/Gm4vxfLsjh+nf72ZC/7YM7lI++2hCdDgf3D+QFGcUD7yzlmlLMrxaplKqbrNlAgjEUUBniQh3DmzHzuyTLEs/4pVzHj5ZyKPzN9IpNpLHr+3ys33NwkOZfffFXNstjqcWbePJj7ZQonckK6WwaQIIxFFAFY3sEUd0RCgzlnk+JNQYw//M30BeQTHPj+9Z5YI1YSFBvDi+F/dc1p63Vu7lnrfWcKqwxOOylVJ1my0TQKBrEBzELRe35evtOez2cJTOrOV7+GZHLk9c24XOLRtXe5zDITx+bRf+PLob3+zI4eapK8jJK/CobKVU3aYJwCK39GtDSJAw04NawPZDefzfZ9u5snMMt/dv69ZrbuvXlukTLiIj9xRjXl7Ozuz8WpevlKrbNAFYJCYyjJE9WjF/TSYnzhTX+PUFxaU8PGcdjcNC+PvY7oi4v/D8FZ1jmHdPf4pLy7jhleUsTz9c4/KVUnWfJgAL3TWwHaeLSpmftr/Gr/3rom3szD7JP2/qQXREgxq/vpszig8fGEhcVBi3z1jNe2sya3wOpVTd5tcEICKjRWSaiLwrIlf7s2w76uaMom9CM2Yu30NpmfvDM7/als2sFXuZeEk7Lu/YotblO5s05L37BnBx+2Y8On8Dz325U4eJKlWPuJ0ARGSGiOSIyOZK24eKyA4RSReRx853DmPMAmPMJOBe4ObahRxY7hyYQOaxM/xna7Zbx+fkFfA/722kS1xjfjO0k8flNw4L4Y07+jK2TzzPfbmLR+dvpKhEh4kqVR/UpAYwExhacYOIBAEvAcOAZGC8iCSLSIqIfFLpEVPhpb9zva7euyo5FmeThrzhRmdwWZnh1/M3cLqohBfH96RB8LlDPmsjNNjBP8Z251dXdeT9tZnc8cbqWvVLKKXqFrcTgDFmCXC00ua+QLoxJsMYUwTMBUYZYzYZY0ZUeuRIuaeBz4wxa733a9RdwUEOJgxoy6rdR9ly4Pw3vs1Ytpuluw7z+xHJJMZEejUOEeHhwUk8e1MPfthzlJtfW8HpIr1XQKlA5mkfgBOo2IOZ6dpWnYeAIcBYEbm3qgNEZLKIpIlIWm5urofh1Q03p7ahYUgQb5xnxbDNWSd4evF2rk6O5Rd92/gslut7xzP19lR2ZOfzl0+3+awcpZT1/NoJbIx5wRjTxxhzrzHm1WqOmQr8EVgbGhrqz/AsE9UohLF94lm4/gCHTxaes/90UQlT5q6jWXgoT99QsyGftXFFpxgmX9qed1bt44sth3xallLKOp4mgCygdYXn8a5tHgn0qSCqcsfABIpKy5i9ct85+/78yTYyDp/i2Zt60jTcP0nxV1d3pGurxjz2wSa9Y1ipAOVpAvgBSBKRdiISCowDFnoaVCBPBledDi0iGNSpBW+v2kthyU9rBSzefIg5q/cx+bL2DEyM9ls8DYKDeH5cT04XlfDr+Rsoq8EwVaVU3VCTYaBzgBVAJxHJFJGJxpgS4EHgc2AbMM8Ys8XToOpjDQDgzoHtyM0v5NON5XP3Hzxxhsc+2EiKM4pfX+X5kM+aSoyJ5HfDk1m66zAzl+/xe/lKKd8KdvdAY8z4arYvAhZ5LSLKawDAyMTERG+e1vYuS4omMSaCGct2M6qnk1+9u4GikjKeH9eT0GBrbtq+5eI2fLsjh78t3s6AxObnnXBOKVW32HIqiPpaAxAR7hiQwOasPB6as5YVGUf4w8iutG8RYWlMf7uhO43DQpgyZz0Fxd5fylIpZQ1bJoD62Adw1vW9nUQ1DGHRpkMMT4njxtR4q0MiOqIBz9zYnR3Z+Ty9eLvV4SilvMSWCaC+1gAAGoUG88AVHejcMpL/G5Pi8yGf7hrUKYY7BiTwxrI9fLezftyfoVSgEztP/pWammrS0tKsDsMSxhjbfPifVVBcynX//p5jp4tZPOVSmtdiFlKllO+JyBpjTOqFjrNlDaA+NwGdZbcPfyhfWvL5cb04cbqY376/SWcOVaqOs2UCqM9NQHZ3dhbSL7dlM2d1zdcxUErZhy0TgLK3uwa249KkaP70yRbSc05aHY5SqpZsmQC0CcjeHA7hmRt70DAkiEfeXafrByhVR9kyAWgTkP3FNg7jbzd0Z3NWHv/6cqfV4SilasGWCUDVDdd0bcn4vq159bsfWfHjEavDUUrVkCYA5ZHfj0imXfNwfjVvPSdO6ypiStUlmgCURxqFBvPcuJ7k5hfy/xbo0FCl6hJbJgDtBK5busc34ZdXdeTTjQf5YK3Hy0EopfzElglAO4Hrnnsv70Dfds148qPN7Dty2upwlFJusGUCUHVPkEP41809cTiER95dR0mpDg1Vyu40ASivcTZpyFNjUli77zh//3wHR08VWR2SUuo83F4QRil3XNejFd/uyGHqkgymLsnA2aQhKc4oUuKjyv91RvltXWOl1PlpAlBe94+xPRjbJ57NWSfYmHmCzVknWLzl0H/3xzdtSPf4KLo5o+jubEKKM4qoRiEWRqxU/WTLBFBfl4QMFEEOYUCHaAZ0+GkR+xNnitmSdYKNWSfYlHWCTZknWLTpp6TQplmjn9UUusdHERmmSUEpX9L1AJRljp8uYnNWXnlCyDrOpqwT7D96BihfheyLX15GM20uUqrG3F0PwJY1AFU/NGkUyiVJ0VyS9FNN4dipIlbtPsL9s9fy2nc/8vi1XSyMUKnApqOAlK00DQ9laLc4RvV0MmvFHnLyC6wOSamApQlA2dKUwUkUlxpe/uZHq0NRKmBpAlC2lBAdzo194nln1T4OHD9jdThKBSS/JQAR6SIir4rIeyJyn7/KVXXXQ4OTAHjx63SLI1EqMLmVAERkhojkiMjmStuHisgOEUkXkcfOdw5jzDZjzL3ATcDA2oes6gtnk4aM79ua+Wn72XvklNXhKBVw3K0BzASGVtwgIkHAS8AwIBkYLyLJIpIiIp9UesS4XnMd8CmwyGu/gQpoD1yRSJBDeP6rXVaHolTAcSsBGGOWAEcrbe4LpBtjMowxRcBcYJQxZpMxZkSlR47rPAuNMcOAW7z5S6jAFdM4jAkDEliwLov0nHyrw1EqoHjSB+AE9ld4nunaViURGSQiL4jIa5ynBiAik0UkTUTScnNzPQhPBYp7LmtPWEgQ//pSawFKeZPfbgQzxnwLfOvGcVNF5CAwMjQ0tI+v41L21zyiAXcNbMe/v0nnwSvy6BLX2OqQlAoIntQAsoDWFZ7Hu7Z5TBeEUZVNurQ9kWHBPPufnVaHolTA8CQB/AAkiUg7EQkFxgELvRGULgmpKotqFMLkS9vzn63ZbNh/3KdlzV29jxEvLmVXtvY5qMDm7jDQOcAKoJOIZIrIRGNMCfAg8DmwDZhnjNnijaC0BqCqcucl7WjaKIR/+rAW8N3OXJ5YsJnNWXnc9NoKNmXqlxAVuNwdBTTeGBNnjAkxxsQbY6a7ti8yxnQ0xnQwxjzlraC0BqCqEtEgmPsGdWDJzlxW7648KM1zO7PzeXD2WpJiIvj04UtoFBrM+GkrWZVxxOtlKWUHtpwKQmsAqjq39UugRWQDnvliB96cyvzwyULumvkDYaFBzLjjIrq2iuK9+/oT27gBt89YzTfbc7xWllJ2YcsEoDUAVZ2GoUE8MKgDq3cfZVm6d76ZFxSXMvnNNHLzC5l2eyqtmjQEIC6qIfPu6U9SbAST3kzjk40HvFKeUnZhywSgNQB1PuMvbkOrqDCv1AKMMfz2/Y2s3XecZ2/qSc/WTX62v3lEA96Z1I9ebZrw0Jx1zF29z6PylLITWyYApc6nQXAQDw1OYv3+43yzw7OmmRe+Suej9Qd49OqODO8eV+UxjcNCePOui7ksqQWPfbCJ15dmeFSmUnZhywSgTUDqQsb2iadNs0b884udlJXVrhawcMMB/vXlTq7v5eSBK86//nTD0CCm3Z7K8JQ4/vLpNp71ch+EUlawZQLQJiB1ISFBDh4ZksSWA3l8vuXQhV9Qydp9x3h0/gYuSmjKX29IQUQu+JrQYAcvjO/FTanxvPB1On/8eGutk49SdmDLBKCUO0b1dNKhRTjP/mcnpTX4IM48dprJb6bRsnEYr92WSoPgILdfG+QQnr6hOxMvacfM5Xv4zfsbKSktq034SlnOlglAm4CUO4Icwi+v6siunJN8vMG9ETr5BcVMnJlGYUkZM+5IpVl4aI3LFRF+N7wLvxzSkffWZPLgO+soLCmt8XmUspotE4A2ASl3Xdstjs4tI3nuy50UX+CbeElpGQ/NWUd67kleuaUPiTGRtS5XRJgyJIknRySzeMsh7p6VxumiklqfTykr2DIBKOUuh0P49dWd2HPkNB+szTzvsX/5dBvf7sjlj9d15ZKkaK+Uf9cl7fj72O4sSz/MbdNXc+JMsVfOq5Q/aAJQdd6QLjH0iI/iha/Sq22KeWvFHmYu38NdA9txa7+2Xi3/ptTWvPSL3mzMPM74qSs5fLLQq+dXylc0Aag6T6S8FpB1/Azzfth/zv7vdubyh4+3cmXnGJ4Y3sUnMQxLieP1CReRcfgkN726ggPHz/ikHKW8yZYJQDuBVU1dmhRN34RmvPh1OgXFP9UCKk7w9sL4XgQ5Ljzcs7Yu79iCtydeTG5+ITe+uoKThdonoOzNlglAO4FVTZXXAjqSk1/I2yv3AnDENcFbg5Agpt9xERENfL8AXmpCM2bceRFZx88wa/ken5enlCdsmQCUqo2L2zfnksRoXv72R46eKmLyW2vIzS/k9QmpOF0TvPnDRQnNuLJzDNOWZpBfoJ3Cyr40AaiA8qurO3L0VBHDX1jKmr3H+OdNPc6Z4M0fHhmSxPHTxVoLULamCUAFlN5tmjK4cwwHTxTw66s6MqJ7K0vi6B7fhCFdYpi2dDd5WgtQNqUJQAWcv96QwnM39+TBK88/wZuvTRnckRNnipm1bI+lcShVHVsmAB0FpDwRExnG6F5OtyZ486WU+CiGdIll2tIMrQUoW7JlAtBRQCpQPDIkibyCEt74fo/VoSh1DlsmAKUCRTdnFFcnx/L69xk6TYSyHU0ASvnYlCFJ5BeU8May3VaHotTPaAJQyse6torimq6xTP9+t9YClK1oAlDKD6YM7kh+QQnTv9dagLIPvyYAEQkXkTQRGeHPcpWyWnKrxgzt2pI3vt/NidNaC1D24FYCEJEZIpIjIpsrbR8qIjtEJF1EHnPjVL8F5tUmUKXquilDksgvLGH69xlWh6IU4H4NYCYwtOIGEQkCXgKGAcnAeBFJFpEUEfmk0iNGRK4CtgI5XoxfqTqjS1xjrk1pyYxlezh+usjqcJRyLwEYY5YARytt7gukG2MyjDFFwFxglDFmkzFmRKVHDjAI6Af8ApgkItr/oOqdhwcncbKwhNeXal+Asp4nH8JOoOLqG5mubVUyxjxhjHkEeAeYZoypcgFXEZns6idIy83N9SA8peync8vGDE+JY+byPRw75ZtawJmiUtbsrfx9Talz+f1buDFmpjHmk/Psnwr8EVgbGhrqv8CU8pOHBydxqqiE133QF5BXUMyt01dxwysr+PfXu7x+fhVYPEkAWUDrCs/jXds8plNBqEDWqWUk16bEMXPZHo56sRZw7FQRt0xbxcbM4/Rv35xnvtjJ60u1w1lVz5ME8AOQJCLtRCQUGAcs9EZQOhmcCnSPDE7idHEp07z0AZ2TX8C4qSvZkZ3P1NtSeWtiX4anxPGXT7f9d4U0pSpzdxjoHGAF0ElEMkVkojGmBHgQ+BzYBswzxmzxRlBaA1CBLik2khHdWzFruee1gAPHz3DzayvZf+w0M++4iCs6xxAc5OBfN/dkcOcYfrdgM++vyfRS5CqQuDsKaLwxJs4YE2KMiTfGTHdtX2SM6WiM6WCMecpbQWkNQNUHD1+ZyJniUqYuqX0tYO+RU9z46goO5xfy1sS+DEiM/u++0GAHL93Sm0sSo/mf9zbw6caD3ghbBRBbDsXUGoCqD5JiIxnZvRVvrtjDkZOFNX59ek4+N722glNFJbwzqR992jY755iwkCCm3t6HPm2bMmXuOr7cmu2FyFWgsGUC0BqAqi8eHpxUq1rA1gN53PzaSkrL4N3J/UmJr/7LUqPQYGbccRFdWzXm/tlrWbpLh1ercrZMAFoDUPVFYkwE1/VoxZsr9nLYzVrAun3HGDd1BaHBDubd049OLSMv+JrIsBBm3dWX9i3CmfRmGqt3630CyqYJQKn65OHBSRSWuFcLWJVxhFtfX0WTRqHMu6c/7VtEuF1Ok0ahvH33xTibNOSumT+wfv9xT8JWAcCWCUCbgFR90qFFBKN6OnlzxR5y86uvBXy3M5cJb6wmrklD5t/bn9bNGtW4rOiIBsy+ux/NwkO5ffoqth7I8yByVdfZMgFoE5Cqbx66MpGikjKmLvmxyv1fbDnEpFlptIuO4N3J/YhtHFbrslpGhTH77ouJaBDMbdNXkZ6TX+tzqbrNlglAqfqmfYsIRvd08tbKveTkF/xs30frs7hv9lqSWzVm7qR+NI9o4HF5rZs1Yvakfjgcwi2vr2LvkVMen1PVPbZMANoEpOqjhwYnUVRSxmvf/dQXMO+H/Tzy7nr6tG3K23dfTFSjEK+V1y46nNl3X0xRSRm/mLaKrONnvHZuVTfYMgFoE5Cqj9pFhzO6l5O3V+4lJ6+AWcv38Jv3N3JJYjSz7uxLRINgr5fZMTaStyZeTF5BMbdMW0lOXsGFX6QChi0TgFL11cNXJlFSZrh9xmr+d+EWrkqO5fUJqTQMDfJZmd2cUcy8sy85+YXc8voqr05Qp+xNE4BSNpIQHc6YXk62H8pnZI9WvHxLbxoE++7D/6w+bZsyfcJF7Dt6mtumr+LEGV23uD7QBKCUzfx+eDLP3dyT527uSUiQ//6L9u/QnNdu68PO7HzueGM1RSVVrtmkAogtE4B2Aqv6LKpRCKN7OQlyiN/LHtQphqdv6M66fcf5docu3x3obJkAtBNYKetc16MV0RGhLFjvlfWdlI3ZMgEopawTHORgZI9WfLktR/sCApwmAKXUOcb0clJUUsbizbqGQCDTBKCUOkeKM4r2LcL5YK02AwUyTQBKqXOICGN6Olm1+6jeIRzAbJkAdBSQUtYb1dMJlM9FpAKTLROAjgJSynptmjcitW1TPlybhTHG6nCUD9gyASil7GF0Lye7ck6y9aCuGxCINAEopao1PCWOkCBhwTptBgpEmgCUUtVqGh7KFZ1i+Gj9AUrL/NMMlJNXwP8t2sbBE9r57GuaAJRS5zWml5Oc/EJW/HjEL+X984udTF2SwbXPL+Xr7dl+KbO+8lsCEJFBIrJURF4VkUH+Klcp5ZkrOscQGRbMh35oBjpw/AwfrMtkWLeWtIxqyF0z0/jLJ1t1YjofcSsBiMgMEckRkc2Vtg8VkR0iki4ij13gNAY4CYQBmbULVynlb2EhQQxPiWPx5oOcKSr1aVnTlmZgDDwxvAsf3j+A2/q15fXvd3Pjq8vZd+S0T8uuj9ytAcwEhlbcICJBwEvAMCAZGC8iySKSIiKfVHrEAEuNMcOA3wJ/9N6voJTytdG9nJwqKuWLrYd8VsaRk4XMWb2P0b2cxDdtRFhIEH8e3Y1XbulNxuFTDH9hKZ9u1KkpvMmtBGCMWQIcrbS5L5BujMkwxhQBc4FRxphNxpgRlR45xpizdbhjgOerWiul/KZvQjNaRYX5dDTQjGW7KSwp497LO/xs+7CUOBY9fCkdYiJ44J21PPHhJgqKfVsTqS886QNwAvsrPM90bauSiFwvIq8BbwH/Ps9xk0UkTUTScnNzPQhPKeUtDocwqpeTJbsOc/hkodfPn1dQzJvL9zKsW0sSYyLO2d+6WSPm39ufey5vz+xV+xj90jLSc056PY76xm+dwMaYD4wx9xhjbjbGfHue46ZS3kS0NjQ01F/hKaUuYEwvJ6Vlhk82HPD6ud9asZf8whLuH5RY7TEhQQ4eH9aFN+68iJz8Qka++D3vrdHuRE94kgCygNYVnse7tnlMp4JQyn46xkaSHNeYD9d7NwGcKSplxve7GdSpBd2cF/4/f0WnGD6bcik9Wkfx6PwN/Greek4Vlng1pvrCkwTwA5AkIu1EJBQYByz0RlA6GZxS9jSml5MN+4+Tkeu95pe5P+zjyKkiHrii+m//lcU2DmP23f14ZEgSC9ZlMfLf37P1gE5XUVPuDgOdA6wAOolIpohMNMaUAA8CnwPbgHnGmC3eCEprAErZ03U9W+EQWOClWkBRSRlTl2TQN6EZFyU0q9FrgxzCI0M6MvvufpwsKGH0y8t4a+VenbiuBtwdBTTeGBNnjAkxxsQbY6a7ti8yxnQ0xnQwxjzlraC0BqCUPcU2DmNgYjQL1nlnhtAF67I4eKKAB650/9t/Zf07NOezKZcyoENzfr9gM/fPXqtLWbrJllNBaA1AKfsa3dPJvqOnWbvvuEfnKS0zvPLdj3RzNuaypGiPztU8ogEzJlzE48M685+t2Qx/YSk7DuV7dM76wJYJQGsAStnXNd1aEhbi8PiegEWbDrL78CkeGJSIiHgcl8Mh3HN5B+bd25+ikjLGvrqc5T8e9vi8gcyWCUBrAErZV0S7j72dAAAMeUlEQVSDYK5ObsknGw/Ueo4eYwwvfZNOhxbhXNO1pVfj692mKR8+MJCWjcOYMGO1rmh2HrZMAEopexvTy8mx08V8t7N2N2t+syOH7YfyuX9QIg6H59/+K3M2ach79w6gd5umTJm7nle+/VE7h6tgywSgTUBK2dslSdE0Dw+tVTOQMYZ/f52Os0lDruvZygfRlYtqFMKbE/syskcrnl68nd9/tNlvaxrUFbZMANoEpJS9hQQ5GNmjFf/Zlk1eQc1G3KzMOMrafce59/L2hAT59iOoQXAQz9/ck3sub8/bK/dxz1trfD6jaV1iywSglLK/0b2cFJWUsXhTzWYIffnbdKIjGnBjausLH+wFDofw+LAu/GlUV77ans34aSs54oP5jOoiWyYAbQJSyv56xEfRLjq8RgvFbNh/nKW7DnP3pe0ICwnyYXTnur1/Aq/e2odtB/O4/pXl7Dl8yq/l25EtE4A2ASllfyLC6J5OVu4+woHj7q3f+9I36TQOC+bWfm19HF3VrunakjmT+5FfUML1ryxn7b5jlsRhF7ZMAEqpumFMLyfGwEI3ZgjdmZ3PF1uzuWNgOyIaBPshuqr1btOU9+8bQGRYML+YtpIvtvhukRu70wSglKq1Ns0b0adtU7dGA73y7Y80Cg3izgEJvg/sAtpFh/P+fQPo1LIx9769hjdX7LE6JEtoAlBKeWR0LyfbD+Wz7WD1s3HuO3KahRsO8Iu+bWgabo91PqIjGjB3Uj+u7BzDkx9t4a+fbaOsng0TtWUC0E5gpeqOESlxBDvkvLWAV5f8SJAIky5r78fILqxhaBCv3ZbKrf3a8Np3GTzy7noKS+rPMFFbJgDtBFaq7mgaHsqgTjEsWJ9V5Y1W2XkFvJeWydjUeGIbh1kQ4fkFOYQ/j+rGb4d2ZuGGA0yYsbrezCZqywSglKpbxvRykp1XyMqMI+fsm7Ykg5KyMu69rEMVr7QHEeG+QR14flxP1uw9xthXlnPsVJHVYfmcJgCllMcGd4khskHwOfcEHDtVxOxV+7iuRyvaNG9kUXTuG9XTyaw7+7L78CmeXOiV9a1sTROAUspjYSFBDEtpyeLNh3421cIby/dwpriU+2uw3KPVBiRG8/DgJD7ecIBFmw5aHY5PaQJQSnnF6F5OThaW8OW2bABOFpYwc9lurk6OpWNspMXR1cx9gzqQ4ozidws2cziAp42wZQLQUUBK1T392jUnLirsv6OB3l65l7yCkjr17f+skCAH/7ypBycLSvjdh5sDdippWyYAHQWkVN3jcAjX9WzFdztzOXD8DK8v3c0lidH0bN3E6tBqpWNsJL+6uiOLtxxy607nusiWCUApVTdd3yuekjLDPW+t4fDJQh6og9/+K5p0aXt6tWnCkx9tITuvwOpwvE4TgFLKazq1jKRLXGM2ZZ2gd5sm9GvfzOqQPBLkEJ65sQcFxaU8/sGmgGsK0gSglPKqMb3KV/m630uLvVutQ4sIfjO0M19vz+G9NZlWh+NV1k3Jp5QKSBMGJJAUG8mgji2sDsVr7hyQwOdbDvGnj7cyMDGaVk0aWh2SV/itBiAiDhF5SkReFJEJ/ipXKeVfDYKDuKJTTEB8+z/L4RCeGduDUmP47fsbA6YpyK0EICIzRCRHRDZX2j5URHaISLqIPHaB04wC4oFiILDqUUqpgNemeSMev7YLS3cd5p3V+6wOxyvcrQHMBIZW3CAiQcBLwDAgGRgvIskikiIin1R6xACdgOXGmF8B93nvV1BKKf+4pW8bBiY256lPt7H/6Gmrw/GYWwnAGLMEOFppc18g3RiTYYwpAuYCo4wxm4wxIyo9cij/1n92/bUyb/0CSinlLw6H8PexPXCI8Oj8DXV+/QBP+gCcwP4KzzNd26rzAXCNiLwIfFfdQSIyWUTSRCQtNzfXg/CUUsr7nE0a8vsRXVi1+2idX0nMb6OAjDGngYluHDdVRA4CI0NDQ/v4PjKllKqZm1Jbs3jzIf62eDuXd4qhXXS41SHViic1gCygdYXn8a5tHtOpIJRSdiYi/PX67oQGOXh0/oYqF8KpCzxJAD8ASSLSTkRCgXHAQm8EpZPBKaXsrmVUGH+4ritr9h5j+vcZVodTK+4OA50DrAA6iUimiEw0xpQADwKfA9uAecYYr6ygoDUApVRdMKaXk6uSY3nmi53sys63OpwaEzve0CAiI4GRiYmJk3bt2mV1OEopVa3c/EKu/td3tGnWiPfvG0BwkPUz7IjIGmNM6oWOsz7SKmgNQClVV7SIbMCfR3djQ+YJXltSt5qCbJkAtA9AKVWXjOjeiuHd43juy51sO5hndThus2UC0BqAUqqu+fOobkQ1DOHX8zZQVFI37nW1ZQJQSqm6pll4KE+NSWHrwTz+/U261eG4xZbTQVfoBLY6FKWUcts1XVsyppeTF77axYzvdyMCQvkUEg4RHFJ+D4FDcD2X8mMqPXeI8OxNPege79vlNG2ZAIwxHwMfp6amTrI6FqWUqok/jepKQvNw8gqKMQbKjMEYQ5nr5zKD6/lP28x/t/30vFFokM9jtWUCUEqpuioyLIQpQ5KsDsMttuwD0FFASinle7ZMADoKSCmlfM+WCUAppZTvaQJQSql6ShOAUkrVU7ZMANoJrJRSvmfLBKCdwEop5Xu2TABKKaV8z5brAZwlIrnA3lq+PBo47MVwvE3j84zG5xmNzzN2j6+tMabFhQ6ydQLwhIikubMgglU0Ps9ofJ7R+Dxj9/jcpU1ASilVT2kCUEqpeiqQE8BUqwO4AI3PMxqfZzQ+z9g9PrcEbB+AUkqp8wvkGoBSSqnzqPMJQESGisgOEUkXkceq2C8i8oJr/0YR6e3H2FqLyDcislVEtojIlCqOGSQiJ0RkvevxpL/ic5W/R0Q2ucpOq2K/ldevU4Xrsl5E8kTkkUrH+PX6icgMEckRkc0VtjUTkf+IyC7Xv02ree1536s+jO8fIrLd9ff7UESqXGbqQu8FH8b3BxHJqvA3vLaa11p1/d6tENseEVlfzWt9fv28zrhWq6mLDyAI+BFoD4QCG4DkSsdcC3xG+cps/YBVfowvDujt+jkS2FlFfIOATyy8hnuA6PPst+z6VfG3PkT5+GbLrh9wGdAb2Fxh29+Bx1w/PwY8XU38532v+jC+q4Fg189PVxWfO+8FH8b3B+BRN/7+lly/Svv/CTxp1fXz9qOu1wD6AunGmAxjTBEwFxhV6ZhRwJum3EqgiYjE+SM4Y8xBY8xa18/5wDbA6Y+yvciy61fJYOBHY0xtbwz0CmPMEuBopc2jgFmun2cBo6t4qTvvVZ/EZ4z5whhT4nq6Eoj3drnuqub6ucOy63eWiAhwEzDH2+Vapa4nACewv8LzTM79gHXnGJ8TkQSgF7Cqit0DXNXzz0Skq18DAwN8KSJrRGRyFfttcf2AcVT/H8/K6wcQa4w56Pr5EBBbxTF2uY53UV6jq8qF3gu+9JDrbzijmiY0O1y/S4FsY8yuavZbef1qpa4ngDpBRCKA94FHjDF5lXavBdoYY7oDLwIL/BzeJcaYnsAw4AERuczP5V+QiIQC1wHzq9ht9fX7GVPeFmDLoXUi8gRQAsyu5hCr3guvUN600xM4SHkzix2N5/zf/m3/f6myup4AsoDWFZ7Hu7bV9BifEZEQyj/8ZxtjPqi83xiTZ4w56fp5ERAiItH+is8Yk+X6Nwf4kPKqdkWWXj+XYcBaY0x25R1WXz+X7LPNYq5/c6o4xur34R3ACOAWV5I6hxvvBZ8wxmQbY0qNMWXAtGrKtfr6BQPXA+9Wd4xV188TdT0B/AAkiUg717fEccDCSscsBG53jWbpB5yoUF33KVeb4XRgmzHm2WqOaek6DhHpS/nf5Iif4gsXkcizP1PeWbi50mGWXb8Kqv3mZeX1q2AhMMH18wTgoyqOcee96hMiMhT4DXCdMeZ0Nce4817wVXwV+5TGVFOuZdfPZQiw3RiTWdVOK6+fR6zuhfb0QfkolZ2UjxB4wrXtXuBe188CvOTavwlI9WNsl1DeHLARWO96XFspvgeBLZSPalgJDPBjfO1d5W5wxWCr6+cqP5zyD/SoCtssu36UJ6KDQDHl7dATgebAV8Au4EugmevYVsCi871X/RRfOuXt52ffg69Wjq+694Kf4nvL9d7aSPmHepydrp9r+8yz77kKx/r9+nn7oXcCK6VUPVXXm4CUUkrVkiYApZSqpzQBKKVUPaUJQCml6ilNAEopVU9pAlBKqXpKE4BSStVTmgCUUqqe+v/YXqNtW3TMMwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"X = list_bisec[0]\n",
"Y = list_bisec[1]\n",
"plt.plot(X, Y)\n",
"\n",
"plt.yscale(\"log\") # y軸を対数目盛に\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 丸め誤差:25点\n",
"\n",
"大きな数どおしのわずかな差は，丸め誤差にとくに影響を受ける．\n",
"1. 23.173-23.094を有効数字がそれぞれ５桁，４桁，３桁，２桁で計算した結果を示せ．\n",
"1. 同様に，0.81321/(23.173-23.094)を有効数字がそれぞれ５桁，４桁，３桁，２桁で計算した結果を示せ．\n",
"\n",
"(E.クライツィグ著「数値解析」(培風館,2003), p.10, 問題1.1-3改)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"from decimal import *\n",
"\n",
"def pretty_p(result,a,b,operator):\n",
" print('context.prec:{}'.format(getcontext().prec))\n",
" print(' %20.14f' % (a))\n",
" print( '%1s%20.14f' % (operator, b))\n",
" print('-----------')\n",
" print( ' %20.14f' % (result))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"context.prec:5\n",
" 23.17300000000000\n",
"- 23.09400000000000\n",
"-----------\n",
" 0.07900000000000\n",
"0.079\n",
"10.294\n"
]
}
],
"source": [
"getcontext().prec = 5\n",
"\n",
"a=Decimal('0.81321')\n",
"b=Decimal('23.173')\n",
"c=Decimal('23.094')\n",
"pretty_p(b-c,b,c,'-')\n",
"\n",
"print(b-c)\n",
"print(a/(b-c))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"context.prec:4\n",
" 23.17000000000000\n",
"- 23.09000000000000\n",
"-----------\n",
" 0.08000000000000\n",
"0.08\n",
"10.16\n"
]
}
],
"source": [
"TWOPLACES = Decimal(10) ** -2 \n",
"getcontext().prec = 4\n",
"a=Decimal('0.81321').quantize(Decimal(10) ** -4)\n",
"b=Decimal('23.173').quantize(Decimal('0.01'))\n",
"c=Decimal('23.094').quantize(Decimal('0.01'))\n",
"\n",
"pretty_p(b-c,b,c,'-')\n",
"\n",
"print(b-c)\n",
"print(a/(b-c))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"context.prec:3\n",
" 23.20000000000000\n",
"- 23.10000000000000\n",
"-----------\n",
" 0.10000000000000\n",
"0.1\n",
"8.13\n"
]
}
],
"source": [
"ONEPLACES = Decimal(10) ** -1\n",
"getcontext().prec = 3\n",
"a=Decimal('0.81321').quantize(Decimal(10) ** -3)\n",
"b=Decimal('23.173').quantize(ONEPLACES)\n",
"c=Decimal('23.094').quantize(ONEPLACES)\n",
"\n",
"pretty_p(b-c,b,c,'-')\n",
"\n",
"print(b-c)\n",
"print(a/(b-c))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"context.prec:2\n",
" 23.00000000000000\n",
"- 23.00000000000000\n",
"-----------\n",
" 0.00000000000000\n",
"0\n"
]
},
{
"ename": "DivisionByZero",
"evalue": "[]",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mDivisionByZero\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mpretty_p\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'-'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mDivisionByZero\u001b[0m: []"
]
}
],
"source": [
"ZEROPLACES = Decimal(10) ** 0\n",
"getcontext().prec = 2\n",
"a=Decimal('0.81321').quantize(Decimal(10) ** -2)\n",
"b=Decimal('23.173').quantize(ZEROPLACES)\n",
"c=Decimal('23.094').quantize(ZEROPLACES)\n",
"\n",
"pretty_p(b-c,b,c,'-')\n",
"\n",
"print(b-c)\n",
"print(a/(b-c))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"2桁の時には0割となっている．また，3桁でも相当大きなズレが出ていることが確認できる．似た同士の数値の引き算には注意が必要．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Newtonの差分商公式:25点\n",
"\n",
"次の4点の内挿式をNewtonの差分商公式を用いて求める．\n",
"``` python\n",
"X:=[-1, 0, 1, 2];\n",
"Y:=[4., -2., -1.2, -.52];\n",
"```\n",
"最初の３点を用いて求めた2次の補間式は，\n",
"\n",
"$$\n",
"- 2.00 - 6.00\\,x+ 3.40\\, \\left( x+1 \\right) x\n",
"$$\n",
"\n",
"である．４点目を取り入れて３次の補間式を求めよ．"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# https://stackoverflow.com/questions/14823891/newton-s-interpolating-polynomial-python\n",
"# by Khalil Al Hooti (stackoverflow)\n",
"\n",
"\n",
"def _poly_newton_coefficient(x,y):\n",
" \"\"\"\n",
" x: list or np array contanining x data points\n",
" y: list or np array contanining y data points\n",
" \"\"\"\n",
"\n",
" m = len(x)\n",
"\n",
" x = np.copy(x)\n",
" a = np.copy(y)\n",
" for k in range(1,m):\n",
" a[k:m] = (a[k:m] - a[k-1])/(x[k:m] - x[k-1])\n",
"\n",
" return a\n",
"\n",
"def newton_polynomial(x_data, y_data, x):\n",
" \"\"\"\n",
" x_data: data points at x\n",
" y_data: data points at y\n",
" x: evaluation point(s)\n",
" \"\"\"\n",
" a = _poly_newton_coefficient(x_data, y_data)\n",
" n = len(x_data) - 1 # Degree of polynomial\n",
" p = a[n]\n",
" for k in range(1,n+1):\n",
" p = a[n-k] + (x -x_data[n-k])*p\n",
" return p"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 4. -6. 3.4 -1.15333333]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9xvHvYhJCZFAwKFNUBkURMBQRBwhQRUSsszRUsUq0jvVqHUpbr49yHapWsVLF4daWIVerKHJVBAnigDLJIKKIlEEEJ0QNUab87h8rXhASEnJ2ss7Z5/08z35ycs4++/xWdnjZWXvttZ2ZISIi8VErdAEiIhItBbuISMwo2EVEYkbBLiISMwp2EZGYUbCLiMSMgl1EJGYU7CIiMaNgFxGJmTohPrRZs2aWnZ1dpfdu2rSJhg0bRltQIGpL8olLO0BtSVaJtGXevHlfmlnzitYLEuzZ2dnMnTu3Su+dMWMGffr0ibagQNSW5BOXdoDakqwSaYtzblVl1lNXjIhIzCjYRURiRsEuIhIzCnYRkZhRsIuIxExkwe6cq+2ce9c5Nzmqbf7EuHGQnU3vvn0hO9t/LyIiu4lyuOM1wFKgUYTb9MaNg/x8KC7GAaxa5b8HyMuL/ONERFJZJEfszrlWwKnAY1FsbzcjRkBxMS9zMndyo3+uuNg/LyIiP+GiuOepc+5fwB3AvsD1ZjaojHXygXyArKysnIKCgkpvv3ffvjgzfsfdjOJqNrAfDSnGnOO16dMTrj+UoqIiMjMzQ5cRibi0JS7tALUlWSXSltzc3Hlm1r3CFc0soQUYBIwufdwHmFzRe3JycmyvtG1rBjaVfgZm/8spZuCfT2GFhYWhS4hMXNoSl3aYqS3JKpG2AHOtErkcRVfMccBg59xKoADo65wbG8F2dxg5EjIyOJ43aEAxUzgZMjL88yIi8hMJB7uZ3WxmrcwsGzgfmG5mQxOubGd5eTBmDPXbtqA3r/FKnVNhzBidOBURKUPqjGPPy4OVKznk8vZ8sK0dq09QqIuIlCXSYDezGVbGidMode++AYBXXqnOTxERSV2pc8ReKju7mJYtYcqU0JWIiCSnlAt25+Ckk2DaNNi+PXQ1IiLJJ+WCHeDkk2HjRpgzJ3QlIiLJJyWDvX9/f+SufnYRkd2lZLDvvz90765+dhGRsqRksIPvZ3/nHfj669CViIgkl5QN9lNO8SdPp00LXYmISHJJ2WA/5hho0gReeil0JSIiySVlg71OHd8d8/LLEMEElSIisZGywQ6+O2bdOli4MHQlIiLJI6WDfcAA//XFF8PWISKSTFI62Fu0gG7d1M8uIrKzlA528N0xs2b5K1FFRCQmwb59O0ydGroSEZHkkPLB3rOnhj2KiOws5YNdwx5FRH4q5YMddgx7XLAgdCUiIuHFJtidg8mTQ1ciIhJeLII9Kwt69FCwi4hATIIdYNAgmD0bPvssdCUiImHFJthPPdV/1egYEUl3sQn2rl3hoIPUHSMiEptgd853x0yZAlu2hK5GRCSchIPdOVffOTfbObfQObfEOXdrFIVVxaBBUFQEM2eGqkBEJLwojtg3A33NrAvQFRjgnOsZwXb3Wr9+UL++umNEJL0lHOzmFZV+W7d0CXINaEYG5Ob6YNdVqCKSrpxFkIDOudrAPKAd8JCZ3VjGOvlAPkBWVlZOQUFBlT6rqKiIzMzMcl9/7rmDeOCBDjz55GzatCmu0mfUlIrakkri0pa4tAPUlmSVSFtyc3PnmVn3Clc0s8gWoAlQCBy5p/VycnKsqgoLC/f4+urVZmB2551V/ogaU1FbUklc2hKXdpipLckqkbYAc60SWRzpqBgz21ga7AOi3O7eaN0ajj4aJk0KVYGISFhRjIpp7pxrUvq4AfBz4INEt5uIwYP9zTd0FaqIpKMojtgPBAqdc4uAOcBUMws6LuX00/3JU42OEZF0VCfRDZjZIqBbBLVEpksXaNvWd8dcfHHoakREalZsrjzdmXO+O2bqVChO7oExIiKRi2Wwgw/277/XvVBFJP3ENth794bGjTU6RkTST2yDvW5dGDgQXngBtm8PXY2ISM2JbbCDHx3zxRd+6KOISLqIdbCfcgrUqwcTJ4auRESk5sQ62Bs1gp//HJ59VpOCiUj6iHWwA5x5JqxcCQsWhK5ERKRmxD7YBw+GWrX8UbuISDqIfbA3awYnnqhgF5H0EftgB98d8/778EHQqclERGpGWgT7L37hv2p0jIikg7QI9tatoUcPBbuIpIe0CHbw3TFz5sDq1aErERGpXmkV7KCTqCISf2kT7O3b+3nan3oqdCUiItUrbYId4Jxz/Lwxa9aErkREpPqkXbADPPNM2DpERKpTWgV7hw5w1FHw9NOhKxERqT5pFezgj9rfegs++SR0JSIi1SMtgx3UHSMi8ZV2wd6xI3TurO4YEYmvtAt28Eftb74Ja9eGrkREJHoJB7tzrrVzrtA5975zbolz7pooCqtO557rv+qoXUTiKIoj9m3AdWbWCegJXOGc6xTBdqtNx47QtSsUFISuREQkegkHu5mtM7P5pY+/A5YCLRPdbnU7/3x45x1YsSJ0JSIi0XIW4c1AnXPZwEzgSDP7dpfX8oF8gKysrJyCKh4uFxUVkZmZmVihwPr1+zBkyLFccskK8vLCzAwWVVuSQVzaEpd2gNqSrBJpS25u7jwz617himYWyQJkAvOAMytaNycnx6qqsLCwyu/dVa9eZp07R7a5vRZlW0KLS1vi0g4ztSVZJdIWYK5VIo8jGRXjnKsLPAOMM7OUmT/x/PNh8WJYsiR0JSIi0YliVIwDHgeWmtl9iZdUc84919/oWidRRSROojhiPw74FdDXObegdBkYwXarXVYW9O0LEyZAhKcaRESCimJUzBtm5szsKDPrWrq8GEVxNWHIEPj4Y5g7N3QlIiLRSMsrT3d2xhlQrx6MGxe6EhGRaKR9sDdtCoMG+X72bdtCVyMikri0D3aAoUPhs89g2rTQlYiIJE7BDgwc6I/c//nP0JWIiCROwQ7ss48f+jhxInz3XehqREQSo2AvNXQofP+9D3cRkVSmYC913HGQnQ1jx4auREQkMQr2Us75o/ZXX4VPPw1djYhI1SnYdzJ0KJSUwPjxoSsREak6BftOOnaEY46BJ5/UFAMikroU7LsYNgzeew/mzQtdiYhI1SjYd3H++X7449//HroSEZGqUbDvokkTP3/M+PGweXPoakRE9p6CvQzDhsHXX8OkSaErERHZewr2MvTvDy1bqjtGRFKTgr0MtWvDBRfAyy9rTLuIpB4FezmGDfNj2jUxmIikGgV7OTp08NMMPPGExrSLSGpRsO/B8OGwbBnMnBm6EhGRylOw78E550CjRvDoo6ErERGpPAX7HmRkQF4e/OtffvijiEgqULBXYPhwf6GSpvMVkVShYK9At26Qk+O7Y3QSVURSQSTB7px7wjn3uXPuvSi2l2yGD4fFi2H27NCViIhULKoj9r8DAyLaVtIZMsT3t48ZE7oSEZGKRRLsZjYT2BDFtpJRo0Y+3CdMgI0bQ1cjIrJnziLqOHbOZQOTzezIcl7PB/IBsrKycgoKCqr0OUVFRWRmZlaxyqpbtiyTSy/tzpVXfsRZZ62NZJuh2lId4tKWuLQD1JZklUhbcnNz55lZ9wpXNLNIFiAbeK8y6+bk5FhVFRYWVvm9ierRw+zww81KSqLZXsi2RC0ubYlLO8zUlmSVSFuAuVaJjNWomL1w+eWwdCm89lroSkREyqdg3wvnngtNm8Lo0aErEREpX1TDHScAs4COzrlPnHMXR7HdZNOgAVx0EUycCOvWha5GRKRsUY2KGWJmB5pZXTNrZWaPR7HdZHTZZbBtGzwe2xaKSKpTV8xeat8eTjoJHn4Ytm4NXY2IyO4U7FVw1VWwdq3vkhERSTYK9ioYOBAOPRQeeCB0JSIiu1OwV0GtWv6o/a23YO7c0NWIiPyUgr2Khg2DzEx48MHQlYiI/JSCvYoaN/ZDHwsK4LPPQlcjIrKDgj0BV14JW7bAI4+ErkREZAcFewI6dIBTTvFXov7wQ+hqREQ8BXuCrrvOd8WMHx+6EhERT8GeoL59oUsXuOceKCkJXY2IiII9Yc7B9df7WR9ffjl0NSIiCvZInHcetGwJ994buhIREQV7JOrWhWuugenTYf780NWISLpTsEdk+HB/wdI994SuRETSnYI9Ik2aQH4+PPUU/PvfoasRkXSmYI/Qtdf6eWT+/OfQlYhIOlOwR6hVK7jwQnjiCVi/PnQ1IpKuFOwRu+EGfwOO++8PXYmIpCsFe8Tat4dzzvHTDGzcGLoaEUlHCvZqcNNN8N138NBDoSsRkXSkYK8GXbv6uyzdfz8UFYWuRkTSjYK9mvzxj/Dll/C3v4WuRETSjYK9mvTsCSed5Ic+btoUuhoRSSeRBLtzboBz7kPn3HLn3E1RbDMObrkFvvhCR+0iUrMSDnbnXG3gIeAUoBMwxDnXKdHtxkGvXv6o/e67ddQuIjUniiP2HsByM1thZluAAuD0CLYbCzpqF5Ga5swssQ04dzYwwMwuKf3+V8AxZnblLuvlA/kAWVlZOQUFBVX6vKKiIjIzMxOquab97ndHsXx5JuPHv0ODBtv///lUbEt54tKWuLQD1JZklUhbcnNz55lZ9wpXNLOEFuBs4LGdvv8V8Nc9vScnJ8eqqrCwsMrvDWXWLDMwGznyp8+nYlvKE5e2xKUdZmpLskqkLcBcq0QuR9EVsxZovdP3rUqfk1I9e8Jpp/m+9g0bQlcjInEXRbDPAdo75w52ztUDzgcmRbDdWLn9dvj2W838KCLVL+FgN7NtwJXAFGAp8JSZLUl0u3Fz1FEwZAg88ACsWxe6GhGJs0jGsZvZi2bWwcwONbORUWwzjm691c/8OFI/IRGpRrrytAa1awcXXwxjxsDHH4euRkTiSsFew265xd/8+uabQ1ciInGlYK9hBx4Iv/sdPP00LFnSKHQ5IhJDdUIXkI6uvx4eeQT+9rdDufxycC50RSLxsm0brFoFa9bAJ5/45auv/M1vNm6EH37w62zbBnXqQP36fmncGA44AJo3h5Yt4dBD/ZJq10Yp2APIzITbboPhwxvz7LNw1lmhKxJJXZ9+CvPn+2XBAvjgA1i+3A9U2Fn9+tC0qQ/vBg18l2jt2rB9uw/677/3of/VV1BS8tP3tmwJ3brB0UfDz34Gxx8PTZrUXBv3loI9kIsugv/6ryJuvDGT006DevVCVySS/Mzgo49g+nR44w2/rFrlX3PO35qyUycYPBg6dIDsbH+T+ZYtoWHDyn3G9u3+QsI1a/wgh+XL4f33/X8cL77oQ985H/R9+/rP6tXL/yeRLBTsgdSuDZddtoIbbzyKUaN894yI7G7TJnj99WaMHQtTp8Lq1f75Aw/0R86//a0/iu7SJZouk9q1fVdM8+b+CH1nxcUwZw7MmOGXUaPgnnugWTMf8EOHQu/eUCvw2UsFe0A9emzg1FP9+Pa8PP+LKiLw9dfw3HPwzDMwbRps3nwkjRtDv35+RFn//r7vu6bPT2Vk+ODu3duPcPv2W5gyBZ5/3g+IeOIJaNsWLrgA8vP9XwshaFRMYH/5C2zerOGPIt9/DxMmwKBBkJUFv/41LFkCv/kN3HffAr74wgf9ZZf5a0KSYdBBo0ZwzjkwdiysXw/jx8Nhh/mLEA8+2F9t/vbbpSuPGwfZ2fTu29f3EY0bV211KdgDa98e/uM/4Mkn4Z13QlcjUrPM4K234JJLoEUL+OUvYdEiuOYa3+WxYoU/+OnWbSN164auds8yMnyQv/yy75e/+mp46SU49ljod8R6Xr/477BqFc7MnxjIz6+2cFewJ4ERI3w3zFVX7X42XiSOvv7a90937gzHHQcFBXDmmVBYCCtX+snyundPjqPyqjj4YLj3Xj/M8t57YckHtTlx81T6MY1FdPYrFRf7f/zVQMGeBPbd10/pO2cOPP546GpEqs9778Gll/q+52uu8Ue5jz3muzH++7+hT5/wJx6jlJnp/yJfUZLNfVzL+3SiZOfY/fFMcMRi9CNMbXl5/oTMjTfC55+HrkYkOma+S6J/f3+E/o9/wPnnw7x5MHu2nz8p1S4A2lsZbZtzLfezmjZ0ZeGOF9q0qZbPU7AnCef8fVGLiuC660JXI5K4LVv8UXjnzjBwoL9w6I47fPfE44/vPpQw1kaOhIwM6rJtx3MZGdU21auCPYkcfjjccIM/w/7qq6GrEama4mLff96unR/ZUqcO/POf/kToTTfB/vuHrjCAvDw/rWvbtphzfkzkmDH++WqgYE8yI0b48bm/+Y2/zFkkVWza5E96Zmf7/vPsbN8F8+67/sKdtL+6Oi8PVq7ktenT/Rniagp1ULAnnQYNYPRof9n07beHrkakYps2+ZP/2dn+L85u3WDmTL8MGJC6I1tSmYI9CZ10Elx4Idx5pz/aEUlGP/zgb/V4yCH+pH9Ojh+TPmUKnHBC6OrSm4I9Sd13n5+r4qKL/EkokWSxbZs/+dm+vZ+n5cgj4c03/YU5xx4bujoBBXvS2m8/ePhhWLgQ7rordDUiftjis8/6US6XXOJnTHz1Vb/06hW6OtmZgj2JnX66v0T5tttg8eLQ1Ug6e/NNH95nneX7zCdOhFmz/LS1knwU7Elu1Ch/c4ChQ/1kYSI1adkyf6n/8cf76U0efdTP5fKLX+ikaDJTsCe5Zs38VKCLFsEf/xi6GkkXGzb4/vMjjvBzoN92mx+pdcklfly6JLeEgt05d45zbolzrsQ51z2qouSnTj3Vz69xzz3w2muhq5E427rVj3Rp1w4efNBfYLR8OfzhD5W/A5GEl+gR+3vAmcDMCGqRPbj3Xv+P7YIL/H0ZRaL24ov+xOhvf+tnVlywwN90PSsrdGWytxIKdjNbamYfRlWMlK9hQz/VwNq1/ujdLHRFEhcffujncjn1VD9t9Asv+LHonTuHrkyqSn3sKaRHD3816lNP+SMpkUR8842/1+6P49DvvddPqztokE6MpjpnFRz6OeemAS3KeGmEmT1fus4M4Hozm7uH7eQD+QBZWVk5BQUFVSq4qKiIzJjM8VmVtpSUwM03d+bdd5syevR82rUrqqbq9k5c9ktc2gHlt6WkBF55pQVjxhzCxo11GThwHb/+9b/Zb7+tAaqsnHTYL5WRm5s7z8wqPp9pZgkvwAyge2XXz8nJsaoqLCys8nuTTVXb8vnnZgcdZNa+vdk330RbU1XFZb/EpR1mZbdlzhyznj3NwH+dM6fm66qKuO+XygLmWiUyVl0xKah5c38rsY8/9jcpUH+7VOTLL/0tNnv0gH//299j9803/UlSiZ9Ehzue4Zz7BDgW+F/n3JRoypKKnHCCn2rgX//SlAOyi3HjIDub3n37sr3tIYweNpsOHfz1ENde6y86uuCCeN2CTn4qoUsNzGwiMDGiWmQvXXedv73Y738PXbv6KVIlzY0b5w/Ni4t5m55csfoh3n3yaHI7refB11twxBGhC5SaoP+zU5hz/kbAnTv7OWU+/jh0RRLciBF8XtyQi3iCXszicw7gfziXV4t6KtTTiII9xTVs6CdkqlULTjtNFy+ls23b4MFVg+nAMsaRx43cyQccxrk8jVuzOnR5UoMU7DFwyCHwzDP+0u+zz/aXhUt6+fFE6NWMogezWUxn7uRmMtnkV2jTJmyBUqMU7DHRp4+fee/VV/39UjVSJj2sX+/vtnX88X7irqevnsmUBmfQkWU7VsrIgJEjwxUpNU7BHiMXXuhngHz8cX9bPYmvrVvh/vuhY0eYMAFuvhmWLoWzHzgR9+gYaNsWcw7atoUxY6r1xsmSfDQBZ8zceiusWOFHyjRv7qdZlXiZMQOuuspf/j9ggJ+NsUOHnVbIy4O8PF6bMYM+ffoEqlJCUrDHjHN+vPKGDX6ysKZN/V1vJPWtWePndnnqKcjOhueeg8GDNa+L7E5dMTFUr56/cKlnT/jlL2HatNAVSSJ++MF3kR92GEyaBP/5n/D++/7WiQp1KYuCPaYyMmDyZN8He/rp/s93SS1m/qi8Uyd/o4uTT/b96LfcAg0ahK5OkpmCPcaaNvW3NcvO9nNt6+5LqWPJEh/kZ5zh/5OeNg2efdbvS5GKKNhjLisLpk/3gTBwoMI92X31FVx5JXTpAnPm+BOjCxZAv36hK5NUomBPAz+Ge9u2cMop/hZokly2bPHDF9u3h4cfhssu8xecXX21bh4te0/BniaysvzR+uGH+z73ceNCVySwox/9iCP8zIs9esDChfDXv8L++4euTlKVgj2NNG8OhYX+KsWhQ2HUqNAVpbd33oHevX0/er168NJL8PLLaLIuSZiCPc00auQD5Iwz4Jpr/LJtW+iq0svy5XDeeX446rJlvutl4UJNuyzRUbCnofr14emn/Z/+o0b5i1y+/TZ0VfG3fj1ccYXvDps8Gf70J/joI38hmfrRJUoK9jRVuzbcd58/WnzlFejVyx89SvS+/tpP8XDooX7aluHD/dz5t94K++4bujqJIwV7mrv0UpgyxR9Ndu/ux0pLNL77Dm67DQ4+GO64w/9ltHQpjB4NLVqErk7iTMEu9OsH8+f7LoKzzvK33NuyJXRVqeu773yQZ2f77pY+fXwf+oQJ0K5d6OokHSjYBfD3YZg50/cB33cfHHOMv/pRKm/jRrj9dh/ov/89HHsszJ7thzMedVTo6iSdKNjl/+2zjx8//fzzsHYt5OT4i2ZKSkJXltw++8zPh96mjZ8P/8dAnzwZfvaz0NVJOlKwy24GD4bFi+HnP/cjZ447DhYtCl1V8vngA38itG1buOsuf1Xvu+8q0CU8BbuUKSvLTxE7dqwfwXH00XDDDb7/OJ2Z+YnVTjvNn5MYOxaGDfMh/z//A127hq5QJMFgd8792Tn3gXNukXNuonOuSVSFSXjO+ZvxLF3qw+vPf/ZzmYwZk34XNX3zDTz0kL8q9KSTfFfLn/4Eq1b5IaM/uYORSGCJHrFPBY40s6OAZcDNiZckyWb//eGxx+Dtt/2ojksv9UemzzwDJf8cB9nZ9O7b1581jNEkNGY+wIcPh4MO8rMuZmTAP/4Bq1f7cegHHBC6SpHdJRTsZvaKmf147PY20CrxkiRZHXMMvP66D/StW+Hss+GoC7tRsKon262WP3zNz0/5cP/0U7j7bn90fswxMH68vxPVnDkwdy786lf+RLNIsoqyj/3XwEsRbk+SkHNw5pn+1mzjm12NmTGEAg7lY+7iBr4sbgAjRoQuc699+SU88ghce20XWrWCG2/0NyoZM8YH/aOP+gu4RFKBM7M9r+DcNKCs6+RGmNnzpeuMALoDZ1o5G3TO5QP5AFlZWTkFBQVVKrioqIjMzMwqvTfZpHpbevftixlMYjAPchXT6cc+/MAZTKTLHUfQvfvX1Kmz59+vkD79tD5vvdWMN95oxuLFjSkpcbRsWUT//l/Sr99ntG79fegSE5Lqv187U1u83NzceWZW8SGGmSW0AMOAWUBGZd+Tk5NjVVVYWFjl9yablG9L27ZmvivaDGwxR9jl/NX2q7XBwKx5c7P8fLMXXjArLg5drNlXX5k9/7zZFVeYtW+/o/QjjzT7wx/M5s83mz69MHSZkUn536+dqC0eMNcqkbEJzSnnnBsA3AD0NrPiRLYlKWjkSN+nXux3/ZEs4aGMG/jL6P14qckQxo3z/dNjxvibL5944o4lJ6d6b8i8davvLpo3z/eNv/EGvPeefy0jA3Jz/cnQQYPgkEN2vE83/ZY4SHSy0L8C+wBTnXMAb5vZZQlXJakhL89/HTECW70a16YNjBxJvbwhnI6/U9PmzT4sJ0/2X3/sfq9VCzp29Pf27NTJT5R18MHQqhU0a+bD1/9Kla2kxA9B/OILWLPGj1JZudKPJ1+61M9UuXmzX3ffff3VoOedByec4OdB18lPibOEgt3MNKVRusvLg7w8Xpsxgz59+uz28j77wMkn+wX8zZrffNOPLlm4EGbNgrJOt9SvD40b+/fXq+enGd6yxS/FxX5ell3P5tSq5f9zOPxw/3nduvm/DNq396+JpAtN7y81av/9/ZQFgwfveO777/1IyRUr/AiUr77yyzff+CDfvBm2b/cBX6+e78LZbz+/NGsGrVv7pWVLHYmLgIJdkkCDBnDYYX4RkcTpD1QRkZhRsIuIxIyCXUQkZhTsIiIxo2AXEYkZBbuISMwo2EVEYkbBLiISMxVO21stH+rcF8CqKr69GfBlhOWEpLYkn7i0A9SWZJVIW9qaWfOKVgoS7Ilwzs21ysxHnALUluQTl3aA2pKsaqIt6ooREYkZBbuISMykYrCPCV1AhNSW5BOXdoDakqyqvS0p18cuIiJ7lopH7CIisgdJH+zOuXOcc0uccyXOuXLPJDvnBjjnPnTOLXfO3VSTNVaWc24/59xU59xHpV+blrPeSufcYufcAufc3JquszwV/YydN6r09UXOuaND1FkZlWhLH+fcN6X7YIFz7k8h6qyIc+4J59znzrn3ynk9lfZJRW1JlX3S2jlX6Jx7vzS7riljnerdL5W543XIBTgc6AjMALqXs05t4GPgEKAesBDoFLr2Muq8G7ip9PFNwF3lrLcSaBa63r39GQMDgZcAB/QE3glddwJt6QNMDl1rJdpyInA08F45r6fEPqlkW1JlnxwIHF36eF9gWU3/W0n6I3YzW2pmH1awWg9guZmtMLMtQAFwevVXt9dOB54sffwk8IuAteytyvyMTwf+Yd7bQBPn3IE1XWglpMrvS4XMbCawYQ+rpMo+qUxbUoKZrTOz+aWPvwOWAi13Wa1a90vSB3sltQTW7PT9J+z+g0wGWWa2rvTxeiCrnPUMmOacm+ecy6+Z0ipUmZ9xquyHytbZq/TP5Jecc0fUTGmRS5V9UlkptU+cc9lAN+CdXV6q1v2SFPc8dc5NA1qU8dIIM3u+putJxJ7asvM3ZmbOufKGJB1vZmudcwcAU51zH5QezUjNmQ+0MbMi59xA4DmgfeCa0l1K7RPnXCbwDPBbM/u2Jj87KYLdzPonuIm1QOudvm9V+lyN21NbnHOfOecONLN1pX92fV7ONtaWfv3cOTcR33UQOtgr8zM2EWJKAAABUElEQVROmv1QgQrr3Pkfopm96Jwb7ZxrZmapNl9JquyTCqXSPnHO1cWH+jgze7aMVap1v8SlK2YO0N45d7Bzrh5wPjApcE1lmQRcWPr4QmC3v0accw2dc/v++Bg4CShzlEANq8zPeBJwQekZ/57ANzt1PSWTCtvinGvhnHOlj3vg/618VeOVJi5V9kmFUmWflNb4OLDUzO4rZ7Xq3S+hzyBX4gzzGfj+p83AZ8CU0ucPAl7c5SzzMvxohxGh6y6nLfsDrwIfAdOA/XZtC36kxsLSZUkytaWsnzFwGXBZ6WMHPFT6+mLKGcWUDEsl2nJl6c9/IfA20Ct0zeW0YwKwDtha+u/k4hTeJxW1JVX2yfH482SLgAWly8Ca3C+68lREJGbi0hUjIiKlFOwiIjGjYBcRiRkFu4hIzCjYRURiRsEuIhIzCnYRkZhRsIuIxMz/ATlv2JnTkJHKAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"from scipy import interpolate\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# もとの点\n",
"x = np.array([-1,0,1,2])\n",
"y = np.array([4,-2,-1.2,-0.52])\n",
"for i in range(0,4):\n",
" plt.plot(x[i],y[i],'o',color='r')\n",
"\n",
"\n",
"print(_poly_newton_coefficient(x,y))\n",
"\n",
"xx = np.linspace(-1,2, 100)\n",
"yy = newton_polynomial(x, y, xx)\n",
"plt.plot(xx, yy, color = 'b')\n",
"\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Newtonの差分商公式の利点であるデータ点の追加に伴う煩雑さの回避は，このコードでは実感できない．pythonは，やっぱり便利なのかな．"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ページランク:25点\n",
"\n",
"次のようなリンクが張られたページ群のページランクを求めよ．\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([[ 0.000, 0.500, 0.000, 0.000, 0.500],\n",
" [ 0.333, 0.000, 0.000, 0.000, 0.500],\n",
" [ 0.333, 0.500, 0.000, 0.500, 0.000],\n",
" [ 0.333, 0.000, 1.000, 0.000, 0.000],\n",
" [ 0.000, 0.000, 0.000, 0.500, 0.000]])\n"
]
}
],
"source": [
"from pprint import pprint\n",
"from numpy import array, zeros, diagflat, dot, transpose, set_printoptions\n",
"from scipy.linalg import eig\n",
"\n",
"A = array([[0,1,1,1,0],\n",
" [1,0,1,0,0],\n",
" [0,0,0,1,0],\n",
" [0,0,1,0,1],\n",
" [1,1,0,0,0]])\n",
"\n",
"n = 5\n",
"diag = []\n",
"for i in range(0,n):\n",
" tmp = 0.0\n",
" for j in range(0,n):\n",
" tmp += A[i,j]\n",
" diag.append(1.0/tmp)\n",
"\n",
"D = diagflat(diag)\n",
"tA = dot(transpose(A),D)\n",
"\n",
"set_printoptions(formatter={'float': '{: 0.3f}'.format}) \n",
"pprint(tA)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"初期ベクトルを等分の値にして，3度ほどホップさせた結果．"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([ 0.125, 0.111, 0.269, 0.328, 0.167])\n"
]
}
],
"source": [
"x = array([1/5,1/5,1/5,1/5,1/5])\n",
"pprint(dot(tA,dot(tA,dot(tA,x))))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"固有値を求める．固有値がソートされているか自信がないので，表示させてみた．"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([ 1.00000000+0.j , 0.10242347+0.50115386j,\n",
" 0.10242347-0.50115386j, -0.81317748+0.j , -0.39166946+0.j ])\n"
]
}
],
"source": [
"l, V = eig(tA)\n",
"pprint(l)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"l[0]に対応する最大固有値のベクトルを取り出す．\n",
"さらに，初期ベクトルからのホップと比較するために値を揃えている．\n",
"だいたい一致しているが，ホップ数が少ないので一致はそれほど高くない．"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"array([ 0.29497881-0.j, 0.26220338-0.j, 0.55718219-0.j, 0.65550846-0.j,\n",
" 0.32775423-0.j])\n"
]
}
],
"source": [
"v0 = V[:,0]\n",
"pprint(v0[0]/0.294*v0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"これより，ページランクは，\n",
"[4, 3, 5, 1, 2]\n",
"の順になる．"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"?eig"
]
},
{
"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
}
},
"nbformat": 4,
"nbformat_minor": 2
}