{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#概要:補間と近似\" data-toc-modified-id=\"概要:補間と近似-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>概要:補間と近似</a></span></li><li><span><a href=\"#多項式補間(polynomial-interpolation)\" data-toc-modified-id=\"多項式補間(polynomial-interpolation)-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>多項式補間(polynomial interpolation)</a></span></li><li><span><a href=\"#Lagrange(ラグランジュ)-の内挿公式\" data-toc-modified-id=\"Lagrange(ラグランジュ)-の内挿公式-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Lagrange(ラグランジュ) の内挿公式</a></span><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><span><a href=\"#python-code\" data-toc-modified-id=\"python-code-3.0.1\"><span class=\"toc-item-num\">3.0.1&nbsp;&nbsp;</span>python code</a></span></li></ul></li></ul></li><li><span><a href=\"#Newton(ニュートン)-の差分商公式\" data-toc-modified-id=\"Newton(ニュートン)-の差分商公式-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Newton(ニュートン) の差分商公式</a></span><ul class=\"toc-item\"><li><span><a href=\"#Newton補間と多項式補間の一致の検証\" data-toc-modified-id=\"Newton補間と多項式補間の一致の検証-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Newton補間と多項式補間の一致の検証</a></span></li></ul></li><li><span><a href=\"#数値積分-(Numerical-integration)\" data-toc-modified-id=\"数値積分-(Numerical-integration)-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>数値積分 (Numerical integration)</a></span><ul class=\"toc-item\"><li><span><a href=\"#中点則-(midpoint-rule)\" data-toc-modified-id=\"中点則-(midpoint-rule)-5.1\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span>中点則 (midpoint rule)</a></span></li><li><span><a href=\"#台形則-(trapezoidal-rule)\" data-toc-modified-id=\"台形則-(trapezoidal-rule)-5.2\"><span class=\"toc-item-num\">5.2&nbsp;&nbsp;</span>台形則 (trapezoidal rule)</a></span></li><li><span><a href=\"#Simpson(シンプソン)則\" data-toc-modified-id=\"Simpson(シンプソン)則-5.3\"><span class=\"toc-item-num\">5.3&nbsp;&nbsp;</span>Simpson(シンプソン)則</a></span></li></ul></li><li><span><a href=\"#数値積分のコード\" data-toc-modified-id=\"数値積分のコード-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>数値積分のコード</a></span><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><span><a href=\"#Midpoint-rule(中点法)\" data-toc-modified-id=\"Midpoint-rule(中点法)-6.0.0.1\"><span class=\"toc-item-num\">6.0.0.1&nbsp;&nbsp;</span>Midpoint rule(中点法)</a></span></li><li><span><a href=\"#Trapezoidal-rule(台形公式)\" data-toc-modified-id=\"Trapezoidal-rule(台形公式)-6.0.0.2\"><span class=\"toc-item-num\">6.0.0.2&nbsp;&nbsp;</span>Trapezoidal rule(台形公式)</a></span></li><li><span><a href=\"#Simpson's-rule(シンプソンの公式)\" data-toc-modified-id=\"Simpson's-rule(シンプソンの公式)-6.0.0.3\"><span class=\"toc-item-num\">6.0.0.3&nbsp;&nbsp;</span>Simpson's rule(シンプソンの公式)</a></span></li></ul></li></ul></li></ul></li><li><span><a href=\"#課題\" data-toc-modified-id=\"課題-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>課題</a></span><ul class=\"toc-item\"><li><span><a href=\"#補間法\" data-toc-modified-id=\"補間法-7.1\"><span class=\"toc-item-num\">7.1&nbsp;&nbsp;</span>補間法</a></span></li><li><span><a href=\"#対数関数のニュートンの差分商補間(2014期末試験，25点)\" data-toc-modified-id=\"対数関数のニュートンの差分商補間(2014期末試験，25点)-7.2\"><span class=\"toc-item-num\">7.2&nbsp;&nbsp;</span>対数関数のニュートンの差分商補間(2014期末試験，25点)</a></span><ul class=\"toc-item\"><li><span><a href=\"#差分商補間の表中の開いている箇所[-XXX-]を埋めよ．\" data-toc-modified-id=\"差分商補間の表中の開いている箇所[-XXX-]を埋めよ．-7.2.1\"><span class=\"toc-item-num\">7.2.1&nbsp;&nbsp;</span>差分商補間の表中の開いている箇所[ XXX ]を埋めよ．</a></span></li><li><span><a href=\"#ニュートンの二次多項式\" data-toc-modified-id=\"ニュートンの二次多項式-7.2.2\"><span class=\"toc-item-num\">7.2.2&nbsp;&nbsp;</span>ニュートンの二次多項式</a></span></li><li><span><a href=\"#ニュートンの三次多項式の値を求めよ．\" data-toc-modified-id=\"ニュートンの三次多項式の値を求めよ．-7.2.3\"><span class=\"toc-item-num\">7.2.3&nbsp;&nbsp;</span>ニュートンの三次多項式の値を求めよ．</a></span></li></ul></li><li><span><a href=\"#数値積分(I)\" data-toc-modified-id=\"数値積分(I)-7.3\"><span class=\"toc-item-num\">7.3&nbsp;&nbsp;</span>数値積分(I)</a></span></li></ul></li><li><span><a href=\"#解答例[7.3]\" data-toc-modified-id=\"解答例[7.3]-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>解答例[7.3]</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br />\n",
    "\n",
    "<div style=\"text-align: center;\">\n",
    "<font size=\"7\">補間(interpolation)と数値積分(Integral)</font>\n",
    "</div>\n",
    "<br />\n",
    "<div style=\"text-align: right;\">\n",
    "<font size=\"4\">file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/interpolationintegral</font>\n",
    "<br />\n",
    "<font size=\"4\">https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python</font>\n",
    "<br />\n",
    "<font size=\"4\">cc by Shigeto R. Nishitani 2017 </font>\n",
    "</div>\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 概要:補間と近似\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "単純な2次元データについて補間と近似を考える．補間はたんに点を「滑らかに」つなぐことを，近似はある関数にできるだけ近くなるように「フィット」することを言う．補間はIllustratorなどのドロー系ツールで曲線を引くときの，ベジエやスプライン補間の基本となる．本章では補間とそれに密接に関連した積分について述べる．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAEICAYAAACUIhp/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWXElEQVR4nO3dfYxldX3H8c9HForD0BKFToRl5tqK1pVW7Wx9og3MRpvlSarRCl4xNLaTNNLQRAqYtVVjV1pjTDU+NFvRtTo6Ep+KC4VidhaK1equomVdaYnOLCtYBEQZrmKQb/84Z9mZ6ezOvb9z7jn3zLxfyYlzzpz7m++ey378zH1aR4QAAADQuyfVPQAAAEBTUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaRQK9u7bP9p4m1Hbc/bPqrsuQCsfnl+/EYP55M5+H8oUquU7VnbL+3ivOQiU7Wlf6aI2B8RwxHxyzrnAjD48vz4WV6E5m3PS3pmRHwv//5223+7zG3IHBwRRQqF8JsZgAY5Py9CB7d76h4IzUeRWuVsX2L7Ntvvtv1j29+3fXb+va2S/kDS+/Pf0N6fH/8t2zfbftD2nbb/eMF6221/yPYNth+RNJEf+8f8Ng/bvsX22ILbvMT2123/JP/flxxm1t+0vdP2A7bvtz1l+4T8ex+XNCrpi/msV9hu2Q7b6/JzTrZ9XT73Xbb/bMHab7N9re1/zmfca3tj2dcbQHPk+fEM25OS2pKuyPPli11mzi7b77D95TxX/s32iQvWf73tuTzT/rrbZwrQLBSpteGFku6UdKKkd0m6xrYjYoukf5d0af7b2aW2j5N0s6RPSvp1SRdK+qDtDQvWe62krZKOl3Rbfqwt6R35z7hd0pQk2X6KpOslvU/SUyW9R9L1tp+6zJyWdLWkkyU9W9Kpkt4mSRFxsaT9OvQb5buWuf20pAP57V8l6Z22Ny34/svzc06QdJ2k9x/pogFYGyJim7LMeleeL+d3mTlSlod/oiwvj5F0uSTlmflBZdn4NEm/JumU/v5JUAeK1NowFxH/lD+v/zFlf6lHDnPueZJmI+KjEfFYRHxT0mclvXrBOf8SEV+OiMcj4uf5sesj4taIeFTSFkkvtn2qpHMl/U9EfDxf71OSvivp/KU/OCLuioibI+LRiPiRstJ1Zjd/wPxnnSHpyoj4eUTcLunDkl6/4LTbIuKG/Dp8XNJzu1kbwKrxBdsP5dsXSlrzoxHx3xHxM0nXSnpefvxVkr4YEbdFxC8k/Y0k/nHbVWhd3QOgEj88+EVEdGxL0vBhzh2T9ELbDy04tk5Z8Tjo7mVu98SxiJi3/aCyR4ZOljS35Nw5LfObme0RSe9V9nTj8cqK/o8PM+dSJ0t6MCIeXvJzFj5998MFX3ckHWt7XUQ81uXPANBsfxQRXzq4Y7uMYrM0Vw5m68lanIsd2w+U8PMwYHhECkuD5G5Jt0TECQu24Yj48yPcRsqehpMk2R6W9BRJ9+Tb2JJzRyX9YJk13pmv/dsR8auSXqfs6b4j/dyD7pH0FNvHd/FzAGCp5fKlSNG6V9L6gzu2n6zs5Q1YZShS+F9JCz9HZYekZ9q+2PbR+fZ7tp+9wjrn2P5928coe63UVyPibkk35Ou91vY626+RtCH/OUsdL2le0k9snyLpr1aY9Qn5z/oPSVfbPtb270h6g6RPrDA3AEjL58thM6cLn5F0fv5mm2OUvd7TR74JmogihfdKelX+jr735U+N/aGyF5nfo+xh67+X9CsrrPNJSW+V9KCkcWWPJikiHlD2uqs3SXpA0hWSzouI+5dZ4+2SflfST5S9QP1zS75/taS35K9vuHyZ218kqZXP/XlJb134MD4AHME1kjYsef3USplzWBGxV9JfKHuDy73Kfkm8T9Kj5Y2MQeAIXvuGYmxvl3QgIt5S9ywAMIjylzw8JOm0iPh+zeOgRDwiBQBAH9g+3/ZQ/rEy75b0X5Jm650KZaNIAQDQHxfo0JtuTpN0YfA00KrDU3sAAACJeEQKAAAgUS0fyHniiSdGq9Xq6txHHnlExx13XH8H6pOmzs7c1Vorc+/Zs+f+iDipjyNVopf8ktbO/TsomLtaa2XuI+ZXRFS+jY+PR7dmZma6PnfQNHV25q7WWplb0u6oIW/K3nrJr4i1c/8OCuau1lqZ+0j5xVN7AAAAiShSAAAAiShSAAAAiShSAAAAiShSAAAAiUorUraPsv1N2zvKWhMAqkKGAUhR5iNSl0naV+J6AFAlMgxAz0opUrbXSzpX0ofLWA8AqkSGAUhVyr+1Z/szkq6WdLykyyPivGXOmZQ0KUkjIyPj09PTXa09Pz+v4eHhwjPW4dN75/Wa5zRv9qZec+auVq9zT0xM7ImIjX0cKdlKGZaaX1Jz71/yq1rMXa1S8+twn9TZ7SbpPEkfzL8+S9KOlW6zVj7ZfOzKHXWPkKSp15y5q7VaPtm81wxbK59sTn5Vi7mrNWifbH6GpJfbnpU0LWmT7U+UsC4AVIEMA5CscJGKiDdHxPqIaEm6UNLOiHhd4ckAoAJkGIAi+BwpAACARKUWqYjYFcu80Bzoq6kpqdXSmZs2Sa1Wtg8kIMMA9Gpd3QMAhUxNSZOTUqcjS9LcXLYvSe12nZMBANYAntpDs23ZInU6i491OtlxAAD6jCKFZtu/v7fjAACUiCKFZhsd7e04AAAlokih2bZulYaGFh8bGsqOAwDQZxQpNFu7LW3bJo2NKWxpbCzb54XmAIAKUKTQfO22NDurW3bulGZnKVEAgMpQpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpHDI1JTUaunMTZukVivbR/9wvQE0Ffn1hHV1D4ABMTUlTU5KnY4sSXNz2b4ktdt1TrY6cb0BNBX5tQiPSCGzZYvU6Sw+1ulkx1E+rjeApiK/FqFIIbN/f2/HUQzXG0BTkV+LUKSQGR3t7TiK4XoDaCryaxGKFDJbt0pDQ4uPDQ1lx1E+rjeApiK/FqFIIdNuS9u2SWNjClsaG8v21+ALByvB9QbQVOTXIhQpHNJuS7OzumXnTml2ds3+pagM1xtAU5FfTyhcpGwfa/trtr9le6/tt5cxGABUgQwDUEQZnyP1qKRNETFv+2hJt9n+14j4aglrA0C/kWEAkhUuUhERkubz3aPzLYquCwBVIMMAFOEsQwouYh8laY+kZ0j6QERcucw5k5ImJWlkZGR8enq6q7Xn5+c1PDxceMY6XHLjI9q++bi6x+hZU685c1er17knJib2RMTGPo6UbKUMS80vqbn3L/lVLeauVqn5FRGlbZJOkDQj6fQjnTc+Ph7dmpmZ6frcQTN25Y66R0jS1GvO3NXqdW5Ju6PEvOnH1k2G9ZJfEc29f8mvajF3tcrMr1LftRcRD+UhtLnMdQGgCmQYgF6V8a69k2yfkH/9ZEkvk/TdousCQBXIMABFlPGuvadJ+lj+GoMnSbo2InaUsC4AVIEMA5CsjHftfVvS80uYBQAqR4YBKIJPNgcAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhEkQIAAEhUuEjZPtX2jO3v2N5r+7IyBgOAKpBhAIpYV8Iaj0l6U0R8w/bxkvbYvjkivlPC2gDQb2QYgGSFH5GKiHsj4hv51w9L2ifplKLrAkAVyDAARTgiylvMbkm6VdLpEfHTJd+blDQpSSMjI+PT09NdrTk/P6/h4eHSZqzSJTc+ou2bj6t7jJ419Zozd7V6nXtiYmJPRGzs40iFHS7DUvNLau79S35Vi7mrVWp+RUQpm6RhSXskvXKlc8fHx6NbMzMzXZ87aMau3FH3CEmaes2Zu1q9zi1pd5SUN/3Yus2wXvIrorn3L/lVLeauVpn5Vcq79mwfLemzkqYi4nNlrAkAVSHDAKQq4117lnSNpH0R8Z7iIwFAdcgwAEWU8YjUGZIulrTJ9u35dk4J6wJAFcgwAMkKf/xBRNwmySXMAgCVI8NQq6kpacsWnbl/vzQ6Km3dKrXbdU+FHpTxOVIAAKBXU1PS5KTU6WRNfm4u25coUw3CPxEDAEAdtmyROp3Fxzqd7DgagyIFAEAd9u/v7TgGEkUKAIA6jI72dhwDiSIFAEAdtm6VhoYWHxsayo6jMShSAADUod2Wtm2TxsYUtjQ2lu3zQvNGoUgBAFCXdluandUtO3dKs7OUqAaiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSAEAACSiSPXD1JTUamVft1rZPgCgf/LcPXPTJnIXlVpX9wCrztSUNDkpdTrZ/txcti9J7XZ9cwHAarUgdy2Ru6gUj0iVbcuWQyXqoE4nOw4AKB+5ixpRpMq2f39vxwEAxZC7qBFFqmyjo70dBwAUQ+6iRhSpsm3dKg0NLT42NJQdBwCUj9xFjShSZWu3pW3bpLGxbH9sLNvnBY8A0B8LcjdscheVokj1Q7stzc5mX8/O8pcZAPotz91bdu4kd1GpUoqU7Y/Yvs/2HWWsBwBVIb8AFFHWI1LbJW0uaS0AqNJ2kV8AEpVSpCLiVkkPlrEWAFSJ/AJQhCOinIXslqQdEXH6Yb4/KWlSkkZGRsanp6e7Wnd+fl7Dw8OlzFi1S258RNs3H1f3GD1r6jVn7mr1OvfExMSeiNjYx5GS9Su/pObev+RXtZi7WqXmV0SUsklqSbqjm3PHx8ejWzMzM12fO2jGrtxR9whJmnrNmbtavc4taXeUlDdlb/3Kr4jm3r/kV7WYu1pl5hfv2gMAAEhEkQIAAEhU1scffErSVyQ9y/YB228oY10A6DfyC0AR68pYJCIuKmMdAKga+QWgCJ7aAwAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAASESRAgAcMjUltVrZ161Wtg/gsNbVPQAAYEBMTUmTk1Knk+3PzWX7ktRu1zcXMMB4RAoAkNmy5VCJOqjTyY4DWBZFCgCQ2b+/t+MAKFIAgNzoaG/HAVCkAAC5rVuloaHFx4aGsuMAlkWRAgBk2m1p2zZpbCzbHxvL9nmhOXBYFCkAwCHttjQ7m309O0uJAlZQSpGyvdn2nbbvsn1VGWsCQFXIMACpChcp20dJ+oCksyVtkHSR7Q1F1wWAKpBhAIoo4xGpF0i6KyK+FxG/kDQt6YIS1gWAKpBhAJKV8cnmp0i6e8H+AUkvXHqS7UlJk5I0MjKiXbt2dbX4p/fO65Ibry8+ZU1aVzV09qZec+auzNmnhqRddY9RhhUzLDW/pGZnGPlVMeauTJn5Vdk/ERMR2yRtk6SNGzfGWWed1eUtd+lDb+z23MHSuup6zf7duXWP0bNdu3ap+/tncDB3tZo6d4r0/JKammHkV7WYu1plzl3GU3s/kHTqgv31+TEAaAIyDECyMorU1yWdZvvpto+RdKGk60pYFwCqQIYBSFb4qb2IeMz2pZJuknSUpI9ExN7CkwFABcgwAEWU8hqpiLhB0g1lrAUAVSPDAKTik80BAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASUaQAAAASFSpStl9te6/tx21vLGsoAKgCGQagqKKPSN0h6ZWSbi1hFgCoGhkGoJB1RW4cEfskyXY50wBAhcgwAEU5IoovYu+SdHlE7D7COZOSJiVpZGRkfHp6uqu15+fnNTw8XHjGOnx677xe85zmzd7Ua87c1ep17omJiT0RMZBPn62UYan5JTX3/iW/qsXc1So1vyLiiJukLyl7+HvpdsGCc3ZJ2rjSWge38fHx6NbMzEzX5w6aps7O3NVaK3NL2h1dZkSZW9kZ1kt+pVynQcHc1WLuapWZXys+tRcRL+26sgHAgCHDAPQTH38AAACQqOjHH7zC9gFJL5Z0ve2byhkLAPqPDANQVNF37X1e0udLmgUAKkWGASiKp/YAAAASUaQAAAASUaQAAAASUaQAAAASlfLJ5j3/UPtHkua6PP1ESff3cZx+aurszF2ttTL3WESc1K9hqtJjfklr5/4dFMxdrbUy92Hzq5Yi1Qvbu2NA/1mJlTR1duauFnOvbk29TsxdLeauVplz89QeAABAIooUAABAoiYUqW11D1BAU2dn7mox9+rW1OvE3NVi7mqVNvfAv0YKAABgUDXhESkAAICBRJECAABI1IgiZfvVtvfaftz2wL/N0vZm23favsv2VXXP0y3bH7F9n+076p6lW7ZPtT1j+zv5fyOX1T1Tt2wfa/trtr+Vz/72umfqlu2jbH/T9o66Z2mCJmUY+VWtpmZYk/NLKjfDGlGkJN0h6ZWSbq17kJXYPkrSBySdLWmDpItsb6h3qq5tl7S57iF69JikN0XEBkkvkvTGBl3vRyVtiojnSnqepM22X1TvSF27TNK+uodokEZkGPlVi6ZmWJPzSyoxwxpRpCJiX0TcWfccXXqBpLsi4nsR8QtJ05IuqHmmrkTErZIerHuOXkTEvRHxjfzrh5X9xTil3qm6E5n5fPfofBv4d3/YXi/pXEkfrnuWpmhQhpFfFWtqhjU1v6TyM6wRRaphTpF094L9A2rAX4rVwHZL0vMl/WfNo3Qtf3j5dkn3Sbo5Ipow+z9IukLS4zXPgfKRXzVqWoY1NL+kkjNsYIqU7S/ZvmOZrRG/DaFetoclfVbSX0bET+uep1sR8cuIeJ6k9ZJeYPv0mkc6ItvnSbovIvbUPcugIcNQRBMzrGn5JfUnw9aVtVBREfHSumcoyQ8knbpgf31+DH1i+2hlATQVEZ+re54UEfGQ7Rllr/EY5BfLniHp5bbPkXSspF+1/YmIeF3Nc9VulWQY+VWDpmdYg/JL6kOGDcwjUqvI1yWdZvvpto+RdKGk62qeadWybUnXSNoXEe+pe55e2D7J9gn510+W9DJJ3611qBVExJsjYn1EtJT9t72TErWqkF8Va2qGNTG/pP5kWCOKlO1X2D4g6cWSrrd9U90zHU5EPCbpUkk3KXvR4LURsbfeqbpj+1OSviLpWbYP2H5D3TN14QxJF0vaZPv2fDun7qG69DRJM7a/rez/wG6OCD5OYBVqSoaRX7VoaoaRXzn+iRgAAIBEjXhECgAAYBBRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABJRpAAAABL9HyA2Y9pdfSL8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "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",
    "\n",
    "fig, (axL, axR) = plt.subplots(ncols=2, figsize=(10,4))\n",
    "\n",
    "x = np.array([0,1,2,3])\n",
    "y = np.array([1.2,3.2,3.8,3.2])\n",
    "for i in range(0,4):\n",
    "    axL.plot(x[i],y[i],'o',color='r')\n",
    "axL.hlines(0, -1, 4, linewidth=1)\n",
    "axL.vlines(0, -1, 4, linewidth=1)\n",
    "axL.set_title('Interpolation')\n",
    "axL.grid(True)\n",
    "\n",
    "x = np.array([0,1,2,3])\n",
    "y = np.array([0.2,1.2,1.8,3.2])\n",
    "for i in range(0,4):\n",
    "    axR.plot(x[i],y[i],'o',color='r')\n",
    "axR.hlines(0, -1, 4, linewidth=1)\n",
    "axR.vlines(0, -1, 4, linewidth=1)\n",
    "axR.set_title('Fitting')\n",
    "axR.grid(True)\n",
    "\n",
    "# fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 多項式補間(polynomial interpolation)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "データを単純に多項式で補間する方法を先ず示そう．$N+1$点を$N$次の多項式でつなぐ．この場合の補間関数は，\n",
    "\n",
    "$$\n",
    "F \\left(x \\right)={\\sum_{i=0}^{N} } a _{i }x ^{i }=a_{0}+a_{1}x +a_{2}x^{2}+\\cdots +a_{N}x^{N}\n",
    "$$\n",
    "である．データの点を$(x_{i},\\,y_{i}),i=0..N$とすると\n",
    "\n",
    "$$\n",
    "\\begin{array}{cl}\n",
    "a _{0}+a _{1}x _{0}+a _{2}x _{0}^{2}+\\cdots +a _{N }x _{0}^{N }& =y _{0}\\\\\n",
    "a _{0}+a _{1}x _{1}+a _{2}x _{1}^{2}+\\cdots +a _{N }x _{1}^{N }& =y _{1}\\\\\n",
    "\\vdots& \\\\\n",
    "a _{0}+a _{1}x _{N}+a _{2}x _{N}^{2}+\\cdots +a _{N }x _{N}^{N }& =y _{N}\n",
    "\\end{array}\n",
    "$$\n",
    "が，係数 $a_i$を未知数と見なした線形の連立方程式となっている．\n",
    "連立方程式の係数行列$\\mathbf{A}$(ヴァンデルモンド行列と呼ばれる)は\n",
    "$$\n",
    "A=\\left[\n",
    "\\begin{array}{ccccc}\n",
    "1&x_0&x_0^2&\\cdots&x_0^N \\\\\n",
    "1&x_1&x_1^2&\\cdots&x_1^N \\\\\n",
    "\\vdots& & & & \\vdots \\\\\n",
    "1&x_N&x_N^2&\\cdots&x_N^N \n",
    "\\end{array} \\right]\n",
    "$$\n",
    "となる．\n",
    "係数$a_i$および\n",
    "データ$y_i$をそれぞれベクトルとみなすと\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "$A$と$y$から$a_i$を導出：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "により未知数ベクトル$a_i$が求まる．これは単純に，前に紹介した Gauss の消去法や LU 分解で解ける．\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lagrange(ラグランジュ) の内挿公式\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "多項式補間は手続きが簡単であるため，計算間違いが少なく，プログラムとして組むのに適している．しかし，あまり\"みとうし\"のよい方法とはいえない．その点，Lagrange(ラグランジュ)の内挿公式は見通しがよい．これは\n",
    "\n",
    "$$\n",
    "F(x)= \\sum_{k=0}^N \\frac{ \\prod_{j \\ne k} (x-x_j)}\n",
    "{ \\prod_{j \\ne k} (x_k-x_j)} y_k\n",
    "=\\sum_{k=0}^N \\frac{ \\frac{ (x-x_0)(x-x_1)\\cdots(x-x_N)}{ (x-x_k)}}\n",
    "{\\frac{ (x_k-x_0)(x_k-x_1)\\cdots(x_k-x_N)}{ (x_k-x_k)}} y_k\n",
    "$$\n",
    "と表わされる．数学的に 2つ目の表記は間違っているが，先に割り算を実行すると読み取って欲しい．これは一見複雑に見えるが，単純な発想から出発している．求めたい関数$F(x)$を\n",
    "\n",
    "$$\n",
    "F(x) = y_0 L_0(x)+y_1 L_1(x)+y_2 L_2(x)\n",
    "$$\n",
    "とすると\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "L_0(x_0) = 1 &L_0(x_1) = 0 &L_0(x_2) = 0 \\\\\n",
    "L_1(x_0) = 0 &L_1(x_1) = 1 &L_1(x_2) = 0 \\\\\n",
    "L_2(x_0) = 0 &L_2(x_1) = 0 &L_2(x_2) = 1 \n",
    "\\end{array}\n",
    "$$\n",
    "となるように関数$L_i(x)$を決めればよい．これを以下のようにとればLagrangeの内挿公式となる．\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "$L_0(x)$：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "$L_1(x)$：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "$L_2(x)$：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "$F(x)$：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### python code\n",
    "pythonではLagrange補間はinterpolate.lagrangeで用意されている．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    3     2\n",
      "-1 x + 3 x - 1 x + 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhS0lEQVR4nO3deXhV1bnH8e+CADJYBJE4QVCvgigiBhHFgakIXgQURL1RsUrjUKvV1kqNs1KvWq21+qi0jhVFigiKKCgELCgUcJ6vAwgVCygOAZnf+8cbK0KAQ84+2Wef8/s8z3kynez9Lk74ZWftNQQzQ0REkqtW3AWIiEh6FOQiIgmnIBcRSTgFuYhIwinIRUQSriCOkzZr1sxatWpVre9dsWIFDRs2jLagmKgt2SdX2gFqS7ZKpy3z5s1bZma7bPr5WIK8VatWzJ07t1rfO23aNLp27RptQTFRW7JPrrQD1JZslU5bQggLqvq8ulZERBJOQS4iknAKchGRhFOQi4gknIJcRCThFOQiIgmnIBcRSbhYxpGL5LOvvoL582HpUn988QWsWgVr1vijbl3YYQeoXx923hkKC2HXXaGoyD8nsikFuUiGfPMNzJkDr7/uj7fego8/9iCvrpYtoXVraN8eDjsMOneGPfeMrGRJqLSDPISwA/AiUK/yeGPM7Op0jyuSNF9/DVOmwNSpMGMGvPkmbNjgX9t1V2jXzoN3r72gVSu/0t5lF7/qrl8f6tWDggJYt86v0Feu9Kv1zz+HxYvho4/g/ffhvffgjjv86h38eL17+6N7d2g0fiSUlXHMp5968g8fDiUlsf27SOZFcUW+GuhuZhUhhDrAjBDCs2Y2K4Jji2S1//s/GDsWnn4aZs2C9euhUSMP7KuugsMPh4MPhubNUz9mnTr+2HFHD/u2bTd/zurVfpU/a5b/8nj4Ybj7bqhfdx391tfl1PXt6c1i6i1YAKWl/k0K85yVdpCb7xVXUflhncqH9o+TnDV/PjzyCIwe7VfdAIccAsOGwbHHeojXqZPZGurVg06d/HHhhX51PmMGPNH/UUZX9OFxTqIZSzmb+zh35T20KitTkOewEMWenSGE2sA84L+Au8zssiqeUwqUAhQWFhaPGjWqWueqqKigUaNGaVSbPdSW7LOldqxaVYvy8uY899yuvPHGTgC0a/cVRx+9lKOOWkZh4eoarrRqx3TvzjqrzRR68Bd+znj6s4Fa9OMpet/dgjZtvo27xGrJlZ8vSK8t3bp1m2dmHTf7gplF9gB2AsqBA7f2vOLiYquu8vLyan9vtlFbss+m7XjnHbNf/tKscWMzMNtvP7MbbjD75JM4qktBUZEXWvlYyB5WxvXWpNZyA7M+fcxefjnuIrdfrvx8maXXFmCuVZGpkY4jN7OvKoO8d5THFalJZjB5st88bNsW7r0XjjsOpk/3G41lZX6zMisNHw4NGvznwz35Fzc0uJH590zixht9FM3hh8PgwT6CRnJD2kEeQtglhLBT5fv1gZ8C76V7XJGatm4dTJnSnPbtva/79dfh+uth0SJ49FE4+mgIIe4qt6GkBEaMgKIiLAQffD5iBD/5+ckMG+b9+9deC888A/vvD7/9LaxYEXfRkq4orsh3A8pDCG8Ac4DnzWxCBMcVqRFr18J990GbNnDDDW1Zvx4efNBD74orfIhgopSUwPz5TJ861Rux0U3Ohg19NM0HH8Cpp8Itt8CBB8KkSfGVK+lLO8jN7A0z62BmB5nZgWZ2XRSFiWTaunXw0EMe4EOHwk47wXXXvcWbb8KQIT4yJFftsYf/snrxRW9n795w+unpTVaS+GitFck7Zj72+8AD4cwzoXFjmDDB+4+POmoZtfLof8VRR8Frr8GVV8KoUT5j9B//iLsq2V559CMr4mOtu3SBgQO9v/uJJ2DePPjv/05A/3eG7LADXHcdzJzp49+7doWrr/bJTZIMCnLJC/Pn+0iNo47y90eM8Mk8J56YvwG+qU6d4NVXvYvluuugTx9fIkCyn4JcctqKFd5t0KaNd59cfbVPq//5z31dE/mxHXf0vvO//tWHWx56qI/ekeymIJec9H0/+P77ww03eFfK++/DNdf4yA3ZurPP9r7yNWt83Pm4cXFXJFujIJec89FHPoFn4EBo0sQDaeRIaNEi7sqSpVMnv3/Qvr13Qf35z3FXJFuiIJecsWYN/P73Phpl5kz40588iI48Mu7Kkquw0FdX7NfPF+f6zW9+WJpXsod6CSUnvPyy93u//TYMGuQhvvvucVeVGxo08NE9F10Et94Ky5b5BKrateOuTL6nIJdEq6jw2Zd33OE75Tz9NPTtG3dVuad2be9aad7cbxivWgV/+1vml+uV1CjIJbGmTvWbcvPnwy9+ATfe6KMuJDNC8On99ev7Gi2rV/skolyeAZsU6iOXxKmo8ODu0cOvCF98Ee68UyFeUy691K/Ox42Dk0/2tWokXgpySZTp0+Ggg3xbs4sv9unlRx0Vd1X554ILPMzHj/dlDjQLNF7qWpFEWLXK+8Jvuw323tsDXQEerwsu8L+Ofvc7vyE6YoRmycZFQS5Z77XX4LTTfETKuef60qs5sutX4g0b5mE+fDg0bQo33RR3RflJQS5Za8MGH+5WVgbNmsHEib7+h2SX66+HL7+Em2+Gli39/oXULAW5ZKWFC+GMM2DaNJ9VOGIE7Lxz3FVJVULw4Z+LFvmkoRYtfAKR1Bzd7JSsM3asTwufM8cnnowZoxDPdgUF8NhjcMghcMop8M9/xl1RflGQS9ZYuRJKS32NlH328SVVzzpLN9CSomFDX2GysBAGDIDPPou7ovyhIJes8Oab0LGjL5962WW+Vsq++8ZdlWyvwkJ46in45hvvElu1Ku6K8oOCXGJlBvfc4yvtLV8OkyfD//4v1K0bd2VSXe3awcMPw+zZPsrILO6Kcp+CXGLz1Ve+a89558Exx/gGBj17xl2VROHEE31Nloce0vK3NUFBLrGYM8dvjI0b58PWJk70BZkkd1x1FfTvD7/+NcyaFXc1uU1BLjXKDG6/3TdAXr/e10m59FLyauf6fFGrFjzwgK9KefLJPtZcMkP/faTGLF/uf3JffLHv4PPqq76NmOSuJk3g73+Hzz+HIUO0KUWmKMilRnzflTJhAvzxj/Dkkz6lW3Jfx44+Q3fCBPjDH+KuJjcpyCWjzHyJ2S5d/Gpsxgz41a80Njzf/OIXPj+grMy335NoKcglY775xvtGf/lLOPZY70o57LC4q5I4hODLLBQWQkmJT/6S6CjIJSPeeMP/pB471keljB+vrpR817SpD0d8/33fxFmioyCXyD3wgF95V1RAeblGpcgPevTw4Yh33+195hKNtP97hRBahBDKQwjvhBDeDiFcFEVhkjzffed7aJ51FhxxhHelaPMH2dTw4b7L09Ch8MUXcVeTG6K4TloH/NrM2gKdgV+EENpGcFxJgpEjoVUrdu9WSufG73L//b6Tz+TJ3h8qsql69XwK/xdf+I1vSV/aQW5mi83slcr3vwXeBfZI97iSACNHQmkpTywopiNzWbS2ORPrncD1bUZSu3bcxUk2a98eLr8cHnlEXSxRiLTnMoTQCugAzI7yuJKd1l5+NZesvJ5BPMH+vMurdKDP6nE+xkxkG8rKfIGtc87xdXek+oJFtDRZCKERMB0YbmZjq/h6KVAKUFhYWDxq1KhqnaeiooJGObJhY5LbsnRpPW4fvIGX6MIvuYM/8BvqshYAC4HpU6fGXGH1JPk12VQS2vL++404//xievX6nMsue3+Lz0tCW1KVTlu6des2z8w6bvYFM0v7AdQBJgGXpPL84uJiq67y8vJqf2+2SWpbJk82a9bMrFH41h7nJDOf9/PDo6go7hKrLamvSVWS0pZhw/zHZmvlJqUtqUinLcBcqyJToxi1EoD7gHfN7LZ0jyfZa/16uOYan9yz664w96apDG7wzI+f1KCBD0sQSdGVV8Jee/lyxqtXx11NMkXRR94FOB3oHkJ4rfJxXATHlSyydKkvdHXttXD66b5pQOtL+/l0vaIiLAQoKvKPS0riLlcSpEEDuOsueO89rcVSXQXpHsDMZgBaOSOHzZjhG+ouWwZ/+YuPFf/PWiklJVBSwvRp0+jatWucZUqC9ekDJ50EN9zgP2v77BN3Rcmi+XayRWZ+hdS1K9Sv75sDDB2qBa8kM26/HerUgfPP1/Zw20tBLlX68kvf3eXSS31H9Llz4eCD465Kctnuu/sV+eTJvjaPpE5BLpuZPdvXDn/uOfjTn3xjgMaN465K8sF558EBB8All8CqVXFXkxwKcvkPM9/04cgjvftk5ky48EJ1pUjNqVPHLx4++QRu0xi4lCnIBfB1L/r18yuhvn3hlVfg0EPjrkryUY8ecMIJPop10aK4q0kGBbkwcyZ06OB9k3fc4WuIN2kSd1WSz2691ectXHZZ3JUkg4I8j61f71c9xxzjf9K+9JLv5qOuFInbXnv55hOPPur3bGTrFOR5avFin6F5xRUweLCvHV5cHHdVIj+47DJo3txHTmk44tYpyPPQhAm+sP/LL8P99/tqtD/5SdxVifzYjjv6khD/+AfMnLlz3OVkNQV5Hlm1ykehHH887Lmn72b+s5+pK0Wy19Ch0Lo1jBixD2vXxl1N9lKQ54k33vBRKH/+s+/KMmsWtGkTd1UiW1enDtx0Eyxc2IC//jXuarKXgjzHbdjgY8MPPdQXvpo40T+uVy/uykRS068fHHTQV1xzjW/oLZtTkOewhQv9huYll0Dv3vDmm744kUiShADnnPMRS5b4ZCHZnII8B5n5Xojt2vkNzXvvhXHjYJdd4q5MpHratv2W44+HW26B5cvjrib7KMhzzJIlvhzo6afDgQfC669DaaluaEryXX89fP211iyvioI8h4wZ4wsOPf003HgjTJ+udZ0ld7RvDyef7N0rS5bEXU12UZDngCVLfDH+k07yTXpeeQWGDYPateOuTCRa114L333nFyryAwV5gpn5ZJ62beHJJ30t55df9qtykVzUujUMGQJ3360FtTamIE+o+fN9lcLTToN99/Up9mVlPu5WJJdddZWvE3TzzXFXkj0U5Amzbp3f7DngAO8Dv/1231Ozbdu4KxOpGa1awRln+D7fixfHXU12UJAnyEsvQceOvohQjx7wzjtw0UXqC5f8c/nlflFzyy1xV5IdFOQJsGSJr4nSpYtvAPHEE76nYcuWcVcmEo999oGSErjnHvj3v+OuJn4K8iy2Zo1Pp99vP7+pOWwYvPcenHiixoWLXH45rF7tm1DkOwV5FjLzpWbbtfPp9Z07+6JXN94IDRvGXZ1Idmjd2ofd3nWXryOUzxTkWWbOHOje3ZeaDQGeecZ3s9dKhSKbKyuDlSt9Vc98piDPEh984FcXnTrBW2/53plvvgnHHRd3ZSLZq21b6N8f7rwzv1dGVJDH7JNP/Ebm/vv71PorroCPPvK9MzUmXGTbhg3zhbTyeb1yBXlMPvgAbr65NfvtB4895sMIP/7YFwbStmsiqevc2TcQv/VWHyCQjxTkNWzuXF/4p00bmDKlOeed51fgt90GhYVxVyeSTJdd5lP2H3ss7krioSCvAevX+3rgRx/tO/U8+6z/4D322CzuuAP22CPuCkWSrXdv31D8ppt8V6x8E0mQhxDuDyEsCSG8FcXxqjRyJLRqxTHdu/sc3ZEjM3aqqCxe7AtZ7b03nHACfPqp//m3cKEPJWzaVLvJikQhBL84evddH+mVb6K6In8Q6B3RsTY3cqTvjrBgAcEMFizwj7MwzFet8nXB+/f3mZdXXukTesaMgQ8/9HHhjRvHXaVI7hk82P/P3XZb3JXUvEiC3MxeBL6M4lhVqhwsOptOPE1f1lDHB4+WlWXslNtj9Wq/CjjrLNhtN18XfM4cuPhiv6n5/PMwcCAUFMRdqUjuKijw0V7TpvlqoPkkmFk0BwqhFTDBzA7cwtdLgVKAwsLC4lGjRqV87GO6dyeYMYQHeZghNOFLBvIEgxmNTbqCunWjacP2WLq0Hv/8ZxPmzGnK3LlNWbGigIYN13HEEcvo1evfdOiwfJuLWVVUVNCoUaOaKTjDcqUtudIOyM+2VFQUMHhwZ448chmXX/5eDVS2/dJ5Xbp16zbPzDpu9gUzi+QBtALeSuW5xcXFtl2KiszA1lBgz9DHTuNha8i3BmYNG5odf7zZHXeYzZljtnr19h06FWvXmr3zjtn995udfbZZ69ZmPpHebPfdzYYONZs4cfvPXV5eHn2xMcmVtuRKO8zyty0XXmhWUGC2aFHm6klHOq8LMNeqyNRk/LE/fDiUllJn5UqO41mO41lW1t+ZF84dw6TVXXn2WZ9MA7DDDnDwwT7BpnVr75/ebTfYdVdo3hzq1998wakNG3xCwdKl/vj0Ux/T/dFH8PbbPtNy1Sp/btOmcMQRMHSo3yk/4AAtYCWSTS66yKfs33UX/P73cVdTM5IR5CUl/rasDPv0U0LLljQYPpx+JV3ph18bf/opzJ7tj3nzYOJEeOCBqg+3ww5Qr56vZ7xmDazdwuCR3XbzKcDnn+8bv3bs6OO/a2nQpkjW2ntvGDAA7r3Xb6Plw0JzkQR5COExoCvQLISwCLjazO6L4tj/UVICJSVMnzaNrl27bnJ+33S4qMjvXH/v6699pMjnn/uaxUuW+Matq1b5o6AA6tb1UG/SBHbZBZo18zvfRUV+9S4iyXPJJb6P7d/+BueeG3c1mRdJkJvZqVEcJ2qNG0NxcdxViEhN69IFOnTwxbTOOSf3uz/VSSAiOScEuOACv8c1fXrc1WSeglxEctKpp/rghDvvjLuSzFOQi0hOql/fR5eNG+fLYuQyBbmI5KzzzvPhxffeG3clmaUgF5Gc1aqVb5s4YoQvpZGrFOQiktMuuMAn+o0eHXclmaMgF5Gc1rOnz/C+5564K8kcBbmI5LQQfNXrl17y5TZykYJcRHLekCE+iztXb3oqyEUk5zVrBoMG+ZT9lSvjriZ6CnIRyQvnnOPrLz3+eNyVRE9BLiJ54aijfHnrXOxeUZCLSF74/qbn7Nnw+utxVxMtBbmI5I0zzvD9CEaMiLuSaCnIRSRvNG0KJ54Ijz7qexPkCgW5iOSVs8+Gr77yxbRyhYJcRPJK166w115wX7R7mMVKQS4ieaVWLfjZz2DKFPjkk7iriYaCXETyzpAhPorlwQfjriQaCnIRyTstW8JPfwoPPADr18ddTfoU5CKSl84+23cOmjIl7krSpyAXkbzUv78PR7z//rgrSZ+CXETyUr16vkHz+PG+BkuSKchFJG+dcQasWgVjxsRdSXoU5CKStw49FFq3hocfjruS9CjIRSRvheBX5S++mOwx5QpyEclrp53mbx95JN460qEgF5G81rIldOvm3StmcVdTPQpyEcl7Z5wBH34Is2bFXUn1KMhFJO8NHAj16yf3pmckQR5C6B1CeD+E8GEIYVgUxxQRqSk77ggDBsDo0bBmTdzVbL+0gzyEUBu4C+gDtAVODSG0Tfe4IiI1qaQEvvwSnn8+7kq2XxRX5J2AD83sYzNbA4wC+kdwXBGRGtOrF+y8s+8elDTB0rxNG0IYBPQ2s6GVH58OHGZmF2zyvFKgFKCwsLB41KhR1TpfRUUFjRo1SqvmbKG2ZJ9caQeoLdXxxz/uy+TJuzJ27Ezq19+QkXOk05Zu3brNM7OOm36+IO2qUmRmI4ARAB07drSuXbtW6zjTpk2jut+bbdSW7JMr7QC1pTpq14annoLly4+mT5/MnCMTbYmia+VfQIuNPt6z8nMiIonSpQu0aJG87pUognwOsG8IYa8QQl3gFOCpCI4rIlKjatXyFREnTYJly+KuJnVpB7mZrQMuACYB7wKjzeztdI8rIhKH//kfWLcuWSsiRjKO3Mwmmtl+ZraPmQ2P4pgiInE46CBo2zZZ3Sua2SkispEQ4OSTYcYM+OyzuKtJjYJcRGQTgwf7AlpJ6V5RkIuIbKJNG+9iGT067kpSoyAXEanC4MEwcyYsWhR3JdumIBcRqcLgwf7273+Pt45UKMhFRKqw777QoUMyulcU5CIiWzB4sG82sWBB3JVsnYJcRGQLvu9eyfbRKwpyEZEt2Htv6NgRHn887kq2TkEuIrIVgwbBnDnw6adxV7JlCnIRka048UR/++ST8daxNQpyEZGt2HdfaNcOnngi7kq2TEEuIrINAwf62iuffx53JVVTkIuIbMPAgb72yrhxcVdSNQW5iMg2HHCAd7GMHRt3JVVTkIuIbEMIflVeXg5ffhl3NZtTkIuIpGDgQN856Kks3MhSQS4ikoLiYmjZMjtHryjIRURSEIKPKX/+eaioiLuaH1OQi4ikaMAAWL0aJk2Ku5IfU5CLiKSoSxdo2jT7hiEqyEVEUlRQAH37wjPPwNq1cVfzAwW5iMh2GDAAli/3mZ7ZQkEuIrIdevWCHXbIru4VBbmIyHZo2BB69oTx433afjZQkIuIbKcBA3z7t9dfj7sSpyAXEdlOffv6uPLx4+OuxCnIRUS2U2EhHHFE9vSTK8hFRKqhf3947TVYuDDuShTkIiLV0revv33mmXjrgDSDPIRwUgjh7RDChhBCx6iKEhHJdm3awN57w9NPx11J+lfkbwEnAi9GUIuISGKEAMcfD1OmwIoV8daSVpCb2btm9n5UxYiIJEnfvr6I1tSp8dYRLIIR7SGEacBvzGzuVp5TCpQCFBYWFo8aNapa56qoqKBRo0bV+t5so7Zkn1xpB6gtNWHt2sCAAV3o3n0Jv/71Byl9Tzpt6dat2zwz27wb28y2+gBewLtQNn303+g504CO2zrW94/i4mKrrvLy8mp/b7ZRW7JPrrTDTG2pKYMGme2+u9mGDak9P522AHOtikwt2NZvADPrWa1fHSIieaBvXxgzBl59FQ45JJ4aNPxQRCQNxx3nNz4nTIivhnSHH54QQlgEHA48E0LIsn0zREQya5ddoHPnBAe5mT1pZnuaWT0zKzSzY6MqTEQkKfr2hTlz4PPP4zm/ulZERNJ03HH+9rnn4jm/glxEJE3t28Nuu8Gzz8ZzfgW5iEiaQoA+fWDyZFi3rubPryAXEYlAnz7w1Vcwa1bNn1tBLiISgZ49oXbteLpXFOQiIhHYaSfo0gUmTqz5cyvIRUQi0qePbzbx2Wc1e14FuYhIROIahqggFxGJSLt2sMceNd9PriAXEYlICNC7tw9DXLu25s6rIBcRiVCfPvDNNzU7DFFBLiISoR49fBji5Mk1d04FuYhIhHbaCQ47DCbV4FqwCnIRkYgdeyzMnQvLltXM+RTkIiIR69ULzGDKlJo5n4JcRCRihx4KTZrUXPeKglxEJGK1a/vaK5Mm+ZV5pinIRUQyoFcvn6r/zjuZP5eCXEQkA46t3PiyJrpXFOQiIhnQogXsv7+CXEQk0Y49Fl58Eb77LrPnUZCLiGRIr16wahXMmJHZ8yjIRUQy5OijoU4deOGFzJ5HQS4ikiENG8IRRyjIRUQSrWdPePXVzE7XV5CLiGRQz54+Kai8PHPnUJCLiGRQx47wk59ktntFQS4ikkEFBdCtm4JcRCTRevaEjz/2RyYoyEVEMqxnT3+bqavytII8hHBLCOG9EMIbIYQnQwg7RVSXiEjOaN0a9my6ghcunsAx3btDq1YwcmRkx0/3ivx54EAzOwj4APhd+iWJiOSW8OhIen49likrD/dlbRcsgNLSyMI8rSA3s8lmtq7yw1nAnumXJCKSY8rK6Ln+Ob5kZ17jYP/cypVQVhbJ4YNFtOp5COFp4HEze2QLXy8FSgEKCwuLR40aVa3zVFRU0KhRo2rXmU3UluyTK+0AtSWbHNO9O/+25vycv3A119KReQBYCEyfOjXl43Tr1m2emXXc7AtmttUH8ALwVhWP/hs9pwx4kspfDNt6FBcXW3WVl5dX+3uzjdqSfXKlHWZqS1YpKjLzeUE/fhQVbddhgLlWRaYWbOs3gJn13NrXQwhnAn2BHpUnEhGRjQ0f7n3iK1f+8LkGDfzzEUh31Epv4LdAPzNbua3ni4jkpZISGDECioqwEKCoyD8uKYnk8OmOWrkT2BF4PoTwWgjhnghqEhHJPSUlMH++94nPnx9ZiAPb7lrZGjP7r6gKERGR6tHMThGRhFOQi4gknIJcRCThFOQiIgkX2czO7TppCEuBBdX89mZABjdNqlFqS/bJlXaA2pKt0mlLkZntsuknYwnydIQQ5lpVU1QTSG3JPrnSDlBbslUm2qKuFRGRhFOQi4gkXBKDfETcBURIbck+udIOUFuyVeRtSVwfuYiI/FgSr8hFRGQjCnIRkYTL2iAPIfQOIbwfQvgwhDCsiq/XCyE8Xvn12SGEVjGUmZIU2nJmCGFp5QqSr4UQhsZR57aEEO4PISwJIby1ha+HEMIdle18I4RwSE3XmIoU2tE1hPD1Rq/HVTVdY6pCCC1CCOUhhHdCCG+HEC6q4jlJeV1SaUvWvzYhhB1CCP8MIbxe2Y5rq3hOtPlV1W4TcT+A2sBHwN5AXeB1oO0mzzkfuKfy/VPwbeZir72abTkTuDPuWlNoy9HAIcBbW/j6ccCzQAA6A7Pjrrma7egKTIi7zhTbshtwSOX7O+KboG/685WU1yWVtmT9a1P579yo8v06wGyg8ybPiTS/svWKvBPwoZl9bGZrgFFA/02e0x94qPL9MUCPEEKowRpTlUpbEsHMXgS+3MpT+gMPm5sF7BRC2K1mqktdCu1IDDNbbGavVL7/LfAusMcmT0vK65JKW7Je5b9zReWHdSofm44qiTS/sjXI9wAWbvTxIjZ/Qf/zHDNbB3wN7Fwj1W2fVNoCMLDyz94xIYQWNVNa5FJtaxIcXvmn8bMhhAPiLiYVlX+ed8CvADeWuNdlK22BBLw2IYTaIYTXgCXA82a2xdckivzK1iDPN08DrczsIOB5fvhNLfF4BV/Toj3wZ2BcvOVsWwihEfAE8Csz+ybuetKxjbYk4rUxs/VmdjCwJ9AphHBgJs+XrUH+L2Djq9I9Kz9X5XNCCAVAY+CLGqlu+2yzLWb2hZmtrvzwr0BxDdUWtVRet6xnZt98/6exmU0E6oQQmsVc1haFEOrgwTfSzMZW8ZTEvC7bakvSXhsz+wooB3pv8qVI8ytbg3wOsG8IYa8QQl38ZsBTmzznKWBI5fuDgKlWeecgy2yzLZv0V/bD+waT6CngjMpREp2Br81scdxFba8Qwq7f91eGEDrh/0+y8SKByjrvA941s9u28LREvC6ptCUJr00IYZcQwk6V79cHfgq8t8nTIs2vtPbszBQzWxdCuACYhI/6uN/M3g4hXAfMNbOn8Bf8byGED/EbV6fEV/GWpdiWC0MI/YB1eFvOjK3grQghPIaPGmgWQlgEXI3fyMHM7gEm4iMkPgRWAj+Lp9KtS6Edg4DzQgjrgO+AU7L0IgGgC3A68GZlnyzA5UBLSNbrQmptScJrsxvwUAihNv6LZrSZTchkfmmKvohIwmVr14qIiKRIQS4iknAKchGRhFOQi4gknIJcRCThFOQiIgmnIBcRSbj/BzEl/0U7TcVsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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([0,1,2,3])\n",
    "y = np.array([1,2,3,-2])\n",
    "for i in range(0,4):\n",
    "    plt.plot(x[i],y[i],'o',color='r')\n",
    "\n",
    "\n",
    "# Lagrange補間\n",
    "f = interpolate.lagrange(x,y)\n",
    "print(f)\n",
    "x = np.linspace(0,3, 100)\n",
    "y = f(x)\n",
    "plt.plot(x, y, color = 'b')\n",
    "\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton(ニュートン) の差分商公式\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "もう一つ有名なNewton(ニュートン) の内挿公式は，\n",
    "\n",
    "$$\n",
    "\\begin{array}{rc}\n",
    "F (x )&=F (x _{0})+\n",
    "(x -x _{0})f _{1}\\lfloor x_0,x_1\\rfloor+\n",
    "(x -x _{0})(x -x _{1})\n",
    "f _{2}\\lfloor x_0,x_1,x_2\\rfloor + \\\\\n",
    "& \\cdots +  \\prod_{i=0}^{n-1} (x-x_i) \\, \n",
    "f_n \\lfloor x_0,x_1,\\cdots,x_n \\rfloor\n",
    "\\end{array}\n",
    "$$\n",
    "となる．ここで$f_i \\lfloor\\, \\rfloor$ は次のような関数を意味していて，\n",
    "\n",
    "$$\n",
    "\\begin{array}{rl}\n",
    "f _{1}\\lfloor x_0,x_1\\rfloor &=  \\frac{y_1-y_0}{x_1-x_0} \\\\\n",
    "f _{2}\\lfloor x_0,x_1,x_2\\rfloor &=  \\frac{f _{1}\\lfloor x_1,x_2\\rfloor-\n",
    "f _{1}\\lfloor x_0,x_1\\rfloor}{x_2-x_0} \\\\\n",
    "\\vdots & \\\\\n",
    "f _{n}\\lfloor x_0,x_1,\\cdots,x_n\\rfloor &=  \\frac{f_{n-1}\\lfloor x_1,x_2\\cdots,x_{n}\\rfloor-\n",
    "f _{n-1}\\lfloor x_0,x_1,\\cdots,x_{n-1}\\rfloor}{x_n-x_0} \\\\\n",
    "\\end{array}\n",
    "$$\n",
    "差分商と呼ばれる．得られた多項式は，Lagrange の内挿公式で得られたものと当然一致する．Newtonの内挿公式の利点は，新たなデータ点が増えたときに，新たな項を加えるだけで，内挿式が得られる点である．\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1  1  0 -1]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhS0lEQVR4nO3deXhV1bnH8e+CADJYBJE4QVCvgigiBhHFgakIXgQURL1RsUrjUKvV1kqNs1KvWq21+qi0jhVFigiKKCgELCgUcJ6vAwgVCygOAZnf+8cbK0KAQ84+2Wef8/s8z3kynez9Lk74ZWftNQQzQ0REkqtW3AWIiEh6FOQiIgmnIBcRSTgFuYhIwinIRUQSriCOkzZr1sxatWpVre9dsWIFDRs2jLagmKgt2SdX2gFqS7ZKpy3z5s1bZma7bPr5WIK8VatWzJ07t1rfO23aNLp27RptQTFRW7JPrrQD1JZslU5bQggLqvq8ulZERBJOQS4iknAKchGRhFOQi4gknIJcRCThFOQiIgmnIBcRSbhYxpGL5LOvvoL582HpUn988QWsWgVr1vijbl3YYQeoXx923hkKC2HXXaGoyD8nsikFuUiGfPMNzJkDr7/uj7fego8/9iCvrpYtoXVraN8eDjsMOneGPfeMrGRJqLSDPISwA/AiUK/yeGPM7Op0jyuSNF9/DVOmwNSpMGMGvPkmbNjgX9t1V2jXzoN3r72gVSu/0t5lF7/qrl8f6tWDggJYt86v0Feu9Kv1zz+HxYvho4/g/ffhvffgjjv86h38eL17+6N7d2g0fiSUlXHMp5968g8fDiUlsf27SOZFcUW+GuhuZhUhhDrAjBDCs2Y2K4Jji2S1//s/GDsWnn4aZs2C9euhUSMP7KuugsMPh4MPhubNUz9mnTr+2HFHD/u2bTd/zurVfpU/a5b/8nj4Ybj7bqhfdx391tfl1PXt6c1i6i1YAKWl/k0K85yVdpCb7xVXUflhncqH9o+TnDV/PjzyCIwe7VfdAIccAsOGwbHHeojXqZPZGurVg06d/HHhhX51PmMGPNH/UUZX9OFxTqIZSzmb+zh35T20KitTkOewEMWenSGE2sA84L+Au8zssiqeUwqUAhQWFhaPGjWqWueqqKigUaNGaVSbPdSW7LOldqxaVYvy8uY899yuvPHGTgC0a/cVRx+9lKOOWkZh4eoarrRqx3TvzjqrzRR68Bd+znj6s4Fa9OMpet/dgjZtvo27xGrJlZ8vSK8t3bp1m2dmHTf7gplF9gB2AsqBA7f2vOLiYquu8vLyan9vtlFbss+m7XjnHbNf/tKscWMzMNtvP7MbbjD75JM4qktBUZEXWvlYyB5WxvXWpNZyA7M+fcxefjnuIrdfrvx8maXXFmCuVZGpkY4jN7OvKoO8d5THFalJZjB5st88bNsW7r0XjjsOpk/3G41lZX6zMisNHw4NGvznwz35Fzc0uJH590zixht9FM3hh8PgwT6CRnJD2kEeQtglhLBT5fv1gZ8C76V7XJGatm4dTJnSnPbtva/79dfh+uth0SJ49FE4+mgIIe4qt6GkBEaMgKIiLAQffD5iBD/5+ckMG+b9+9deC888A/vvD7/9LaxYEXfRkq4orsh3A8pDCG8Ac4DnzWxCBMcVqRFr18J990GbNnDDDW1Zvx4efNBD74orfIhgopSUwPz5TJ861Rux0U3Ohg19NM0HH8Cpp8Itt8CBB8KkSfGVK+lLO8jN7A0z62BmB5nZgWZ2XRSFiWTaunXw0EMe4EOHwk47wXXXvcWbb8KQIT4yJFftsYf/snrxRW9n795w+unpTVaS+GitFck7Zj72+8AD4cwzoXFjmDDB+4+POmoZtfLof8VRR8Frr8GVV8KoUT5j9B//iLsq2V559CMr4mOtu3SBgQO9v/uJJ2DePPjv/05A/3eG7LADXHcdzJzp49+7doWrr/bJTZIMCnLJC/Pn+0iNo47y90eM8Mk8J56YvwG+qU6d4NVXvYvluuugTx9fIkCyn4JcctqKFd5t0KaNd59cfbVPq//5z31dE/mxHXf0vvO//tWHWx56qI/ekeymIJec9H0/+P77ww03eFfK++/DNdf4yA3ZurPP9r7yNWt83Pm4cXFXJFujIJec89FHPoFn4EBo0sQDaeRIaNEi7sqSpVMnv3/Qvr13Qf35z3FXJFuiIJecsWYN/P73Phpl5kz40588iI48Mu7Kkquw0FdX7NfPF+f6zW9+WJpXsod6CSUnvPyy93u//TYMGuQhvvvucVeVGxo08NE9F10Et94Ky5b5BKrateOuTL6nIJdEq6jw2Zd33OE75Tz9NPTtG3dVuad2be9aad7cbxivWgV/+1vml+uV1CjIJbGmTvWbcvPnwy9+ATfe6KMuJDNC8On99ev7Gi2rV/skolyeAZsU6iOXxKmo8ODu0cOvCF98Ee68UyFeUy691K/Ox42Dk0/2tWokXgpySZTp0+Ggg3xbs4sv9unlRx0Vd1X554ILPMzHj/dlDjQLNF7qWpFEWLXK+8Jvuw323tsDXQEerwsu8L+Ofvc7vyE6YoRmycZFQS5Z77XX4LTTfETKuef60qs5sutX4g0b5mE+fDg0bQo33RR3RflJQS5Za8MGH+5WVgbNmsHEib7+h2SX66+HL7+Em2+Gli39/oXULAW5ZKWFC+GMM2DaNJ9VOGIE7Lxz3FVJVULw4Z+LFvmkoRYtfAKR1Bzd7JSsM3asTwufM8cnnowZoxDPdgUF8NhjcMghcMop8M9/xl1RflGQS9ZYuRJKS32NlH328SVVzzpLN9CSomFDX2GysBAGDIDPPou7ovyhIJes8Oab0LGjL5962WW+Vsq++8ZdlWyvwkJ46in45hvvElu1Ku6K8oOCXGJlBvfc4yvtLV8OkyfD//4v1K0bd2VSXe3awcMPw+zZPsrILO6Kcp+CXGLz1Ve+a89558Exx/gGBj17xl2VROHEE31Nloce0vK3NUFBLrGYM8dvjI0b58PWJk70BZkkd1x1FfTvD7/+NcyaFXc1uU1BLjXKDG6/3TdAXr/e10m59FLyauf6fFGrFjzwgK9KefLJPtZcMkP/faTGLF/uf3JffLHv4PPqq76NmOSuJk3g73+Hzz+HIUO0KUWmKMilRnzflTJhAvzxj/Dkkz6lW3Jfx44+Q3fCBPjDH+KuJjcpyCWjzHyJ2S5d/Gpsxgz41a80Njzf/OIXPj+grMy335NoKcglY775xvtGf/lLOPZY70o57LC4q5I4hODLLBQWQkmJT/6S6CjIJSPeeMP/pB471keljB+vrpR817SpD0d8/33fxFmioyCXyD3wgF95V1RAeblGpcgPevTw4Yh33+195hKNtP97hRBahBDKQwjvhBDeDiFcFEVhkjzffed7aJ51FhxxhHelaPMH2dTw4b7L09Ch8MUXcVeTG6K4TloH/NrM2gKdgV+EENpGcFxJgpEjoVUrdu9WSufG73L//b6Tz+TJ3h8qsql69XwK/xdf+I1vSV/aQW5mi83slcr3vwXeBfZI97iSACNHQmkpTywopiNzWbS2ORPrncD1bUZSu3bcxUk2a98eLr8cHnlEXSxRiLTnMoTQCugAzI7yuJKd1l5+NZesvJ5BPMH+vMurdKDP6nE+xkxkG8rKfIGtc87xdXek+oJFtDRZCKERMB0YbmZjq/h6KVAKUFhYWDxq1KhqnaeiooJGObJhY5LbsnRpPW4fvIGX6MIvuYM/8BvqshYAC4HpU6fGXGH1JPk12VQS2vL++404//xievX6nMsue3+Lz0tCW1KVTlu6des2z8w6bvYFM0v7AdQBJgGXpPL84uJiq67y8vJqf2+2SWpbJk82a9bMrFH41h7nJDOf9/PDo6go7hKrLamvSVWS0pZhw/zHZmvlJqUtqUinLcBcqyJToxi1EoD7gHfN7LZ0jyfZa/16uOYan9yz664w96apDG7wzI+f1KCBD0sQSdGVV8Jee/lyxqtXx11NMkXRR94FOB3oHkJ4rfJxXATHlSyydKkvdHXttXD66b5pQOtL+/l0vaIiLAQoKvKPS0riLlcSpEEDuOsueO89rcVSXQXpHsDMZgBaOSOHzZjhG+ouWwZ/+YuPFf/PWiklJVBSwvRp0+jatWucZUqC9ekDJ50EN9zgP2v77BN3Rcmi+XayRWZ+hdS1K9Sv75sDDB2qBa8kM26/HerUgfPP1/Zw20tBLlX68kvf3eXSS31H9Llz4eCD465Kctnuu/sV+eTJvjaPpE5BLpuZPdvXDn/uOfjTn3xjgMaN465K8sF558EBB8All8CqVXFXkxwKcvkPM9/04cgjvftk5ky48EJ1pUjNqVPHLx4++QRu0xi4lCnIBfB1L/r18yuhvn3hlVfg0EPjrkryUY8ecMIJPop10aK4q0kGBbkwcyZ06OB9k3fc4WuIN2kSd1WSz2691ectXHZZ3JUkg4I8j61f71c9xxzjf9K+9JLv5qOuFInbXnv55hOPPur3bGTrFOR5avFin6F5xRUweLCvHV5cHHdVIj+47DJo3txHTmk44tYpyPPQhAm+sP/LL8P99/tqtD/5SdxVifzYjjv6khD/+AfMnLlz3OVkNQV5Hlm1ykehHH887Lmn72b+s5+pK0Wy19Ch0Lo1jBixD2vXxl1N9lKQ54k33vBRKH/+s+/KMmsWtGkTd1UiW1enDtx0Eyxc2IC//jXuarKXgjzHbdjgY8MPPdQXvpo40T+uVy/uykRS068fHHTQV1xzjW/oLZtTkOewhQv9huYll0Dv3vDmm744kUiShADnnPMRS5b4ZCHZnII8B5n5Xojt2vkNzXvvhXHjYJdd4q5MpHratv2W44+HW26B5cvjrib7KMhzzJIlvhzo6afDgQfC669DaaluaEryXX89fP211iyvioI8h4wZ4wsOPf003HgjTJ+udZ0ld7RvDyef7N0rS5bEXU12UZDngCVLfDH+k07yTXpeeQWGDYPateOuTCRa114L333nFyryAwV5gpn5ZJ62beHJJ30t55df9qtykVzUujUMGQJ3360FtTamIE+o+fN9lcLTToN99/Up9mVlPu5WJJdddZWvE3TzzXFXkj0U5Amzbp3f7DngAO8Dv/1231Ozbdu4KxOpGa1awRln+D7fixfHXU12UJAnyEsvQceOvohQjx7wzjtw0UXqC5f8c/nlflFzyy1xV5IdFOQJsGSJr4nSpYtvAPHEE76nYcuWcVcmEo999oGSErjnHvj3v+OuJn4K8iy2Zo1Pp99vP7+pOWwYvPcenHiixoWLXH45rF7tm1DkOwV5FjLzpWbbtfPp9Z07+6JXN94IDRvGXZ1Idmjd2ofd3nWXryOUzxTkWWbOHOje3ZeaDQGeecZ3s9dKhSKbKyuDlSt9Vc98piDPEh984FcXnTrBW2/53plvvgnHHRd3ZSLZq21b6N8f7rwzv1dGVJDH7JNP/Ebm/vv71PorroCPPvK9MzUmXGTbhg3zhbTyeb1yBXlMPvgAbr65NfvtB4895sMIP/7YFwbStmsiqevc2TcQv/VWHyCQjxTkNWzuXF/4p00bmDKlOeed51fgt90GhYVxVyeSTJdd5lP2H3ss7krioSCvAevX+3rgRx/tO/U8+6z/4D322CzuuAP22CPuCkWSrXdv31D8ppt8V6x8E0mQhxDuDyEsCSG8FcXxqjRyJLRqxTHdu/sc3ZEjM3aqqCxe7AtZ7b03nHACfPqp//m3cKEPJWzaVLvJikQhBL84evddH+mVb6K6In8Q6B3RsTY3cqTvjrBgAcEMFizwj7MwzFet8nXB+/f3mZdXXukTesaMgQ8/9HHhjRvHXaVI7hk82P/P3XZb3JXUvEiC3MxeBL6M4lhVqhwsOptOPE1f1lDHB4+WlWXslNtj9Wq/CjjrLNhtN18XfM4cuPhiv6n5/PMwcCAUFMRdqUjuKijw0V7TpvlqoPkkmFk0BwqhFTDBzA7cwtdLgVKAwsLC4lGjRqV87GO6dyeYMYQHeZghNOFLBvIEgxmNTbqCunWjacP2WLq0Hv/8ZxPmzGnK3LlNWbGigIYN13HEEcvo1evfdOiwfJuLWVVUVNCoUaOaKTjDcqUtudIOyM+2VFQUMHhwZ448chmXX/5eDVS2/dJ5Xbp16zbPzDpu9gUzi+QBtALeSuW5xcXFtl2KiszA1lBgz9DHTuNha8i3BmYNG5odf7zZHXeYzZljtnr19h06FWvXmr3zjtn995udfbZZ69ZmPpHebPfdzYYONZs4cfvPXV5eHn2xMcmVtuRKO8zyty0XXmhWUGC2aFHm6klHOq8LMNeqyNRk/LE/fDiUllJn5UqO41mO41lW1t+ZF84dw6TVXXn2WZ9MA7DDDnDwwT7BpnVr75/ebTfYdVdo3hzq1998wakNG3xCwdKl/vj0Ux/T/dFH8PbbPtNy1Sp/btOmcMQRMHSo3yk/4AAtYCWSTS66yKfs33UX/P73cVdTM5IR5CUl/rasDPv0U0LLljQYPpx+JV3ph18bf/opzJ7tj3nzYOJEeOCBqg+3ww5Qr56vZ7xmDazdwuCR3XbzKcDnn+8bv3bs6OO/a2nQpkjW2ntvGDAA7r3Xb6Plw0JzkQR5COExoCvQLISwCLjazO6L4tj/UVICJSVMnzaNrl27bnJ+33S4qMjvXH/v6699pMjnn/uaxUuW+Matq1b5o6AA6tb1UG/SBHbZBZo18zvfRUV+9S4iyXPJJb6P7d/+BueeG3c1mRdJkJvZqVEcJ2qNG0NxcdxViEhN69IFOnTwxbTOOSf3uz/VSSAiOScEuOACv8c1fXrc1WSeglxEctKpp/rghDvvjLuSzFOQi0hOql/fR5eNG+fLYuQyBbmI5KzzzvPhxffeG3clmaUgF5Gc1aqVb5s4YoQvpZGrFOQiktMuuMAn+o0eHXclmaMgF5Gc1rOnz/C+5564K8kcBbmI5LQQfNXrl17y5TZykYJcRHLekCE+iztXb3oqyEUk5zVrBoMG+ZT9lSvjriZ6CnIRyQvnnOPrLz3+eNyVRE9BLiJ54aijfHnrXOxeUZCLSF74/qbn7Nnw+utxVxMtBbmI5I0zzvD9CEaMiLuSaCnIRSRvNG0KJ54Ijz7qexPkCgW5iOSVs8+Gr77yxbRyhYJcRPJK166w115wX7R7mMVKQS4ieaVWLfjZz2DKFPjkk7iriYaCXETyzpAhPorlwQfjriQaCnIRyTstW8JPfwoPPADr18ddTfoU5CKSl84+23cOmjIl7krSpyAXkbzUv78PR7z//rgrSZ+CXETyUr16vkHz+PG+BkuSKchFJG+dcQasWgVjxsRdSXoU5CKStw49FFq3hocfjruS9CjIRSRvheBX5S++mOwx5QpyEclrp53mbx95JN460qEgF5G81rIldOvm3StmcVdTPQpyEcl7Z5wBH34Is2bFXUn1KMhFJO8NHAj16yf3pmckQR5C6B1CeD+E8GEIYVgUxxQRqSk77ggDBsDo0bBmTdzVbL+0gzyEUBu4C+gDtAVODSG0Tfe4IiI1qaQEvvwSnn8+7kq2XxRX5J2AD83sYzNbA4wC+kdwXBGRGtOrF+y8s+8elDTB0rxNG0IYBPQ2s6GVH58OHGZmF2zyvFKgFKCwsLB41KhR1TpfRUUFjRo1SqvmbKG2ZJ9caQeoLdXxxz/uy+TJuzJ27Ezq19+QkXOk05Zu3brNM7OOm36+IO2qUmRmI4ARAB07drSuXbtW6zjTpk2jut+bbdSW7JMr7QC1pTpq14annoLly4+mT5/MnCMTbYmia+VfQIuNPt6z8nMiIonSpQu0aJG87pUognwOsG8IYa8QQl3gFOCpCI4rIlKjatXyFREnTYJly+KuJnVpB7mZrQMuACYB7wKjzeztdI8rIhKH//kfWLcuWSsiRjKO3Mwmmtl+ZraPmQ2P4pgiInE46CBo2zZZ3Sua2SkispEQ4OSTYcYM+OyzuKtJjYJcRGQTgwf7AlpJ6V5RkIuIbKJNG+9iGT067kpSoyAXEanC4MEwcyYsWhR3JdumIBcRqcLgwf7273+Pt45UKMhFRKqw777QoUMyulcU5CIiWzB4sG82sWBB3JVsnYJcRGQLvu9eyfbRKwpyEZEt2Htv6NgRHn887kq2TkEuIrIVgwbBnDnw6adxV7JlCnIRka048UR/++ST8daxNQpyEZGt2HdfaNcOnngi7kq2TEEuIrINAwf62iuffx53JVVTkIuIbMPAgb72yrhxcVdSNQW5iMg2HHCAd7GMHRt3JVVTkIuIbEMIflVeXg5ffhl3NZtTkIuIpGDgQN856Kks3MhSQS4ikoLiYmjZMjtHryjIRURSEIKPKX/+eaioiLuaH1OQi4ikaMAAWL0aJk2Ku5IfU5CLiKSoSxdo2jT7hiEqyEVEUlRQAH37wjPPwNq1cVfzAwW5iMh2GDAAli/3mZ7ZQkEuIrIdevWCHXbIru4VBbmIyHZo2BB69oTx433afjZQkIuIbKcBA3z7t9dfj7sSpyAXEdlOffv6uPLx4+OuxCnIRUS2U2EhHHFE9vSTK8hFRKqhf3947TVYuDDuShTkIiLV0revv33mmXjrgDSDPIRwUgjh7RDChhBCx6iKEhHJdm3awN57w9NPx11J+lfkbwEnAi9GUIuISGKEAMcfD1OmwIoV8daSVpCb2btm9n5UxYiIJEnfvr6I1tSp8dYRLIIR7SGEacBvzGzuVp5TCpQCFBYWFo8aNapa56qoqKBRo0bV+t5so7Zkn1xpB6gtNWHt2sCAAV3o3n0Jv/71Byl9Tzpt6dat2zwz27wb28y2+gBewLtQNn303+g504CO2zrW94/i4mKrrvLy8mp/b7ZRW7JPrrTDTG2pKYMGme2+u9mGDak9P522AHOtikwt2NZvADPrWa1fHSIieaBvXxgzBl59FQ45JJ4aNPxQRCQNxx3nNz4nTIivhnSHH54QQlgEHA48E0LIsn0zREQya5ddoHPnBAe5mT1pZnuaWT0zKzSzY6MqTEQkKfr2hTlz4PPP4zm/ulZERNJ03HH+9rnn4jm/glxEJE3t28Nuu8Gzz8ZzfgW5iEiaQoA+fWDyZFi3rubPryAXEYlAnz7w1Vcwa1bNn1tBLiISgZ49oXbteLpXFOQiIhHYaSfo0gUmTqz5cyvIRUQi0qePbzbx2Wc1e14FuYhIROIahqggFxGJSLt2sMceNd9PriAXEYlICNC7tw9DXLu25s6rIBcRiVCfPvDNNzU7DFFBLiISoR49fBji5Mk1d04FuYhIhHbaCQ47DCbV4FqwCnIRkYgdeyzMnQvLltXM+RTkIiIR69ULzGDKlJo5n4JcRCRihx4KTZrUXPeKglxEJGK1a/vaK5Mm+ZV5pinIRUQyoFcvn6r/zjuZP5eCXEQkA46t3PiyJrpXFOQiIhnQogXsv7+CXEQk0Y49Fl58Eb77LrPnUZCLiGRIr16wahXMmJHZ8yjIRUQy5OijoU4deOGFzJ5HQS4ikiENG8IRRyjIRUQSrWdPePXVzE7XV5CLiGRQz54+Kai8PHPnUJCLiGRQx47wk59ktntFQS4ikkEFBdCtm4JcRCTRevaEjz/2RyYoyEVEMqxnT3+bqavytII8hHBLCOG9EMIbIYQnQwg7RVSXiEjOaN0a9my6ghcunsAx3btDq1YwcmRkx0/3ivx54EAzOwj4APhd+iWJiOSW8OhIen49likrD/dlbRcsgNLSyMI8rSA3s8lmtq7yw1nAnumXJCKSY8rK6Ln+Ob5kZ17jYP/cypVQVhbJ4YNFtOp5COFp4HEze2QLXy8FSgEKCwuLR40aVa3zVFRU0KhRo2rXmU3UluyTK+0AtSWbHNO9O/+25vycv3A119KReQBYCEyfOjXl43Tr1m2emXXc7AtmttUH8ALwVhWP/hs9pwx4kspfDNt6FBcXW3WVl5dX+3uzjdqSfXKlHWZqS1YpKjLzeUE/fhQVbddhgLlWRaYWbOs3gJn13NrXQwhnAn2BHpUnEhGRjQ0f7n3iK1f+8LkGDfzzEUh31Epv4LdAPzNbua3ni4jkpZISGDECioqwEKCoyD8uKYnk8OmOWrkT2BF4PoTwWgjhnghqEhHJPSUlMH++94nPnx9ZiAPb7lrZGjP7r6gKERGR6tHMThGRhFOQi4gknIJcRCThFOQiIgkX2czO7TppCEuBBdX89mZABjdNqlFqS/bJlXaA2pKt0mlLkZntsuknYwnydIQQ5lpVU1QTSG3JPrnSDlBbslUm2qKuFRGRhFOQi4gkXBKDfETcBURIbck+udIOUFuyVeRtSVwfuYiI/FgSr8hFRGQjCnIRkYTL2iAPIfQOIbwfQvgwhDCsiq/XCyE8Xvn12SGEVjGUmZIU2nJmCGFp5QqSr4UQhsZR57aEEO4PISwJIby1ha+HEMIdle18I4RwSE3XmIoU2tE1hPD1Rq/HVTVdY6pCCC1CCOUhhHdCCG+HEC6q4jlJeV1SaUvWvzYhhB1CCP8MIbxe2Y5rq3hOtPlV1W4TcT+A2sBHwN5AXeB1oO0mzzkfuKfy/VPwbeZir72abTkTuDPuWlNoy9HAIcBbW/j6ccCzQAA6A7Pjrrma7egKTIi7zhTbshtwSOX7O+KboG/685WU1yWVtmT9a1P579yo8v06wGyg8ybPiTS/svWKvBPwoZl9bGZrgFFA/02e0x94qPL9MUCPEEKowRpTlUpbEsHMXgS+3MpT+gMPm5sF7BRC2K1mqktdCu1IDDNbbGavVL7/LfAusMcmT0vK65JKW7Je5b9zReWHdSofm44qiTS/sjXI9wAWbvTxIjZ/Qf/zHDNbB3wN7Fwj1W2fVNoCMLDyz94xIYQWNVNa5FJtaxIcXvmn8bMhhAPiLiYVlX+ed8CvADeWuNdlK22BBLw2IYTaIYTXgCXA82a2xdckivzK1iDPN08DrczsIOB5fvhNLfF4BV/Toj3wZ2BcvOVsWwihEfAE8Csz+ybuetKxjbYk4rUxs/VmdjCwJ9AphHBgJs+XrUH+L2Djq9I9Kz9X5XNCCAVAY+CLGqlu+2yzLWb2hZmtrvzwr0BxDdUWtVRet6xnZt98/6exmU0E6oQQmsVc1haFEOrgwTfSzMZW8ZTEvC7bakvSXhsz+wooB3pv8qVI8ytbg3wOsG8IYa8QQl38ZsBTmzznKWBI5fuDgKlWeecgy2yzLZv0V/bD+waT6CngjMpREp2Br81scdxFba8Qwq7f91eGEDrh/0+y8SKByjrvA941s9u28LREvC6ptCUJr00IYZcQwk6V79cHfgq8t8nTIs2vtPbszBQzWxdCuACYhI/6uN/M3g4hXAfMNbOn8Bf8byGED/EbV6fEV/GWpdiWC0MI/YB1eFvOjK3grQghPIaPGmgWQlgEXI3fyMHM7gEm4iMkPgRWAj+Lp9KtS6Edg4DzQgjrgO+AU7L0IgGgC3A68GZlnyzA5UBLSNbrQmptScJrsxvwUAihNv6LZrSZTchkfmmKvohIwmVr14qIiKRIQS4iknAKchGRhFOQi4gknIJcRCThFOQiIgmnIBcRSbj/BzEl/0U7TcVsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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([0,1,2,3])\n",
    "y = np.array([1,2,3,-2])\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(0,3, 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補間と多項式補間の一致の検証 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "関数$F(x)$を$x$の多項式として展開．その時の，係数の取るべき値と，差分商で得られる値が一致．\n",
    "```maple\n",
    "> restart: F:=x->f0+(x-x0)*f1p+(x-x0)*(x-x1)*f2p;\n",
    "```\n",
    "$$\n",
    "F\\, := \\,x\\mapsto f0+ ( x-x0 ) f1p+ ( x-x0 )  ( x-x1 ) f2p\n",
    "$$\n",
    "```maple\n",
    "> F(x1); \n",
    "  sf1p:=solve(F(x1)=f1,f1p);\n",
    "```\n",
    "$$\n",
    "f0+ ( x1-x0 ) f1p \\notag \\\\\n",
    "sf1p\\, := \\,{\\frac {f0-f1}{-x1+x0}} \\notag\n",
    "$$\n",
    "f20の取るべき値の導出\n",
    "```maple\n",
    "> sf2p:=solve(F(x2)=f2,f2p); \n",
    "  fac_f2p:=factor(subs(f1p=sf1p,sf2p));\n",
    "```\n",
    "$$\n",
    "sf2p\\, :=  -{\\frac {f0+f1p\\,x2-f1p\\,x0-f2}{ ( -x2+x0 ) \n",
    "( -x2+x1 ) }} \\notag \\\\\n",
    "{\\it fac\\_f2p}\\, :=  {\\frac {f0\\,x1-x2\\,f0+x2\\,f1-x0\\,f1-f2\\,x1+f2\\,x0}{ ( -x1+x0 )  ( -x2+x0 )  ( -x2+x1 ) }} \\notag\n",
    "$$\n",
    "ニュートンの差分商公式を変形\n",
    "```maple\n",
    "> ff11:=(f0-f1)/(x0-x1); \n",
    "  ff12:=(f1-f2)/(x1-x2); \n",
    "  ff2:=(ff11-ff12)/(x0-x2);\n",
    "  fac_newton:=factor(ff2);\n",
    "```\n",
    "$$\n",
    "ff11:= { \\frac {f0-f1}{-x1+x0}} \\notag \\\\\n",
    "ff12 := { \\frac {f1-f2}{-x2+x1}} \\notag \\\\\n",
    "ff2 :=  \\frac{ { \\frac {f0-f1}{-x1+x0}}-{ \\frac {f1-f2}{-x2+x1}} }{-x2+x0 }\\notag \\\\\n",
    "{\\it fac\\_newton} := { \\frac {f0\\,x1-x2\\,f0+x2\\,f1-x0\\,f1-f2\\,x1+f2\\,x0}{ ( -x1+x0 )  ( -x2+x0 )  ( -x2+x1 ) }} \\notag \n",
    "$$\n",
    "\n",
    "二式が等しいかどうかをevalbで判定\n",
    "```maple\n",
    "> evalb(fac_f2p=fac_newton);\n",
    "```\n",
    "$$\n",
    "true\n",
    "$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数値積分 (Numerical integration)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAN00lEQVR4nO3dYYgc93nH8d+vkoPjWxe/sLMNlnIXaAjIBjvs4bq4tLeiCefYJDQ02EYxGFLuTQwOJFQOfpUXoX0V+sYtFYnRi6hZlyai4ZzKVfFeTSDUvnNsV7JsMK6cWAQUk4Z4JUhQ8vTFraktnX17M39l5hl/P7Dc7uzMf56HkX4azf131hEhAEBev9d0AQCAeghyAEiOIAeA5AhyAEiOIAeA5HY3sdNrr702FhYWKm177tw5zc3NlS2oIfTSPl3pQ6KXtqrTy8bGxusRcd3FyxsJ8oWFBa2vr1fadm1tTUtLS2ULagi9tE9X+pDopa3q9GL71a2Wc2kFAJIjyAEgOYIcAJIjyAEgOYIcAJIrMmvF9mlJb0j6jaQLEbFYYlwAwPZKTj8cRsTrBccDAMyASysAkJxL3I/c9v9I+l9JIekfI+LQFuusSFqRpH6/PxiNRpX2NZlM1Ov1alTbHo+enOiuG7rRS1eOS1f6kOilrer0MhwON7a8dB0RtR+Srp/+/ICk5yT96butPxgMoqrxeFx527aZP7jadAnFdOW4dKWPCHppqzq9SFqPLTK1yKWViDgz/XlW0lFJt5QYFwCwvdpBbnvO9tVvPpf0CUkn6o4LAJhNiVkrfUlHbb853j9FxLEC4wIAZlA7yCPiFUk3FagFAFAB0w8BIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSI8gBIDmCHACSKxbktnfZ/pHt1VJjAgC2V/KM/AFJpwqOBwCYQZEgt71H0h2SvlFiPADA7BwR9Qex/0XS30i6WtKXI+LOLdZZkbQiSf1+fzAajSrtazKZqNfr1ai2Pe47dk6Hl+eaLqOIrhyXrvQh0Utb1ellOBxuRMTiJW9ERK2HpDsl/f30+ZKk1e22GQwGUdV4PK68bdvMH1xtuoRiunJcutJHBL20VZ1eJK3HFpla4tLKbZI+Zfu0pJGk/ba/VWBcAMAMagd5RHwlIvZExIKkuyU9ERGfq10ZAGAmzCMHuubIEWlhQX+2f7+0sLD5Gp22u+RgEbEmaa3kmAB24MgRaWVFOn9elqRXX918LUkHDjRZGS4jzsiBLnnoIen8+bcvO39+czk6iyAHuuTHP97ZcnQCQQ50yYc+tLPl6ASCHOiSr31Nuuqqty+76qrN5egsghzokgMHpEOHpPl5hS3Nz2++5hednUaQA11z4IB0+rT+84knpNOnCfH3AIIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAGJG00htaI3zQJS4kZTSI4zcoAbTSE5ghzgRlNIjiAHuNEUkiPIAW40heQIcoAbTSE5ghyQuNFUWzEtdCZMPwTQTkwLnRln5ADaiWmhMyPIAbQT00JnRpADaCemhc6MIAfQTkwLnVntILd9pe2nbD9n+6Ttr5YoDMB7HNNCZ1bijPxXkvZHxE2Sbpa0bPvWAuMCeK9jWuhMak8/jIiQNJm+vGL6iLrjAgBm480crjmIvUvShqQ/lPRwRBzcYp0VSSuS1O/3B6PRqNK+JpOJer1ejWrb475j53R4ea7pMoroynHpSh8SvbRVnV6Gw+FGRCxe8kZEFHtIukbSWNKN77beYDCIqsbjceVt22b+4GrTJRTTlePSlT4i6KWt6vQiaT22yNSis1Yi4hfTIF8uOS4A4J2VmLVyne1rps/fL+njkl6sOy4AYDYlzsg/KGls+3lJT0s6HhGrBcZFBtzUCGhciVkrz0v6WIFakA03NQJagU92ojpuagS0AkGO6ripEdAKBDmq46ZGQCsQ5KiOmxoBrUCQozpuagS0AkGOeripEdA4ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASA5ghwAkiPIASC52kFue6/tse0XbJ+0/UCJwgAAs9ldYIwLkr4UEc/YvlrShu3jEfFCgbEBANuofUYeET+NiGemz9+QdErS9XXHBQDMxhFRbjB7QdKTkm6MiF9e9N6KpBVJ6vf7g9FoVGkfk8lEvV6vZqXtcN+xczq8PNd0GUV05bh0pQ+JXtqqTi/D4XAjIhYveSMiijwk9SRtSPrMdusOBoOoajweV962beYPrjZdQjFdOS5d6SOCXtqqTi+S1mOLTC0ya8X2FZK+I+lIRHy3xJgAgNmUmLViSd+UdCoivl6/JADATpQ4I79N0r2S9tt+dvr4ZIFxAQAzqD39MCJ+IMkFagEAVMAnOwEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIjyAEgOYIcAJIrEuS2H7F91vaJEuMBAGZX6oz8sKTlQmMBAHagSJBHxJOSfl5iLADAzjgiygxkL0hajYgb3+H9FUkrktTv9wej0ajSfiaTiXq9XtUyW+W+Y+d0eHmu6TKK6Mpx6UofEr20VZ1ehsPhRkQsXvJGRBR5SFqQdGKWdQeDQVQ1Ho8rb9s28wdXmy6hmK4cl670EUEvbVWnF0nrsUWmMmsFAJIjyAEguVLTD78t6YeSPmr7NdufLzEuAGB7u0sMEhH3lBgHALBzXFoBgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIjiAHgOQIcgBIrkiQ2162/ZLtl20/WGJMAMBsage57V2SHpZ0u6R9ku6xva/uuACA2ZQ4I79F0ssR8UpE/FrSSNKnC4wLAJiBI6LeAPZfSlqOiL+avr5X0h9FxP0XrbciaUWS+v3+YDQaVdrfoycn+refuFbNANCU2/eG7rqhV2nb4XC4ERGLFy/fXbuqGUXEIUmHJGlxcTGWlpYqjrSmf/hC1W3bZeHBx3T6b+9ouowi1tbWVP2YtkdX+pDopa0uRy8lLq2ckbT3La/3TJcBAH4HSgT505I+YvvDtt8n6W5J3yswLgBgBrUvrUTEBdv3S3pc0i5Jj0TEydqVAQBmUuQaeUR8X9L3S4wFANgZPtkJAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQHEEOAMkR5ACQXK0gt/1Z2ydt/9b2YqmiAACzq3tGfkLSZyQ9WaAWAEAFu+tsHBGnJMl2mWoAADvmiKg/iL0m6csRsf4u66xIWpGkfr8/GI1GlfY1mUzU6/Uqbds2j56c6K4butFLV45LV/qQ6KWt6vQyHA43IuKSy9jbnpHb/g9Jf7DFWw9FxL/OWkBEHJJ0SJIWFxdjaWlp1k3fZm1tTVW3bZ/u9NKV49KVPiR6aavL0cu2QR4Rf150jwCAoph+CADJ1Z1++Be2X5P0x5Ies/14mbIAALOqO2vlqKSjhWoBAFTApRUASI4gB4DkCHIASI4gB4Dkinyyc8c7tX8m6dWKm18r6fWC5TSJXtqnK31I9NJWdXqZj4jrLl7YSJDXYXt9q4+oZkQv7dOVPiR6aavL0QuXVgAgOYIcAJLLGOSHmi6gIHppn670IdFLWxXvJd01cgDA22U8IwcAvAVBDgDJpQzy7F/6bHvZ9ku2X7b9YNP11GH7EdtnbZ9oupY6bO+1Pbb9wvTP1gNN11SV7SttP2X7uWkvX226pjps77L9I9urTddSh+3Ttv/b9rO23/Hb1KpIGeRK/KXPtndJeljS7ZL2SbrH9r5mq6rlsKTlposo4IKkL0XEPkm3SvpC4uPyK0n7I+ImSTdLWrZ9a7Ml1fKApFNNF1HIMCJuZh65Nr/0OSJearqOim6R9HJEvBIRv5Y0kvTphmuqLCKelPTzpuuoKyJ+GhHPTJ+/oc3guL7ZqqqJTZPpyyumj5SzGmzvkXSHpG80XUubpQzy5K6X9JO3vH5NSQOjq2wvSPqYpP9quJTKppcjnpV0VtLxiMjay99J+mtJv224jhJC0r/b3ph+GX0xtb5Y4nIq9aXPwE7Y7kn6jqQvRsQvm66nqoj4jaSbbV8j6ajtGyMi1e8xbN8p6WxEbNhearicEv4kIs7Y/oCk47ZfnP6PtrbWBnmHv/T5jKS9b3m9Z7oMDbN9hTZD/EhEfLfpekqIiF/YHmvz9xipglzSbZI+ZfuTkq6U9Pu2vxURn2u4rkoi4sz051nbR7V5mbVIkHNp5XfvaUkfsf1h2++TdLek7zVc03uebUv6pqRTEfH1puupw/Z10zNx2X6/pI9LerHRoiqIiK9ExJ6IWNDm35Mnsoa47TnbV7/5XNInVPAf1pRBnvlLnyPigqT7JT2uzV+o/XNEnGy2qupsf1vSDyV91PZrtj/fdE0V3SbpXkn7p9PDnp2eCWb0QUlj289r88TheESknrrXAX1JP7D9nKSnJD0WEcdKDc5H9AEguZRn5ACA/0eQA0ByBDkAJEeQA0ByBDkAJEeQA0ByBDkAJPd/CMiRuDfdJjYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "x = np.array([1,2,3,4])\n",
    "y = np.array([2.8,3.4,3.8,3.2])\n",
    "for i in range(0,4):\n",
    "    plt.plot(x[i],y[i],'o',color='r')\n",
    "plt.hlines(0, -1, 5, linewidth=1)\n",
    "plt.vlines(0, -1, 5, linewidth=1)\n",
    "plt.grid(True)\n",
    "\n",
    "# 数値積分の模式図"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "積分，\n",
    "\n",
    "$$\n",
    "I = \\int_a^b f(x) dx\n",
    "$$\n",
    "を求めよう．1次元の数値積分法では連続した領域を細かい短冊に分けて，それぞれの面積を寄せ集めることに相当する．分点の数を N とすると，\n",
    "\n",
    "$$\n",
    "\\begin{array}{c}\n",
    "x_i = a+ \\frac{b-a}{N} i = a + h \\times i \\\\\n",
    "h = \\frac{b-a}{N}\n",
    "\\end{array}\n",
    "$$\n",
    "ととれる．そうすると，もっとも単純には，\n",
    "\n",
    "$$\n",
    "I_N = \\left\\{\\sum_{i=0}^{N-1} f(x_i)\\right\\}h =\n",
    "\\left\\{\\sum_{i=0}^{N-1} f(a+i \\times h)\\right\\}h\n",
    "$$\n",
    "となる．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 中点則 (midpoint rule) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "中点法 (midpoint rule) は，短冊を左端から書くのではなく，真ん中から書くことに対応し，\n",
    "\n",
    "$$\n",
    "I_N = \\left\\{\\sum_{i=0}^{N-1}f\\left(a+\\left(i+\\frac{1}{2}\\right) \\times h\\right)\\right\\}h\n",
    "$$\n",
    "となる．\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 台形則 (trapezoidal rule) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さらに短冊の上側を斜めにして，短冊を台形にすれば精度が上がりそうに思う．\n",
    "その場合は，短冊一枚の面積$S_i$は，\n",
    "\n",
    "$$\n",
    "S_i = \\frac{f(x_i)+f(x_{i+1})}{2}h\n",
    "$$\n",
    "で求まる．これを端から端まで加えあわせると，\n",
    "\n",
    "$$\n",
    "i_N =\\sum _{i=0}^{N-1}S_i =h \\left\\{ \\frac{1}{2} f ( x_0 ) +\\sum _{i=1}^{N-1}f ( x_i ) +\\frac{1}{2} f \\left( x_N \\right)  \\right\\} \n",
    "$$\n",
    "が得られる．\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simpson(シンプソン)則 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simpson(シンプソン) 則では，短冊を2次関数，\n",
    "\n",
    "$$\n",
    "f(x) = ax^2+bx+c\n",
    "$$\n",
    "で近似することに対応する．こうすると，\n",
    "\n",
    "$$\n",
    "S_i=\\int _{x_i}^{x_{i+1}}f ( x )\\, {dx}=\\int _{x_i}^{x_{i+1}}(ax^2+bx+c)\\,{dx}\n",
    "$$\n",
    "\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Simpson則の導出（数式変形）：\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "\n",
    "\n",
    "$$\n",
    "\\frac{h}{6} \\left\\{f(x_i)+4f\\left(x_i+\\frac{h}{2}\\right)+f(x_i+h)\\right\\}\n",
    "$$\n",
    "となる．これより，\n",
    "\n",
    "$$\n",
    "I_N=\\frac{h}{6} \\left\\{ f \\left( x_0 \\right) +4\\,\\sum _{i=0}^{N-1}f \\left( x_i+\\frac{h}{2} \\right) +2\\,\\sum_{i=1}^{N-1}f \\left( x_i \\right) +f \\left( x_N \\right)  \\right\\}\n",
    "$$\n",
    "として計算できる．ただし，関数値を計算する点の数は台形則などの倍となっている．\n",
    "\n",
    "教科書によっては，分割数$N$を偶数にして，点を偶数番目 (even) と奇数番目 (odd) に分けて，\n",
    "\n",
    "$$\n",
    "I_{{N}}=\\frac{h}{3} \\left\\{ f \\left( x_{{0}} \\right) +4\\,\\sum _{i={\\it even}}^{N-2}f \\left( x_{{i}}+\\frac{h}{2} \\right) +2\\,\\sum _{i={\\it odd}}^{N-1}f \\left( x_{{i}} \\right) +f \\left( x_{{N}} \\right) \\right\\}\n",
    "$$\n",
    "としている記述があるが，同じ計算になるので誤解せぬよう．\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数値積分のコード\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次の積分を例に，pythonのコードを示す．\n",
    "\n",
    "$$\n",
    "\\int_0^1 \\frac{4}{1+x^2} \\, dx\n",
    "$$\n",
    "先ずは問題が与えられたらできるだけお任せで解いてしまう．答えをあらかじめ知っておくと間違いを見つけるのが容易．プロットしてみる．\n",
    "\n",
    "scipyで積分計算をお任せでしてくれる関数はquadで，これはFortran libraryのQUADPACKを利用している．自動で色々してくれるが，精度は1.49e_08以下．"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiiElEQVR4nO3de3hU1bnH8e8LREWiUkQDQjAKFAVUcBBBEAmKIrXoabHF9ihQhWLrUat4QT1UsRy1rbbeqkXkVG1r8FqRg1ovQav1RkDuooC2Sql3xUBFgff8sSYlxoRMJjPZMzu/z/Psx5nZa4Z3nd3zY7P23muZuyMiIvmvRdQFiIhIZijQRURiQoEuIhITCnQRkZhQoIuIxESrqP7g9u3be0lJSVrf3bhxI23atMlsQRFRX3JTXPoSl36A+lKloqLifXffq7Z9kQV6SUkJCxYsSOu78+fPZ+jQoZktKCLqS26KS1/i0g9QX6qY2d/q2qchFxGRmFCgi4jEhAJdRCQmFOgiIjGhQBcRiYmUA93MWprZIjObW8u+nc1stpmtNrMXzawko1WKiEi9GnKGfg6wso59pwMfuXs34FfANY0tTEREGial+9DNrDPwDWA6cF4tTU4ELk++vg+4yczMszE37/LllMyaBYsWwde+Bu3bw0EHQZcuYJbxP05EJF9YKplrZvcBVwG7AZPd/YQa+5cBI9z97eT7NcDh7v5+jXYTgYkARUVFibKysgYXvFd5OT2vvBKrUfcXu+/Op1//Oh/27897Rx7J5g4dGvzbUaisrKSwsDDqMjJCfck9cekHqC9VSktLK9y9X6073X2HG3AC8Jvk66HA3FraLAM6V3u/Bmi/o99NJBKervInn3T/8EP3NWvc//pX91tucZ8wwb13b3cIW79+7r/9rfumTWn/OU2hvLw86hIyRn3JPXHph7v6UgVY4HXkaipj6IOAUWb2JlAGDDOz39dosw4oBjCzVsAewAep/53TQC1ahOGW/feHgQNh0iSYMQOWLoXXX4drroEvvoAf/hD23Rcuvxw+/DBr5YiI5IJ6A93dp7h7Z3cvAcYAT7n7f9ZoNgcYm3w9OtkmmrXtunWDCy8MY+xPPQX9+8MVV0D37nDjjSHoRURiKO370M1smpmNSr69HdjTzFYTLppenIniGsUMSkth7lxYvBj69oWzz4aDDw5BLyISMw0KdHef78kLou4+1d3nJF9/5u4nu3s3d+/v7muzUWzaDj4YHn8cHnoItmyBo48O4b5pU9SViYhkTPN5UtQMRo0KZ+tnnx2GX/r2hZdfjroyEZGMaD6BXmXXXeH668Owy2efweDBcNttUVclItJozS/Qq5SWhgunQ4fCxIkwYQJs3hx1VSIiaWu+gQ7Qrh3MmweXXAIzZ8KwYfBB9u62FBHJpuYd6AAtW8L06XDPPVBRAYMGwZtvRl2ViEiDKdCrnHwy/PnP8M474WGlRYuirkhEpEEU6NUNGQLPPgsFBWGM/cUXo65IRCRlCvSaevUKod6+PQwfDs89F3VFIiIpUaDXpksXePpp6NABjjsOnnkm6opEROqlQK9Lp04h1IuLYeRIDb+ISM5ToO9Ix47hAaQOHeD442HJkqgrEhGpkwK9Ph07whNPhCdMjz02TM8rIpKDFOipKCkJk3tt3RoulK5fH3VFIiJfoUBP1YEHwqOPwvvvwwknwKefRl2RiMiXKNAbIpEIT5QuXgzf+Y4WyxCRnKJAb6iRI+HWW8PZ+plnhhVMRURyQKuoC8hLZ5wBf/87XHkl9OgBF1wQdUUiIgr0tF1xBaxaBRddBAccAN/8ZtQViUgzpyGXdJnB//5vGFf/3vdg6dKoKxKRZq7eQDezXczsJTNbbGbLzeyKWtqMM7P3zOyV5HZGdsrNMbvuGtYp3X33cIb+/vtRVyQizVgqZ+ibgWHufgjQBxhhZgNqaTfb3fskt5mZLDKn7bNPCPV//hO++92wCLWISATqDXQPKpNvC5Kbbu2orl8/mDEjTBNw0UVRVyMizZR5CrfdmVlLoALoBtzs7hfV2D8OuAp4D3gN+Im7v1XL70wEJgIUFRUlysrK0iq6srKSwsLCtL6bTd1uuIHODz7Iiksv5d1jjknpO7nal3SoL7knLv0A9aVKaWlphbv3q3Wnu6e8AW2BcqB3jc/3BHZOvv4h8FR9v5VIJDxd5eXlaX83qz7/3P3II91bt3Z/5ZWUvpKzfUmD+pJ74tIPd/WlCrDA68jVBt3l4u4fJwN9RI3PP3D3zcm3M4FEQ343NgoK4N57oW1bGD0aPvkk6opEpBlJ5S6XvcysbfJ1a2A48GqNNh2rvR0FrMxgjfmlqAhmz4Y33oAf/EBPkopIk0nlDL0jUG5mS4CXgcfdfa6ZTTOzUck2ZydvaVwMnA2My065eeLII+Gaa+CBB+C666KuRkSaiXqfFHX3JUDfWj6fWu31FGBKZkvLc+edB3/9a7jrZcAAGDQo6opEJOb0pGi2mMGsWWEu9TFj9NCRiGSdAj2b9tgjTLf77rswdixs2xZ1RSISYwr0bDv00DCOPm8eXHtt1NWISIwp0JvCj34UbmOcMgWefz7qakQkphToTcEMZs6ELl3glFPgo4+irkhEYkiB3lT22APuvhvWrYMJE3R/uohknAK9KR1+OEyfDvffHybzEhHJIAV6U5s8GY49Fs49F5Yti7oaEYkRBXpTa9EC7rwzDMGMGUOLzZvr/46ISAoU6FEoKoI77oDly+l6yy1RVyMiMaFAj8pxx8F559HpoYfCikciIo2kQI/S//wPn3bvHmZlXLcu6mpEJM8p0KO0886suOwy+OwzOO00TQ0gIo2iQI/Yv7p0geuvD+uRamoAEWkEBXouOP10+Na34NJLYeHCqKsRkTylQM8FZnDbbbD33mFqgI0bo65IRPKQAj1XtGsHd90Fr78OP/lJ1NWISB5SoOeS0lK48MJwtv7gg1FXIyJ5JpVFoncxs5fMbHFy3dArammzs5nNNrPVZvaimZVkpdrmYNq0MIf6GWfAP/4RdTUikkdSOUPfDAxz90OAPsAIMxtQo83pwEfu3g34FXBNRqtsTnbaCf74x3Aro1Y5EpEGqDfQPahMvi1IbjXnfj0RuCP5+j7gaDOzjFXZ3PToAb/+NTzxBPzqV1FXIyJ5wjyFebnNrCVQAXQDbnb3i2rsXwaMcPe3k+/XAIe7+/s12k0EJgIUFRUlysrK0iq6srKSwsLCtL6ba+rsizu9pk5lzxdfZOFvfkNlt25NX1wDNYvjkmfi0g9QX6qUlpZWuHu/Wne6e8ob0BYoB3rX+HwZ0Lna+zVA+x39ViKR8HSVl5en/d1cs8O+vPeee8eO7gce6L5xY5PVlK5mc1zySFz64a6+VAEWeB252qC7XNz942Sgj6ixax1QDGBmrYA9gA8a8ttSi/btw1S7K1eGu19ERHYglbtc9jKztsnXrYHhwKs1ms0BxiZfjwaeSv5NIo11zDFw/vlw880wd27U1YhIDkvlDL0jUG5mS4CXgcfdfa6ZTTOzUck2twN7mtlq4Dzg4uyU20xNnw59+sD48fDPf0ZdjYjkqFb1NXD3JUDfWj6fWu31Z8DJmS1N/m3nncOtjIkEjBsH8+aFlY9ERKpRKuSLAw+E666Dxx6DG26IuhoRyUEK9Hzywx/CqFFw0UWweHHU1YhIjlGg5xMzmDkzTOT1ve/Bpk1RVyQiOUSBnm/22ivcyrhiBUyeHHU1IpJDFOj5aPjwEOa33AJ/+lPU1YhIjlCg56vp08OsjKefrgWmRQRQoOev6rMynnoqbN0adUUiEjEFej7r0QNuvBHKy+HnP4+6GhGJmAI9340fD9/9Lvz3f8MLL0RdjYhESIGe78zg1luhuDgsMP3JJ1FXJCIRUaDHQdu2YTz9rbfCw0eaF02kWVKgx8XAgWE90tmz4fbbo65GRCKgQI+Tiy8O0+2efTYsXx51NSLSxBTocdKiBdx1F+y2W7hQqqkBRJoVBXrcdOgQQn35cjj33KirEZEmpECPo2OPDcMvt90Gd98ddTUi0kQU6HE1bRoMGgQTJ8Jrr0VdjYg0AQV6XBUUhLPznXeG73wnTBEgIrGWyiLRxWZWbmYrzGy5mZ1TS5uhZvaJmb2S3KbW9lvSxIqLw1S7ixdrPF2kGah3TVFgC3C+uy80s92ACjN73N1X1Gj3F3c/IfMlSqOMHAkXXhjmehkyJCyMISKxVO8Zuruvd/eFydefAiuBTtkuTDLoZz+DwYPDePrKlVFXIyJZ0qAxdDMrAfoCL9aye6CZLTazR8ysVyaKkwwpKICyMth1Vxg9GjZujLoiEckC8xTn/TCzQuBpYLq7P1Bj3+7ANnevNLORwPXu3r2W35gITAQoKipKlJWVpVV0ZWUlhYWFaX031zRlX75WUcHBF1zAO8ccw6tTpoSJvTJIxyX3xKUfoL5UKS0trXD3frXudPd6N6AAeAw4L8X2bwLtd9QmkUh4usrLy9P+bq5p8r5ccYU7uN9yS8Z/Wscl98SlH+7qSxVggdeRq6nc5WLA7cBKd7+ujjYdku0ws/6EoZwPGvo3jzSByy6D448P8728WNvImYjkq1TuchkEnAosNbNXkp9dAnQBcPdbgdHAmWa2BfgXMCb5N4nkmhYt4Pe/h0QijKcvXAh77RV1VSKSAfUGurs/C+xwsNXdbwJuylRRkmXt2sH998MRR4RFMR57DFq2jLoqEWkkPSnaXB16KNx8Mzz5JFx6adTViEgGKNCbs9NPDyscXXMN3Htv1NWISCMp0Ju7668Pqx2NHw/LlkVdjYg0ggK9udt5Z7jvvrAoxkknwUcfRV2RiKRJgS6wzz7hIunf/w5jxsCWLVFXJCJpUKBLcMQR4SLpn/8cFscQkbyTyn3o0lxMmBCm2r32WjjkEDj11KgrEpEG0Bm6fNmvfgVDh4Zw15OkInlFgS5fVlAQbmHcZ59wkfStt6KuSERSpECXr2rfHh5+OEyzO2qUptsVyRMKdKldr14wezYsWRLG0rdti7oiEamHAl3qdvzx4QLpgw/CJZdEXY2I1EN3uciOnXMOrFoVpgfo1g3OOCPqikSkDgp02TEzuPFGePNNmDQJ9t0Xhg+PuioRqYWGXKR+rVqF8fSePcMc6przRSQnKdAlNbvvDv/3f9CmDYwcCevWRV2RiNSgQJfUFRfDvHlhAq+RI2HDhqgrEpFqFOjSMH36hIm8VqyAb38bPv886opEJEmBLg137LEwcyY88QT84Ae6R10kR9Qb6GZWbGblZrbCzJab2Tm1tDEzu8HMVpvZEjM7NDvlSs4YOxamT4c//AEmTwatCS4SuVRuW9wCnO/uC81sN6DCzB539xXV2hwPdE9uhwO3JP8rcTZlCrzzTpjQq0MH6N8/6opEmrV6z9Ddfb27L0y+/hRYCXSq0exE4E4PXgDamlnHjFcrucUshPkpp8BFF9HhkUeirkikWTNvwD+VzawEeAbo7e4bqn0+F7ja3Z9Nvn8SuMjdF9T4/kRgIkBRUVGirKwsraIrKyspLCxM67u5Jg59sS++4KBLLuFrCxeyfOpU3j/qqKhLarQ4HBeITz9AfalSWlpa4e79at3p7iltQCFQAXyrln1zgcHV3j8J9NvR7yUSCU9XeXl52t/NNbHpS2Wlf9yrl3tBgfujj0ZdTaPF5bjEpR/u6ksVYIHXkasp3eViZgXA/cAf3P2BWpqsA4qrve+c/EyaizZtWHr11WGWxv/4D3juuagrEml2UrnLxYDbgZXufl0dzeYApyXvdhkAfOLu6zNYp+SBLYWF8Nhj4QGkkSNhwYL6vyQiGZPKGfog4FRgmJm9ktxGmtkkM5uUbDMPWAusBm4DfpSdciXn7b03PPkk7LlnuF998eKoKxJpNuq9bdHDhU6rp40DP85UUZLnOneGp56CIUPgmGPg6afDxF4iklV6UlSyo6QknKkXFMCwYbByZdQVicSeAl2yp3v3cKZuBqWlCnWRLFOgS3YdcACUl4fXCnWRrFKgS/YdcADMnx9eDx2qBTJEskSBLk2jKtRbtQqhvmhR1BWJxI4CXZrOAQfAM8+EVY+GDYOXXoq6IpFYUaBL0+raNYR6u3bbb2kUkYxQoEvT23ffEOrFxTBiRFirVEQaTYEu0ejUKZyd9+4NJ50Ed98ddUUieU+BLtFp3z48fDRoEHz/+3DTTVFXJJLXFOgSrd13h0cegRNPhP/6L7jsMi1nJ5ImBbpEr3VruPdemDAhrFM6YQJs2RJ1VSJ5J5U1RUWyr1Ur+O1vw9qkV14J69fD7NkQkxVqRJqCztAld5jBtGkh2B97DI46KgS7iKREgS65Z+JEmDMHVq2CAQM0VYBIihTokptGjgy3NX7xBRxxRLhwKiI7pECX3JVIhOkBunaFE06AG2+MuiKRnKZAl9zWuTP85S8h0M8+GyZNgs8/j7oqkZyUyiLRs8zsXTOrdSDTzIaa2SfV1hudmvkypVkrLIQHHoCLLw4XTI85Bt57L+qqRHJOKmfovwNG1NPmL+7eJ7lNa3xZIjW0bAlXXQV//CO8/DL066cpeEVqqDfQ3f0Z4MMmqEWkfqecAs8+C9u2hYuld94ZdUUiOSNTY+gDzWyxmT1iZr0y9JsitUskoKICBg6EsWPhrLM0ri4CmKcwb4aZlQBz3b13Lft2B7a5e6WZjQSud/fudfzORGAiQFFRUaKsrCytoisrKymMyROE6kv6bOtW9rvtNrrMns2GAw9k+U9/yuaiooz8dlyOS1z6AepLldLS0gp371frTnevdwNKgGUptn0TaF9fu0Qi4ekqLy9P+7u5Rn3JgHvvdd9tN/d27dznzcvIT8bluMSlH+7qSxVggdeRq40ecjGzDmZmydf9CcM4HzT2d0VSNnp0GILp3Dk8kDRlSnggSaSZSeW2xbuB54EeZva2mZ1uZpPMbFKyyWhgmZktBm4AxiT/FhFpOt27wwsvhJkar746zAPz5ptRVyXSpOqdbdHdT6ln/02AViaQ6LVuDTNmwNFHh/lg+vSB226Dk0+OujKRJqEnRSV+vvvdcI96jx7wne/A+PHw6adRVyWSdQp0iaf99w/3q192WbhXvU8feP75qKsSySoFusRXQUFYLOPpp2HrVhg8GC65RPesS2wp0CX+Bg+GJUtg3LgwfcBhh8HixVFXJZJxCnRpHnbfHW6/PSyc8c47YS6Yyy/X2brEigJdmpdvfhOWLw8XTq+4IpytL1wYdVUiGaFAl+Znzz3h97+Hhx4K0/D27w8XXACbNkVdmUijKNCl+Ro1ClasgB/8AH75S+jdGx5/POqqRNKmQJfmrW3b8DDS/Pnhrphjj4VTTmGnDzR7heQfBboIhKkCFi8OF0offJD+Y8eGNUy3bIm6MpGUKdBFquyyC/z0p7B0KRsOPDCsYZpIhDVNRfKAAl2kpu7dWfLzn8N998FHH8GQIfD978Pbb0ddmcgOKdBFamMG3/42vPpqmD7g/vvD3DDTpuluGMlZCnSRHdl11zB9wKuvwgknhCGZHj3C/DDbtkVdnciXKNBFUlFSArNnwzPPQMeOYS3TRAKeeCLqykT+TYEu0hBHHhkW0vjjH8P4+vDh4VZHPW0qOUCBLtJQLVrAKaeEYZhrrw3L3yUSMGYMrFoVdXXSjCnQRdK1yy5w3nmwdm24cDp3LvTsGWZ1XLs26uqkGVKgizTWHnuEC6dr18K554ax9h494IwzFOzSpFJZJHqWmb1rZsvq2G9mdoOZrTazJWZ2aObLFMkDe+8dhmDWrIFJk8IEYF//epgrZvXqqKuTZiCVM/TfASN2sP94oHtymwjc0viyRPLYPvuEaQPWrIEf/zhcQO3RIzyctKzW8yKRjKg30N39GeDDHTQ5EbjTgxeAtmbWMVMFiuStTp3g+uvhzTfh/PPD4hoHHRRmeXzuuairkxgyd6+/kVkJMNfde9eyby5wtbs/m3z/JHCRuy+ope1Ewlk8RUVFibKysrSKrqyspLCwMK3v5hr1JTdloy+tNmyg04MP0vmBByjYsIFPevXirTFjeH/gQGjZMqN/VhUdk9zUmL6UlpZWuHu/Wne6e70bUAIsq2PfXGBwtfdPAv3q+81EIuHpKi8vT/u7uUZ9yU1Z7UtlpfuNN7qXlLiDe7du7jfdFD7PMB2T3NSYvgALvI5czcRdLuuA4mrvOyc/E5HatGkDZ50Fr78O99wTVlA66yzo3BkmTw5DNCJpyESgzwFOS97tMgD4xN3XZ+B3ReKtVSs4+WR4/vkwpn7ssfDrX0PXrnDSSWH1pBSGREWqpHLb4t3A80APM3vbzE43s0lmNinZZB6wFlgN3Ab8KGvVisSRGRxxRLh//Y034MILtwf8AQeEkP9wR/cliASp3OVyirt3dPcCd+/s7re7+63ufmtyv7v7j929q7sf5LVcDBWRFBUXw1VXwVtvwV13Qbt28JOfhDtmTjstBL3O2qUOelJUJBftsgv853+G4ZhXXoHx4+FPf4LBg8P0Ar/8Jbz7btRVSo5RoIvkukMOgd/8Bv7xD5g1K5y1X3BBOGs/6SR46CH44ouoq5QcoEAXyReFheFM/bnnYPnyMG/MCy+EUO/UKayBumCBhmSaMQW6SD7q2RN+8YuwzunDD8PQoTBjBhx2WNh35ZWaP6YZUqCL5LNWrcLSePfcA//8Zwj1oiKYOhW6d+fQM8+E667TAtfNhAJdJC7atoUJE2D+fPj73+EXv8C2bg3zyBQXh9WWbrgB1um5v7hSoIvEUXExTJ5MxYwZ8NprYQjm44/hnHPCE6mDBoWpfjVfe6wo0EXirnv3sKLS0qVh2byf/Qw2bQrTDHTtCn36wOWXw6JFuqCa5xToIs1Jjx5w6aUhvNeuDWfphYUwbRoceiiUlIQ53B99FD77LOpqpYEU6CLN1X77hTVRn302XFCdNQv69oXf/Q6OPx7at4cTT4Tf/jY8uSo5r1XUBYhIDth773CP+/jx8K9/hQurDz8M8+aFhTkAeveG446DESPCE6u77BJpyfJVCnQR+bLWrcMZ+vHHhzH1lStDsD/6aFha79prQ5shQ2D48LD17g0t9A/+qCnQRaRuZuFBpZ49w0XUjRvD2fuf/xym9508ObTbay8YNgyOPhpKS8PFVrNIS2+OFOgikro2beAb3wgbhLH1J5/cvs2eHT4vLg7BftRRYdt/fwV8E1Cgi0j6ioth3LiwucOqVVBeDk89BY88AnfeGdp16hQebBoyJPy3Z08N0WSBAl1EMsMsLMhxwAFw5pnbx9+ffhqeeSZsVQvDt20bFvUYNCj897DDwtm/NIoCXUSyo/r4e1XAv/EG/OUvYcbIZ58NF1sBWrYM0wQPHAgDBsDhh0O3bhqmaSAFuog0DbMwlr7//jB2bPjsww/DFMB//WvY7rgDbr457Pva16B//7Addhg7ac73eqUU6GY2ArgeaAnMdPera+wfB/wCqJr15yZ3n5nBOkUkjtq1g5EjwwawdSusWBFC/qWXwjZ9OmzbxhEQxuL79YNEYvtWVBRlD3JKvYFuZi2Bm4HhwNvAy2Y2x91X1Gg6293PykKNItJctGwJBx0UtgkTwmcbN8KiRawuK6PbRx9BRUVYpanKPvuEJ1z79g3z0vTtG6YwaIYXXVM5Q+8PrHb3tQBmVgacCNQMdBGRzGvTBgYP5u0tW+g2dGj4bMOGMB/NokWwcGHYHn00nOED7LYbHHxwGJc/+OCw9e4dPo8x83pmVzOz0cAIdz8j+f5U4PDqZ+PJIZergPeA14CfuPtXJn8ws4nARICioqJEWdUV7waqrKyksLAwre/mGvUlN8WlL3HpB9TflxabN9PmjTcofP11Cteupc2aNRSuXUurjRv/3eZfHTqwcf/92VhSwsaSEjbttx+bunRh2047NUUX/q0xx6W0tLTC3fvVti9TF0UfBu52981m9kPgDmBYzUbuPgOYAdCvXz8fWvW3bQPNnz+fdL+ba9SX3BSXvsSlH5BmX9zhb38LUwcvWULrpUtpvWwZ7e+5B7ZsCW1atAhPtvbqFe7IOfDAsPXoEWaizIJsHZdUAn0dUFztfWe2X/wEwN0/qPZ2JvDzxpcmItJIZmE8vaQEvvnN7Z9//nlY+GP58rCtWBG2uXO3Bz2EB6eq7q3v0SNsX/96WCQkB8foUwn0l4HuZrYfIcjHAN+r3sDMOrr7+uTbUcDKjFYpIpJJO+0UxtR79/7y559/HhbXfvXV8FDUq6+Gp19/9zv49NPt7Vq3DvfJd+8eAr7qdbdu0LFjZPfP1xvo7r7FzM4CHiPctjjL3Zeb2TRggbvPAc42s1HAFuBDYFwWaxYRyY6ddtr+MFR17rB+fTirX7UqbK+/Hs7uH34Yqt8jv+uu4V77bt3CUE7XruF9167QpUv4M7IkpTF0d58HzKvx2dRqr6cAUzJbmohIjjALt0fusw/UHPvesiUsyv3662FbsyZsq1aF+Ww2b97etkUL6NyZzt/4xld/JwP0pKiISGO0arX9Cdjjjvvyvm3bwpn9mjVh2oO1a+GNN/h8zz2zU0pWflVERMIZeadOYRsy5N8fvzt/Pj138LW0/7gs/KaIiERAgS4iEhMKdBGRmFCgi4jEhAJdRCQmFOgiIjGhQBcRiQkFuohITNQ7H3rW/mCz94C/pfn19sD7GSwnSupLbopLX+LSD1Bfquzr7nvVtiOyQG8MM1tQ1wTv+UZ9yU1x6Utc+gHqSyo05CIiEhMKdBGRmMjXQJ8RdQEZpL7kprj0JS79APWlXnk5hi4iIl+Vr2foIiJSgwJdRCQmcjrQzWyEma0ys9VmdnEt+3c2s9nJ/S+aWUkEZaYkhb6MM7P3zOyV5HZGFHXWx8xmmdm7Zrasjv1mZjck+7nEzA5t6hpTlUJfhprZJ9WOydTa2kXNzIrNrNzMVpjZcjM7p5Y2eXFcUuxLvhyXXczsJTNbnOzLFbW0yWyGuXtOboQFqdcA+wM7AYuBnjXa/Ai4Nfl6DDA76rob0ZdxwE1R15pCX4YAhwLL6tg/EngEMGAA8GLUNTeiL0OBuVHXmUI/OgKHJl/vBrxWy/++8uK4pNiXfDkuBhQmXxcALwIDarTJaIbl8hl6f2C1u69198+BMuDEGm1OBO5Ivr4PONrMrAlrTFUqfckL7v4M8OEOmpwI3OnBC0BbM+vYNNU1TAp9yQvuvt7dFyZffwqsBDrVaJYXxyXFvuSF5P+tK5NvC5JbzbtQMpphuRzonYC3qr1/m68e2H+3cfctwCdAdlZfbZxU+gLw7eQ/h+8zs+KmKS3jUu1rvhiY/CfzI2bWK+pi6pP8J3tfwtlgdXl3XHbQF8iT42JmLc3sFeBd4HF3r/O4ZCLDcjnQm5uHgRJ3Pxh4nO1/a0t0FhLmzTgEuBH4U7Tl7JiZFQL3A+e6+4ao62mMevqSN8fF3be6ex+gM9DfzHpn88/L5UBfB1Q/S+2c/KzWNmbWCtgD+KBJqmuYevvi7h+4++bk25lAoolqy7RUjltecPcNVf9kdvd5QIGZtY+4rFqZWQEhAP/g7g/U0iRvjkt9fcmn41LF3T8GyoERNXZlNMNyOdBfBrqb2X5mthPhgsGcGm3mAGOTr0cDT3ny6kKOqbcvNcYzRxHGDvPRHOC05F0VA4BP3H191EWlw8w6VI1nmll/wv+/5NwJQ7LG24GV7n5dHc3y4rik0pc8Oi57mVnb5OvWwHDg1RrNMpphrdL9Yra5+xYzOwt4jHCXyCx3X25m04AF7j6HcODvMrPVhItbY6KruG4p9uVsMxsFbCH0ZVxkBe+Amd1NuMugvZm9DfyUcLEHd78VmEe4o2I1sAkYH02l9UuhL6OBM81sC/AvYEyOnjAMAk4FlibHawEuAbpA3h2XVPqSL8elI3CHmbUk/KVzj7vPzWaG6dF/EZGYyOUhFxERaQAFuohITCjQRURiQoEuIhITCnQRkZhQoIuIxIQCXUQkJv4fZZEL8k+i0voAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3.1415926535897936, 3.4878684980086326e-14)\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "\n",
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "x = np.linspace(0,3, 100)\n",
    "y = func(x)\n",
    "\n",
    "plt.plot(x, y, color = 'r')\n",
    "\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "print(integrate.quad(func, 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "なんでと思うかもしれないが，\n",
    "```maple\n",
    ">int(1/(1+x^2),x);\n",
    "```\n",
    "$$\n",
    "arctan(x)\n",
    "$$\n",
    "となるので，納得できるでしょう．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Midpoint rule(中点法) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.142894729591689\n"
     ]
    }
   ],
   "source": [
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "N, x0, xn = 8, 0.0, 1.0\n",
    "\n",
    "h = (xn-x0)/N\n",
    "S = 0.0\n",
    "for i in range(0, N):\n",
    "    xi = x0 + (i+0.5)*h\n",
    "    dS = h * func(xi)\n",
    "    S = S + dS\n",
    "\n",
    "print(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Trapezoidal rule(台形公式) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "3.138988494491089\n"
     ]
    }
   ],
   "source": [
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "N, x0, xn = 8, 0.0, 1.0\n",
    "\n",
    "h = (xn-x0)/N\n",
    "S = func(x0)/2.0\n",
    "for i in range(1, N):\n",
    "    xi = x0 + i*h\n",
    "    dS = func(xi)\n",
    "    S = S + dS\n",
    "    print(\"{0}\".format(i))\n",
    "\n",
    "S = S + func(xn)/2.0\n",
    "print(h*S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simpson's rule(シンプソンの公式) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "3\n",
      "5\n",
      "7\n",
      "2\n",
      "4\n",
      "6\n",
      "3.141592502458707\n"
     ]
    }
   ],
   "source": [
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "N, x0, xn = 8, 0.0, 1.0\n",
    "\n",
    "M = int(N/2)\n",
    "h = (xn-x0)/N\n",
    "Seven, Sodd = 0.0, 0.0\n",
    "for i in range(1, 2*M, 2): #rangeの終わりに注意\n",
    "    xi = x0 + i*h\n",
    "    Sodd += func(xi)\n",
    "    print(\"{0}\".format(i))\n",
    "for i in range(2, 2*M, 2):\n",
    "    xi = x0 + i*h\n",
    "    Seven += func(xi)\n",
    "    print(\"{0}\".format(i))\n",
    "\n",
    "print(h*(func(x0)+4*Sodd+2*Seven+func(xn))/3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 課題\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  補間法\n",
    "次の4点\n",
    "```maple\n",
    "x y \n",
    "0 1 \n",
    "1 2\n",
    "2 3\n",
    "3 1\n",
    "```\n",
    "を通る多項式を逆行列で求めよ．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 対数関数のニュートンの差分商補間(2014期末試験，25点) \n",
    "\n",
    "2を底とする対数関数(Mapleでは$\\log[2](x)$)の$F(9.2)=2.219203$をニュートンの差分商補間を用いて求める．ニュートンの内挿公式は，\n",
    "\n",
    "$$\n",
    "\\begin{array}{rc}\n",
    "F (x )&=F (x _{0})+\n",
    "(x -x _{0})f _{1}\\lfloor x_0,x_1\\rfloor+\n",
    "(x -x _{0})(x -x _{1})\n",
    "f _{2}\\lfloor x_0,x_1,x_2\\rfloor + \\\\\n",
    "& \\cdots +  \\prod_{i=0}^{n-1} (x-x_i) \\, \n",
    "f_n \\lfloor x_0,x_1,\\cdots,x_n \\rfloor\n",
    "\\end{array}\n",
    "$$\n",
    "である．ここで$f_i \\lfloor\\, \\rfloor$ は次のような関数を意味していて，\n",
    "$$\n",
    "\\begin{array}{rc}\n",
    "f _{1}\\lfloor x_0,x_1\\rfloor &=&  \\frac{y_1-y_0}{x_1-x_0} \\\\\n",
    "f _{2}\\lfloor x_0,x_1,x_2\\rfloor &=&  \\frac{f _{1}\\lfloor x_1,x_2\\rfloor-\n",
    "f _{1}\\lfloor x_0,x_1\\rfloor}{x_2-x_0} \\\\\n",
    "\\vdots & \\\\\n",
    "f _{n}\\lfloor x_0,x_1,\\cdots,x_n\\rfloor &=&  \\frac{f_{n-1}\\lfloor x_1,x_2\\cdots,x_{n}\\rfloor-\n",
    "f _{n-1}\\lfloor x_0,x_1,\\cdots,x_{n-1}\\rfloor}{x_n-x_0} \n",
    "\\end{array}\n",
    "$$\n",
    "差分商と呼ばれる．$x_k=8.0,9.0,10.0,11.0$をそれぞれ選ぶと，差分商補間のそれぞれの項は以下の通りとなる．\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccl|lll}\n",
    "\\hline\n",
    "k  &  x_k & y_k=F_0( x_k) &f_1\\lfloor x_k,x_{k+1}\\rfloor & f_2\\lfloor x_k,x_{k+1},x_{k+2}\\rfloor &  f_3\\lfloor x_k,x_{k+1},x_{k+2},x_{k+3}\\rfloor \\\\\n",
    "\\hline\n",
    "0  &   8.0  &  2.079442  &          &              &\\\\\n",
    "&      &     &     0.117783     &              &\\\\ \n",
    "1  &   9.0  &  2.197225  &           &       [ XXX ]      &\\\\\n",
    "&      &     &     0.105360     &              & 0.0003955000 \\\\\n",
    "2  &  10.0  &  2.302585  &           &       -0.0050250      &\\\\ \n",
    "&      &     &     0.095310    &              &\\\\ \n",
    "3  &  11.0  &  2.397895 &           &              &\\\\ \n",
    "\\hline\n",
    "\\end{array}\n",
    "$$\n",
    "それぞれの項は，例えば，\n",
    "\n",
    "$$\n",
    "f_2\\lfloor x_1,x_2,x_3 \\rfloor =\\frac{0.095310-0.105360}{11.0-9.0}=-0.0050250\n",
    "$$\n",
    "で求められる．ニュートンの差分商の一次多項式の値はx=9.2で\n",
    "\n",
    "$$\n",
    "F(x)=F_0(8.0)+(x-x_0)f_1\\lfloor x_1,x_0\\rfloor =2.079442+0.117783(9.2-8.0)=2.220782\n",
    "$$\n",
    "となる．\n",
    "\n",
    "###  差分商補間の表中の開いている箇所[ XXX ]を埋めよ．\n",
    "###  ニュートンの二次多項式\n",
    "\n",
    "$$\n",
    "F (x )=F (x _{0})+(x -x _{0})f _{1}\\lfloor x_0,x_1\\rfloor+(x -x _{0})(x -x _{1})\n",
    "f _{2}\\lfloor x_0,x_1,x_2\\rfloor\n",
    "$$\n",
    "の値を求めよ．\n",
    "###  ニュートンの三次多項式の値を求めよ．\n",
    "ただし，ここでは有効数字7桁程度はとるように．(E.クライツィグ著「数値解析」(培風館,2003), p.31, 例4改)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数値積分(I)\n",
    "次の関数\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{4}{1+x^2}\n",
    "$$\n",
    "を$x = 0..1$で数値積分する．\n",
    "1.  $N$を2,4,8,…256ととり，$N$個の等間隔な区間にわけて中点法と台形則で求めよ．(15)\n",
    "1.  小数点以下10桁まで求めた値3.141592654との差をdXとする．dXと分割数Nとを両対数プロット(loglogplot)して比較せよ(10)\n",
    "(2008年度期末試験）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 解答例[7.3]\n",
    "\n",
    "以下には，中点則の結果を示した．課題では，台形則を加えて，両者を比較せよ．予測とどう違うか．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def func(x):\n",
    "    return 4.0/(1.0+x**2)\n",
    "\n",
    "def mid(N):\n",
    "    x0, xn = 0.0, 1.0\n",
    "\n",
    "    h = (xn-x0)/N\n",
    "    S = 0.0\n",
    "    for i in range(0, N):\n",
    "        xi = x0 + (i+0.5)*h\n",
    "        dS = h * func(xi)\n",
    "        S = S + dS\n",
    "    return S\n",
    "\n",
    "x, y = [], []\n",
    "for i in range(1,8):\n",
    "    x.append(2**i)\n",
    "    y.append(abs(mid(2**i)-np.pi))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfEElEQVR4nO3deZgU1b3/8feXXdkElDECss0gmyIMa9zAKIKAaGJ+aIxL9AlGE000uUauiTEukVxxh0gMIlETUaNxJZr8FDTXIAJiEIMoLmyREBUQNCrI9/5xuqVn2Bqme6qq6/N6nnqYqu7p/lb38OnqU6fOMXdHRERKX52oCxARkdqhwBcRSQkFvohISijwRURSQoEvIpIS9aIuYGf23Xdf79Chwy7v99FHH9G4cePiF1RESd8H1R+9pO+D6i+c+fPnv+fu+1XfHsvAN7NRwKjy8nLmzZu3y/vPmjWLwYMHF72uYkr6Pqj+6CV9H1R/4ZjZsu1tj2WTjrs/5u5jmzdvHnUpIiIlI5aBLyIihafAFxFJCQW+iEhKKPBFRFJCgS8ikhIKfBGRlCjNwH/4Ybj++qirEBGJldIM/CefhGuvjboKEZFYiWXgm9koM7t9/fr1e/YAFRXw/vthERERIKaBX+Mrbbt0Cf++8UbhihIRSbhYBn6NKfBFRLZRmoHfsSPUqQOvvx51JSIisVGagd+gQQh9Bb6IyBdKM/AhNOuoSUdE5AulG/gVFeEI3z3qSkREYqF0A79LF/joI1i9OupKRERioXQDv6Ii/Kt2fBERoJQDP9s1U4EvIgKUcuC3awcNG+rErYhIRukGft260LmzjvBFRDJKN/AhNOso8EVEgFoOfDM70cx+Y2b3mdnQoj9hly7w5pvw+edFfyoRkbjLO/DNbKqZrTGzRdW2DzOzJWa21Mwu3dljuPvD7v5t4DvAmD0reTdUVMBnn8Hy5UV/KhGRuNudI/xpwLDcDWZWF5gEDAe6A6eaWXczO9jMHq+2tM751Z9kfq+4NIiaiMgX6uV7R3d/zsw6VNvcH1jq7m8BmNl0YLS7XwuMrP4YZmbAeOBP7v7SHledr9y++EOL34IkIhJneQf+DrQBVuSsrwQG7OT+FwDHAM3NrNzdJ1e/g5mNBcYClJWVMWvWrF0WsXHjxu3fz53D99qL1c88w9KePXf5OFHa4T4khOqPXtL3QfXXAnfPewE6AIty1k8GpuSsnw5M3J3H3NlSWVnp+Zg5c+aOb+zTx33YsLweJ0o73YcEUP3RS/o+qP7CAeb5djK1pr10VgHtctbbZrbVSI2nOMyVHURNRCTlahr4c4EKM+toZg2AU4BHa1qU13SKw1xdusA774TeOiIiKbY73TLvBWYDB5nZSjM7x903A98DngIWA/e7+6vFKXUPdekCW7bAW29FXYmISKR2p5fOqTvYPgOYUbCKCE06wKjy8vKaP1huT52uXWv+eCIiCRXLoRUK2qSTDXz1xReRlItl4BdUy5bQqpVO3IpI6sUy8AvaSwc0iJqICDEN/II26YAmNBcRIaaBX3AVFbBqVZjjVkQkpdIR+NlB1JYujbYOEZEIxTLwi9KGD2rHF5FUi2XgF7wNP9ufX4EvIikWy8AvuMaNoU0bBb6IpFo6Ah+gVy+YPTvqKkREIhPLwC94Gz7A8OGha6aO8kUkpWIZ+AVvwwcYMSL8+8QThXtMEZEEiWXgF0XHjtCtmwJfRFIrPYEP4Sj/uedgw4aoKxERqXXpC/xNm+Avf4m6EhGRWhfLwC/KSVuAww6D5s3VrCMiqRTLwC/KSVuA+vVh6FCYMSPMgiUikiKxDPyiGjECVq+GBQuirkREpFalL/CHDwczNeuISOqkL/Bbt4Z+/RT4IpI66Qt8CM06c+fCmjVRVyIiUmvSG/ju8Kc/RV2JiEitiWXgF61bZlbv3vClL6lZR0RSJZaBX7RumVl16sDxx8NTT4ULsUREUiCWgV8rRoyADz+E55+PuhIRkVqR3sA/5phwIZaadUQkJdIb+E2bwlFHKfBFJDXSG/gQmnUWL4a33466EhGRolPgg47yRSQV0h34FRVhUeCLSAqkO/AhHOXPnAkffRR1JSIiRaXAHzECPv0Unnkm6kpERIoqloFf9Cttcx15JDRpomYdESl5sQz8ol9pm6tBAzj22BD47sV/PhGRiMQy8GvdiBGwciUsXBh1JSIiRaPAhzCuDqhZR0RKmgIfwsiZffoo8EWkpCnws0aOhBdegPffj7oSEZGiUOBnjRgBW7bAk09GXYmISFEo8LP69g3z3apZR0RKlAI/q04dGD48HOFv3hx1NSIiBafAzzViBKxdC3/7W9SViIgUnAI/13HHQYsWcOWVughLREqOAj9Xs2Zw9dXw9NPw0ENRVyMiUlAK/OrGjoVDDoGLL4aPP466GhGRgqm1wDezbmY22cz+YGbn1dbz7rZ69eCWW2D5crjuuqirEREpmLwC38ymmtkaM1tUbfswM1tiZkvN7NKdPYa7L3b37wD/Dzhsz0uuBUcdBWPGwPjxsGxZ1NWIiBREvkf404BhuRvMrC4wCRgOdAdONbPuZnawmT1ebWmd+Z0TgCeAGQXbg2K57jowgx/9KOpKREQKwjzP3ihm1gF43N17ZtYHAVe4+3GZ9XEA7n5tHo/1hLuP2MFtY4GxAGVlZZXTp0/fZW0bN26kSZMmee3H7mh/9910nDqVl6+/nnV9+hT88XMVax9qi+qPXtL3QfUXzpAhQ+a7e99tbnD3vBagA7AoZ/1kYErO+unAxJ38/mDgFuDXwHfzec7KykrPx8yZM/O63277z3/cO3Z079HDfdOm4jxHRtH2oZao/uglfR9Uf+EA83w7mVprJ23dfZa7X+ju57r7pNp63hpp1AhuvBFefRVuuy3qakREaqQmgb8KaJez3jazrcZqdYrDXTnhBBg6FC6/HP7976irERHZYzUJ/LlAhZl1NLMGwCnAo4UoymtzisNdMYObboKNG+Gyy6KuRkRkj+XbLfNeYDZwkJmtNLNz3H0z8D3gKWAxcL+7v1q8UiPUrRtceCFMmQLz50ddjYjIHqmXz53c/dQdbJ9BEbpYmtkoYFR5eXmhH3rPXX453HMPXHABPP98OPIXEUmQWA6tEKsmnazmzcOFWLNnw+9+F3U1IiK7LZaBH1tnngn9+sEll8CGDVFXIyKyW2IZ+LHqpZOrTh249VZ491245pqoqxER2S2xDPxYNulkDRgAZ50FN9wAb7wRdTUiInmLZeDH3rXXhouyLroo6kpERPIWy8CPbZNO1v77w89+FiY816TnIpIQsQz8WDfpZF1wAXTtCj/4AXz6adTViIjsUiwDPxEaNICbb4alS8OVuCIiMafAr4mhQ2H0aLjqKvjnP6OuRkRkpxT4NXXDDbB5M/z4x1FXIiKyU7EM/NiftM3VqVOYFeuee8KQCyIiMRXLwE/ESdtc48ZBmzZhgLXPP4+6GhGR7Ypl4CdO48YwYQK89BJMnRp1NSIi26XAL5QxY+CII+C//xvWro26GhGRbSjwC8UMbrkFPvgArrgi6mpERLYRy8BP1EnbXIceCueeC5MmwaJFUVcjIlJFLAM/cSdtc111VRg7/8ILwT3qakREvhDLwE+0Vq3g6qth5kx48MGoqxER+YICvxjGjoVeveCHP4SPP466GhERQIFfHHXrhhO4y5eHXjtq2hGRGFDgF8uRR8J554UB1r71Lfjss6grEpGUqxd1ASVt0iQoKwvdNN95Bx56CFq2jLoqEUmpWB7hJ7ZbZnVmYaKUe+6B2bNh4EBNiygikYll4Ce6W+b2nHYaPP10uChr4ED461+jrkhEUiiWgV+SDj8c5syB/faDr3wF7r476opEJGUU+LWpc+fQtHP44XDGGXD55erBIyK1RoFf21q0gCefhLPPDlflfuMb8MknUVclIimgwI9CgwYwZQqMHw/Tp8PRR1N/3bqoqxKREqfAj4pZmBbxgQdgwQL6nH8+LF4cdVUiUsIU+FE7+WR49lnqfvIJDBoUevOIiBSBAj8O+vdn/q9+Be3awbBh8JvfRF2RiJQgBX5MfLr//mES9GOOCYOvXXIJbNkSdVkiUkJiGfglc6Xt7mrWDB57LIzBc911oblHo22KSIHEMvBL7krb3VGvXhiD58Yb4eGH4aij4N13o65KREpALAM/9czgBz+ARx4JPXcGDICFC6OuSkQSToEfZ6NGhXF3tmyBww6DGTOirkhEEkyBH3e9e4cxeCoqwgfAxIlRVyQiCaXAT4I2beC552DkSLjggjBB+uefR12ViCSMAj8pmjQJE6hcfDHceiuMHg0bNkRdlYgkiAI/SerWheuvh9tuCwOwHXEErFgRdVUikhAK/CT6znfCCdy33w49eObPj7oiEUkABX5SDR0arsxt0CAc6T/8cNQViUjMKfCTrGfP0IPnkEPgq1+FCRM0oYqI7JACP+nKymDmzDAMw3/9V2ju2bQp6qpEJIbqRV2AFMBee4WJVCoq4Be/gLfeCuPs77NP1JWJSIzU6hG+mTU2s3lmNrI2nzcV6tSBa66BO++EZ5+FL385nNQVEcnIK/DNbKqZrTGzRdW2DzOzJWa21MwuzeOhfgzcvyeFSp7OOgv+/GdYvTr04Jk9O+qKRCQm8j3CnwYMy91gZnWBScBwoDtwqpl1N7ODzezxaktrMzsW+AewpoD1y/YMHhyCvlkzGDIE7rsv6opEJAbM8+zVYWYdgMfdvWdmfRBwhbsfl1kfB+Du1+7g968BGhM+HP4DnOTu28zwYWZjgbEAZWVlldOnT99lbRs3bqRJkyZ57UdcFWMf6q9fT4/LL2efhQt5++yzWfbNb4aROIsg6e9B0uuH5O+D6i+cIUOGzHf3vtvc4O55LUAHYFHO+snAlJz104GJeTzOWcDIfJ6zsrLS8zFz5sy87hdnRduHTz5x/+Y33cH9jDPCehEk/T1Iev3uyd8H1V84wDzfTqbWei8dd59W28+Zag0bwl13QZcucPnl8M47YUyeVq2irkxEallNeumsAtrlrLfNbKux1E5xWCxm8NOfwu9/Dy+8AAMHwuuvR12ViNSymgT+XKDCzDqaWQPgFODRQhTlaZ7isJhOPRWeeQbWrQuhf9NNoA9VkdTIt1vmvcBs4CAzW2lm57j7ZuB7wFPAYuB+d3+1EEXpCL+IDjssDMfQowdcdFEYa/+73w1TKYpIScsr8N39VHf/krvXd/e27n5HZvsMd+/i7p3d/ZpCFaUj/CLr1ClMnThvXhiSYcoU6N4djj0WHn1Uk6uIlCiNpZNmlZUwbVoYU//qq8NR/ujRYYiGCRNg7dqoKxSRAlLgC7RuDZddFoZiuP9+aNs2DMTWti2cey4sWrTrxxCR2Itl4KsNPyL168PXvx7mz12wIJzkvesuOPjgcMXuQw/B5s1RVykieyiWga82/Bg49NDQtr9yJYwfH0bg/NrXoHNn+OUv4f33o65QRHZTLANfYqRVK/jxj+HNN8MRfufOcOmlobnnnHPg5ZejrlBE8hTLwFeTTgzVqwcnnRT68b/yCpx5ZhiDv3dvOOII9ps5UxOviMRcLANfTTox17MnTJ4cmnsmTIBVq+hx5ZXQsWMYk3+NBkQViaNYBr4kRIsW8MMfwhtv8Mo114S+/D/5CbRrF74BzJsXdYUikkOBLzVXty7vf/nLYeKVxYvh29+GBx+Efv1g0KAwhs9nn0VdpUjqKfClsLp2hYkTYdWqMFbPe+/BaadB+/bw85+HmbhEJBKxDHydtC0BzZvD978PS5bAjBnh5O4VV8CBB4YPgBdegDwn3xGRwohl4OukbQmpUweGDw+hv2QJnHcePPZYaOrp3x/uvhs+/TTqKkVSIZaBLyWqSxe4+ebQ3DNxImzcCGecEY76f/rTsF1EikaBL7WvadMwJPM//hFO9A4YELpzdugAY8bA//6vmntEikCBL9Ex2zok89KlcOGF8NRTcMQR0KcP3Hkn/Oc/UVcpUjJiGfg6aZtCnTrB9deHZp3Jk8NVu2efHfr0jxsHy5dHXaFI4sUy8HXSNsUaNw5DMr/yShjG4cgj4X/+J1zF+7WvwaxZau4R2UOxDHwRzLYOyfzWW/CjH4WwHzIEevWC3/wGPv446ipFEkWBL/HXvn0YknnFijBkc506MHbs1ola3nkn6gpFEkGBL8mx995hSOYFC8IkLV/5Ctx4Y2j/P/FEePppNfeI7IQCX5LHLPTkeeCBMC3juHHw/PNwzDFhJM/bbgt9/EWkCgW+JFu7dqEP/4oVYUL2vfaC888PzT0XXRS6e4oIENPAV7dM2W2NGoUhmefOhb/9LQznMHFiuLp3xIjQv3/LlqirFIlULANf3TJlj5mFcXruvReWLQtDNsyfD8OGQbducOut8OGHUVcpEolYBr5IQRxwQBiSedkyuOeeMGHLhReG5p4LLgiDuYmkiAJfSl/DhluHZJ4zB0aPhl//OozdP2wYPPGEmnskFRT4ki7ZIZlXrIArr4SFC2HkSAacfnro4rluXdQVihSNAl/SqawstO8vWwbTp/NZixZw8cWhuef888NIniIlRoEv6Va/PowZw4KJE8PJ3a9/HaZOhR49Qr/+Rx6Bzz+PukqRglDgi2Rlh2ResQJ+8YtwUvfEE6G8HK67Dj74IOoKRWpEgS9S3X77hat33347XM174IFwySWhueekk2D8+DCSp7p3SsLUi7qA7TGzUcCo8vLyqEuRNKtXD04+OSx//3sYsuGZZ+Dhh8PtZqFv/4AB4WTwgAFhaIf69SMtW2RHYhn47v4Y8Fjfvn2/HXUtIkAYknny5PDzBx/Aiy+GZc6cMCn7nXeG2/baKzQN5X4ItG8fPhxEIhbLwBeJtZYtQ//9YcPCunto/sl+AMyZA5MmwQ03hNtbt94a/v37h2WffSIrX9JLgS9SU2ZhiOZOneCUU8K2TZtCH//sh8CLL8Ljj2/9nYMOqvoh0KsXNGgQTf2SGgp8kWKoXx8qK8Ny3nlh2/r1YXC37IfAn/8cLgKDEPa9e1dtCurcWU1BUlAKfJHa0rx56Nt/zDFh3T10Ac1+A5gzJ8zodcst4faWLbdtCtp33+jql8RT4ItExSx0+TzwwHDBF8DmzfDqq1XPBzz11NaZvDp1Ch8A2Q+B3r3D0NAieVDgi8RJvXqhPb9XL/h2ppPahg3hKuDsh8Bf/xqGf4bQdNSrF/TvT1mzZrD//mEOgDq6xEa2pcAXibumTWHw4LBk/fOfVZuC7rqLbhs3hovCmjeHfv2qng8oK4uqeokRBb5IEh1wQLjq96STwvrnn/PiXXfRH7Y2BY0fv3UcoPbtq54PqKwMk8JLqijwRUpB3bp83LFj+BbwrW+FbR9/DC+9VLVr6AMPfHF/evasej6gW7ewXUqWAl+kVO29Nxx+eFiy/vWvqlcJ33cf3H57uK1Jk9AUlPtNoE2baGqXolDgi6RJWRmMGhUWCDN9vfFG1W8BN9wQLhyDEPi5HwB9+4ZzCpJICnyRNKtTJ1z1e9BBcPrpYdsnn8DLL1ftGvrHP4bbzMJcAbkfAj17ht5FEnu19i6Z2WDgKuBVYLq7z6qt5xaR3dCoEQwcGJas99+v2hT0yCNhohgIA8ZVVlY9H3DggbpKOIbyCnwzmwqMBNa4e8+c7cOAm4G6wBR3H7+Th3FgI9AIWLnHFYtI7WvVCoYPDwuEC8Heeqtq19CJE+H668PtZWVbvwUMGBCagjRgXOTyPcKfBkwE7spuMLO6wCTgWEKAzzWzRwnhf2213z8b+Ku7P2tmZcANwGk1K11EImMWxvrp3Bm+8Y2w7bPPth0w7rHHtv5O165Vm4IOOUQDxtUy8+wl27u6o1kH4PHsEb6ZDQKucPfjMuvjANy9ethXf5wGwO/d/eQd3D4WGAtQVlZWOX369F3WtnHjRpo0aZLXfsRV0vdB9UcvjvtQb+NGmr72Gk1fe41mixfTbPFiGqxdC8CW+vXZ0KULG7p25cNu3Vh94IHUKS9PbFNQnF7/IUOGzHf3vtW316QNvw2wImd9JTBgR3c2s68CxwH7EL4tbJe73w7cDtC3b18fnHt14Q7MmjWLfO4XZ0nfB9Ufvdjuw8iRW392h+XL4cUXqTNnDs3nzKH5jBnw4IN0h9B0lNsU1K9f2JYAsX39c9TaSVt3fwh4KJ/7aopDkRJlFq76bd++6oBxixax5O67OWjdutAU9OSTWweMKy+v+iFw6KHQsGFUe5BoNQn8VUC7nPW2mW01pikORVKkXj049FDeXbeOg7JHyBs2wLx5W88HPPss/P734bb69UPo554PqKjQgHF5qEngzwUqzKwjIehPAb5RkKpEJN2aNoUhQ8KStWpV1WsDfvvbMJUkhB5A2TkDsh8CrVtHUnqc5dst815gMLCvma0Efubud5jZ94CnCD1zprr7q4UoSk06IrKNNm22GTCOxYurfghce+3WAeM6dKjaFNS7d+oHjMsr8N391B1snwHMKGhFqElHRPKQHQCuZ084++yw7aOPqg4YN2cO3H//1vsfckjVD4GuXVPVFKTroUWkdDRuDEccEZas1aurXiU8fTr8+tfhtqZNtx0w7oADoqm9FsQy8NWkIyIFs//+cMIJYYEwYNzrr1f9FjBhQugtBNC2bdXJYyorw0iiJSCWga8mHREpmjp1QlNO165wxhlh2yefwIIFVa8SfvDBrffPHTBuwADo3j2RA8Ylr2IRkUJr1AgGDQpL1nvvVW0K+uMf4Y47wm177x3GB8r9EMhz1IIoKfBFRLZn333h+OPDAiHQ33yzalPQLbeEMYSAQS1bhnMH2Q+Bvn3D/MIxEsvAVxu+iMSOWbjqt7y86oBxf/87vPgiax95hP1fey0MHZ29f9euVc8HHHxwuHAsIrEMfLXhi0giNGgQevn068drPXqw/+DBsHYtzJ279ZvAE0/AtGnh/o0aQZ8+VZuCOnSotQHjYhn4IiKJ1aIFDB0aFghNQcuWVZ07YPJkuOmmcPt++217lXCLFkUpTYEvIlJMZuEovkMHGDMmbNu0CRYtqno+YMaMrSd+KyrC0BG5J5ELIJaBrzZ8ESlp9euHoR5694Zzzw3bPvyw6oBxX/pSwZ82loGvNnwRSZ1mzeDoo8NSJOkZREJEJOUU+CIiKaHAFxFJCQW+iEhKxDLwzWyUmd2+fv36qEsRESkZsQx8d3/M3cc2j9k4FCIiSRbLwBcRkcJT4IuIpIR5jMdwNrN/A8vyuOu+wHtFLqfYkr4Pqj96Sd8H1V847d19v+obYx34+TKzee7eN+o6aiLp+6D6o5f0fVD9xacmHRGRlFDgi4ikRKkE/u1RF1AASd8H1R+9pO+D6i+ykmjDFxGRXSuVI3wREdkFBb6ISEokPvDNbJiZLTGzpWZ2adT17IqZtTOzmWb2DzN71cy+n9ne0sz+YmZvZP4tzqSWBWJmdc1sgZk9nlnvaGZzMu/DfWbWIOoad8bM9jGzP5jZa2a22MwGJek9MLOLMn8/i8zsXjNrFPf3wMymmtkaM1uUs227r7kFt2T2ZaGZ9Ymu8i9q3V7912X+hhaa2R/NbJ+c28Zl6l9iZsdFUnQ1iQ58M6sLTAKGA92BU82se7RV7dJm4Ifu3h0YCHw3U/OlwNPuXgE8nVmPs+8Di3PWfwnc6O7lwFrgnEiqyt/NwJPu3hXoRdiXRLwHZtYGuBDo6+49gbrAKcT/PZgGDKu2bUev+XCgIrOMBW6rpRp3Zhrb1v8XoKe7HwK8DowDyPyfPgXokfmdX2XyKlKJDnygP7DU3d9y98+A6cDoiGvaKXd/191fyvy8gRA0bQh1/zZzt98CJ0ZSYB7MrC0wApiSWTfgaOAPmbvEvf7mwJHAHQDu/pm7ryNB7wFhetK9zKwesDfwLjF/D9z9OeCDapt39JqPBu7y4AVgHzMr/CSvu2F79bv7n919c2b1BaBt5ufRwHR3/9Td3waWEvIqUkkP/DbAipz1lZltiWBmHYDewBygzN3fzdy0GiiLqq483ARcAmzJrLcC1uX84cf9fegI/Bu4M9MsNcXMGpOQ98DdVwETgOWEoF8PzCdZ70HWjl7zJP7fPhv4U+bnWNaf9MBPLDNrAjwI/MDdP8y9zUNf2Vj2lzWzkcAad58fdS01UA/oA9zm7r2Bj6jWfBPz96AF4QiyI3AA0JhtmxoSJ86v+a6Y2WWE5trfRV3LziQ98FcB7XLW22a2xZqZ1SeE/e/c/aHM5n9lv7Jm/l0TVX27cBhwgpm9Q2hCO5rQHr5PpnkB4v8+rARWuvuczPofCB8ASXkPjgHedvd/u/sm4CHC+5Kk9yBrR695Yv5vm9lZwEjgNN96YVMs60964M8FKjK9ExoQTpI8GnFNO5Vp774DWOzuN+Tc9ChwZubnM4FHaru2fLj7OHdv6+4dCK/3M+5+GjATODlzt9jWD+Duq4EVZnZQZtNXgH+QkPeA0JQz0Mz2zvw9ZetPzHuQY0ev+aPAGZneOgOB9TlNP7FhZsMIzZsnuPvHOTc9CpxiZg3NrCPh5POLUdRYhbsnegGOJ5wdfxO4LOp68qj3cMLX1oXAy5nleEI7+NPAG8D/B1pGXWse+zIYeDzzcyfCH/RS4AGgYdT17aL2Q4F5mffhYaBFkt4D4OfAa8Ai4G6gYdzfA+BewjmHTYRvWefs6DUHjNAD703gFUKPpDjWv5TQVp/9vzw55/6XZepfAgyPun5319AKIiJpkfQmHRERyZMCX0QkJRT4IiIpocAXEUkJBb6ISEoo8EVEUkKBLyKSEv8Hjft9dM33ijUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, y, color = 'r')\n",
    "plt.yscale('log')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdRklEQVR4nO3deXSV1bnH8e8TRhkuLSiTWFFRJBAGCSCTBlEIDoCICiJOCCLgUFsrILZaxYFqHRjkUlBELIMooyiiNUIVGRQRkItyrVW0ioptjRNS9/1jR5ubMiQ5J9nvOe/vs9ZZy/PmDE8Wr7/std/9Ptucc4iISPrLCF2AiIiUDwW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jERMXQBRzIoYce6ho3bhy6DImQL7/8kurVq4cuQ2SfonJ+vvrqq5865w4rejySgW9mZwFnNWnShA0bNoQuRyIkLy+PnJyc0GWI7FNUzk8z++u+jkdySsc5t9Q5N6xWrVqhSxERSRuRDHwREUk+Bb6ISEwo8EVEYkKBLyISEwp8EZGYSM/A37oV1qwJXYWISKSkZ+CPGwedO8O110J+fuhqREQiIZKBb2Znmdm0f/zjH6X7gFmzYMQIuP9+yMqClSuTW6CISAqKZOAnfONVzZowaRKsWgWVK0OPHnDZZfD558ktVEQkhUQy8JOma1fYtAnGjPGj/sxMePLJ0FWJiASR3oEPULUq3H47rF8PDRrAOedA//7w0UehKxMRKVfpH/g/aNMG1q6FO++EZcugWTOYORO0ibuIxER8Ah+gUiW44QY/zZOVBZdeCj17wl/+EroyEZEyF6/A/0HTppCXB1Om+PX6LVr4FT3/+lfoykREykw8Ax8gIwOuvNLfpHXyyX7Nfpcu8OaboSsTESkT8Q38H/zsZ/DUUzB7Nrz9tp/rv/VW2LMndGUiIkmlwAcwg0GD/Oi+Xz/49a8hO9uv7BERSRMK/MLq1oU5c2DxYvjsMzjxRLj+evjqq9CViYgkTIG/L717+9H+0KFw993QsiW88ELoqkREEqLA359atWDq1H8H/SmnwBVXQGn7+4iIBFaugW9mfc3sD2Y2z8x6lOd3l1pODrzxhp/amT7dt2dYsiR0VSIiJVbswDezh8xsl5ltKXI818y2m9kOMxt9oM9wzi1yzg0FhgPnl67kAKpVgwkT/J26depAnz4wYADs2hW6MhGRYivJCH8mkFv4gJlVACYDvYBMYKCZZZpZlpktK/KoW+it4wrel1qys2HDBr9sc+FC355h9my1ZxCRlFDswHfOrQJ2FzncHtjhnHvHObcHmAv0cc5tds6dWeSxy7y7gKedc68l79coR5Ur+w1WNm70d+wOHgxnnAHvvRe6MhGRA6qY4PsPB94v9Hwn0OEAr78KOBWoZWZNnHNTi77AzIYBwwDq1atHXl5egiWWoVtv5fBFizh6+nTc8cfzzrBhfNi7t7+LV8pEfn5+tM8JibWon5/mSjAdYWaNgWXOuRYFz/sDuc65ywueDwY6OOdGJaO47Oxst2HDhmR8VNl6910YNszvrNWli7+427Rp6KrSUl5eHjk5OaHLENmnqJyfZvaqcy676PFEh6IfAEcUet6o4FhCEt7isLw1bgwrVvh2y1u3QqtWvg3zd9+FrkxE5EeJBv564FgzO8rMKgMDgITXLCa8xWEIZnDxxf6GrbPO8rtstW/v5/pFRCKgJMsy5wBrgKZmttPMhjjn9gKjgBXANmC+c25r2ZSaIurXh8cfhyee8LtqtWvnw//rr0NXJiIxV5JVOgOdcw2cc5Wcc42cczMKji93zh3nnDvGOTc+GUWl3JTOvvTr50f7F1/sp3dat4bVq0NXJSIxFsnlJCk5pbMvP/0pzJjhL+bu2QMnnQQjR8I//xm6MhGJoUgGfto59VTYssVvsvLgg36HreXLQ1clIjETycBPiymdoqpXh3vvhZdfhpo1/c1agwfDp5+GrkxEYiKSgZ82Uzr7cuKJ8NprfpOVuXN9M7Z589SeQUTKXCQDP+1VqQK33OKDv3Fj34itb1/4IOFbGERE9kuBH1JWFqxZA/fc4y/sZmbCtGnw/fehKxORNBTJwE/LOfz9qVABrrsONm/23TivuAK6d4cdO0JXJiJpJpKBn9Zz+PtzzDHw3HO+D8/GjX70f/fdsHdv6MpEJE1EMvBjywyGDPE3bPXs6XfZ6tjR77glIpIgBX4UNWzoN1iZP9/32W/b1q/q+fbb0JWJSAqLZODHag5/f8zg3HP9aP+CC/wuWyecAK+8EroyEUlRkQz8WM7h70+dOvDII/7O3C++gE6d4Oc/hy+/DF2ZiKSYSAa+7EOvXr7X/ogRcN99/qLu88+HrkpEUogCP5XUrAmTJsGqVVCxou/Rc/nl8Pe/h65MRFKAAj8Vde0KmzbB6NF+l63MTFi8OHRVIhJxCvxUdcghcMcdsHYt1K3rWzOcfz58/HHoykQkoiIZ+FqlUwJt28L69TB+PCxa5Ef7s2erGZuI/IdIBr5W6ZRQpUowdiy8/jo0berbLp9xhl/DLyJSIJKBL6XUrJnfRvGBB/yF3ebN/YYrasYmIijw00+FCnDVVX6HrY4d/TLOnBx4663QlYlIYAr8dNW4MaxYAQ8/7DtxtmoFEyaoGZtIjCnw05kZXHKJb8/QqxfccAN06OCXdIpI7EQy8LVKJ8kaNIAnn4QFC/yuWtnZcNNNasYmEjORDHyt0ikj55zjR/uDBsFtt0GbNn5TdRGJhUgGvpSh2rX93bnPPOMbsHXpAtdcA/n5oSsTkTKmwI+rnj39Sp6RI/0yzqwsv6+uiKQtBX6c1awJEyf6tftVqkCPHnDZZfD556ErE5EyoMAXP63z+uswZgzMmuXbMyxcGLoqEUkyBb54VavC7bfDunVQvz706+d33Proo9CViUiSKPDl/zvhBB/6t98OS5f60f6sWWrGJpIGFPjynypV8tM7r7/u+/NcfLG/ceuvfw1dmYgkQIEv+3f88f6C7sSJ8Oc/Q4sWMHmymrGJpKhIBr7utI2QjAwYNcov4ezUyf/3ySfD9u2hKxOREopk4OtO2whq3NjfrDVzpt9MvVUruPNO+O670JWJSDFFMvAlosz8fP6bb8KZZ/p5/g4dYOPG0JWJSDEo8KXk6tf3jdgWLIAPP4R27eDGG+Gbb0JXJiIHoMCX0vuhGdvgwX4ZZ+vW8NJLoasSkf1Q4Etiatf2m6ysWOFH+F27wtVXqxmbSAQp8CU5evTwK3lGjYJJk/wSzmefDV2ViBSiwJfkqVHDd95cvdq3aujZEy69FHbvDl2ZiKDAl7LQubO/S3fsWHj0Ud+e4YknQlclEnsKfCkbVavC+PGwYQM0bAj9+/uHmrGJBKPAl7LVujWsXetv0lq2zI/2Z85UMzaRABT4UvYqVYIbboBNm6B5cz+vn5sL774bujKRWFHgS/lp2hRefNGv4nn5Zb+SZ+JENWMTKSflFvhm1szMpprZAjO7sry+VyImI8Pvo7tli99p6+qr4aST4H/+J3RlImmvWIFvZg+Z2S4z21LkeK6ZbTezHWY2+kCf4Zzb5pwbDpwHdC59yZIWjjwSnn4aHnnE363bqpW/W1fN2ETKTHFH+DOB3MIHzKwCMBnoBWQCA80s08yyzGxZkUfdgvf0Bp4CliftN5DUZQYXXQTbtkHv3r4fT/v2asYmUkaKFfjOuVVA0btn2gM7nHPvOOf2AHOBPs65zc65M4s8dhV8zhLnXC9gUDJ/CUlx9erB44/7tfoffeSbsY0Zo2ZsIklWMYH3Hg68X+j5TqDD/l5sZjlAP6AKBxjhm9kwYBhAvXr1yMvLS6BESSm1a1Nx2jSOmTKFBnfeyVePPcb266/nH1lZP74kPz9f54REVtTPz0QCv0Scc3lAXjFeNw2YBpCdne1ycnLKtC6JoLPOgpUrqTZsGG2uvtpf5L3jDqhZk7y8PHROSFRF/fxMZJXOB8ARhZ43KjiWMG1xKJx2Gmze7FfxTJnil3CuWBG6KpGUlkjgrweONbOjzKwyMABYkoyitMWhAL4Z2/33+w3Uq1WD3FyOv+MONWMTKaXiLsucA6wBmprZTjMb4pzbC4wCVgDbgPnOua1lV6rEVqdOfuXOjTdS9/nnoVkzv9uWiJRIcVfpDHTONXDOVXLONXLOzSg4vtw5d5xz7hjn3PhkFaUpHfkPVavCbbfx2tSp0KgRnHuu33Hrb38LXZlIyohkawVN6cj+5Ddp8u9mbE895ZuxPfywmrGJFEMkA1/kgCpW9M3Y3ngDsrLgssv8ZitqxiZyQJEMfE3pSLEcdxzk5cHkybBmjV/J88AD8K9/ha5MJJIiGfia0pFiy8iAESNg61a/gfo11/hmbNu2ha5MJHIiGfgiJfazn8Hy5TBrlu+82bq133FLzdhEfqTAl/RhBoMH++6bffvCuHG+L89rr4WuTCQSIhn4msOXhNSrB/PmwcKF8PHHvgPn6NHw9dehKxMJKpKBrzl8SYq+ff1o/5JL4K67/DTP6tWBixIJJ5KBL5I0P/0pTJ8OK1fCnj3+gu7IkfDFF6ErEyl3CnyJh1NP9dsqXnstPPig30z96adDVyVSriIZ+JrDlzJRvTrcey+89JJvzHb66X7Hrc8+C12ZSLmIZOBrDl/KVMeOvhnbTTfBnDm+PcPjj6s9g6S9SAa+SJmrUgV++1vYsAGOOALOOw/69VMzNklrCnyJt1at4JVXYMIEeOYZ33r5oYc02pe0pMAXqVgRrr8eNm3yfwCGDIEePeAvfwldmUhSRTLwddFWgjjuOHjhBb+KZ+1aNWOTtBPJwNdFWwkmIwOGD/fN2E4+2Tdj69pVzdgkLUQy8EWCO+IIv8HK7Nnw1ltqxiZpQYEvsj9mMGiQb89w9tm+GVt2Nrz6aujKREpFgS9yMHXrwty5sGgRfPIJdOjgd9xSMzZJMQp8keLq0+ffzdgmTPArelatCl2VSLEp8EVK4ic/8c3YnnsO9u71F3ZHjoR//jN0ZSIHFcnA17JMibzu3WHzZvj5z/0yzhYt1IxNIi+Sga9lmZISqleH3/8eXn4ZatZUMzaJvEgGvkhKOfFEv43ir3/tm7E1awbz56s9g0SOAl8kGapUgVtu8Us2jzwSzj/fN2P78MPQlYn8SIEvkkwtW8KaNfC73/lmbJmZMGOGRvsSCQp8kWSrWBF++Ut/Ubd1a7j8cjjtNHjnndCVScwp8EXKSpMm8Kc/wdSpsG4dZGXBffepGZsEo8AXKUsZGXDFFf6GrW7d/DLOLl38c5FypsAXKQ+NGsHSpfDYY/D229CmDdx6K+zZE7oyiREFvkh5MYMLLvCtlvv188s427Xz2yyKlINIBr7utJW0dthhfr3+4sXw6ae+GduvfqVmbFLmIhn4utNWYqF3bz+XP2SIX8bZsiW8+GLoqiSNRTLwRWKjVi2YNg2efx6+/x5ycuDKK9WMTcqEAl8kCk45xa/bv+46/wegeXNYvjx0VZJmFPgiUVGtGtxzj2/GVqsWnHEGXHihn+cXSQIFvkjUdOjgm7H95je+CVtmJsybp/YMkjAFvkgUVa4MN9/sm7E1bgwDBkDfvmrGJglR4ItEWVaWn+K5+2549lk/2p8+XaN9KRUFvkjUVawIv/jFv5uxDR0Kp56qZmxSYgp8kVTxQzO2//5vWL/eb6t4771qxibFpsAXSSUZGTBsmL9h65RT/DLOzp1h69bQlUkKUOCLpKIfmrH98Y/wv//rm7H99rdqxiYHpMAXSVVmMHCgH+337++XcWZn++kekX0o18A3s+pmtsHMzizP7xVJa4cd5kf6S5bA7t1+U/Xrr4evvgpdmURMsQLfzB4ys11mtqXI8Vwz225mO8xsdDE+6gZgfmkKFZGDOOssP5d/+eV+GWerVpCXF7oqiZDijvBnArmFD5hZBWAy0AvIBAaaWaaZZZnZsiKPumZ2GvAmsCuJ9YtIYbVq+VU8f/qTX6vfrRsMHw5qNS6AuWLewGFmjYFlzrkWBc87Ajc753oWPB8D4Jy7Yz/vHw9Ux/9x+Bo42zn3/T5eNwwYBlCvXr22c+fOLeGvJOksPz+fGjVqhC4jJWR88w1HPfwwjRYsYE/t2my/7jp2d+wYuqy0FpXzs1u3bq8657KLHq+YwGceDrxf6PlOoMP+XuycuxHAzC4BPt1X2Be8bhowDSA7O9vl5OQkUKKkm7y8PHROlEBuLqxbR5UhQ2g5dqzfceu++/y8vyRd1M/Pcl+l45yb6ZxbVt7fKxJb7dv7njw33wyPP+7bM8ydq/YMMZRI4H8AHFHoeaOCYwnTFociSVa5sl+2+dprcPTRfjlnnz7wQVL+l5UUkUjgrweONbOjzKwyMABYkoyitMWhSBlp0cI3Y7vnHnjuOT/a/8MfNNqPieIuy5wDrAGamtlOMxvinNsLjAJWANuA+c453d8tEnUVKviWDJs3Q9u2vlVD9+7+jl1Ja8UKfOfcQOdcA+dcJedcI+fcjILjy51zxznnjnHOjU9WUZrSESkHxxzj99KdNs3P8Wdlwe9/r2ZsaSySrRU0pSNSTsx8u+U33/Qtl3/xC+jUCbZsOfh7JeVEMvBFpJwdfjgsXgxz5vg++yecALfcomZsaSaSga8pHZEAzPxWitu2wbnn+mWcbdvCunWhK5MkiWTga0pHJKBDD4XHHvPtlz//HDp2hF/+Us3Y0kAkA19EIuDMM30ztqFD/TLOrCx44YXQVUkCIhn4mtIRiYhatWDqVB/0Zn6XrSuuUDO2FBXJwNeUjkjE5OTAG2/4qZ3p0/0NW0uXhq5KSiiSgS8iEVStGvzud/DKK1CnDvTu7ZuxffJJ6MqkmBT4IlIy7drBhg1+D90FC6BZM7/jltozRJ4CX0RKrnJluOkm2LgRmjSBQYP8iH/nztCVyQFEMvB10VYkRTRvDi+9BPfe63fZysz0O259v8/tLiSwSAa+LtqKpJAKFeDaa30ztvbt/ZaK3bvDjh2hK5MiIhn4IpKCjj4aVq70q3g2bvTr9u++G/buDV2ZFFDgi0jymMGQIb4ZW8+ecP31vhnb5s2hKxMU+CJSFho2hIULYd48ePdd34ztN7+Bb78NXVmsRTLwddFWJA2YwXnn+WZsAwb4ZZxt28LataEri61IBr4u2oqkkTp14NFH4amnfEuGjh39jltffhm6stiJZOCLSBo6/XTfjG34cL+Ms2VLv5RTyo0CX0TKz3/9F0yZAi++6Jdzdu/uu3H+/e+hK4sFBb6IlL+TToJNm+BXv4KHHvI3cC1ZErqqtKfAF5EwDjkE7rrLX8Q99FDo08df3N21K3RlaUuBLyJhZWf7Zmy33uqXcmZm+h231Iwt6SIZ+FqWKRIzlSrBuHH+Dt1jj4ULL/Q7br3/fujK0kokA1/LMkViKjMT/vxnuO8+yMvzc/sPPqhmbEkSycAXkRirUAGuuQa2bIEOHWDECOjWDd5+O3RlKU+BLyLRdNRR8OyzMGOGX9HTsiVMmKBmbAlQ4ItIdJnBZZf5Zmy5uXDDDXDiif4PgJSYAl9Eoq9hQ3jySZg/31/Izc72O26pGVuJKPBFJDWYwbnn+tH+BRfAbbdBmzawZk3oylKGAl9EUkudOvDII7B8OeTnQ+fOfsctNWM7KAW+iKSmXr18M7YRI+D++6FFC3juudBVRZoCX0RSV82aMGkSrFoFlSvDaaf5HbfUjG2fIhn4utNWREqka1d4/XUYPdpP92RmwqJFoauKnEgGvu60FZESO+QQuOMO34ytbl04+2y/49bHH4euLDIiGfgiIqXWti2sXw/jx8PixX60/+ijasaGAl9E0lGlSjB2rJ/madoULroIzjgD3nsvdGVBKfBFJH01awarV8MDD/gLu82b+x23YtqMTYEvIumtQgW46irfjK1jRxg5EnJy4K23QldW7hT4IhIPjRvDihXw8MOwebNvxnbXXbFqxqbAF5H4MINLLvHtGU4/3S/j7NAhNs3YFPgiEj8NGvhmbAsWwAcf+GZs48bBN9+ErqxMKfBFJL7OOceP9gcN8ss427SBl18OXVWZUeCLSLzVrg0zZ8Izz8BXX0GXLn7Hrfz80JUlnQJfRASgZ0+/kmfkSL+MMysLVq4MXVVSKfBFRH5QsyZMnOjX7lepAj16+B23Pv88dGVJUW6Bb2Y5ZrbazKaaWU55fa+ISIl16eLv0h0zBmbN8u0ZFi4MXVXCihX4ZvaQme0ysy1Fjuea2XYz22Fmow/yMQ7IB6oCO0tXrohIOalaFW6/Hdatg/r1oV8/v+PWRx+FrqzUijvCnwnkFj5gZhWAyUAvIBMYaGaZZpZlZsuKPOoCq51zvYAbgFuS9yuIiJShE07woX/77bB0qR/tz5qVks3YzBWzaDNrDCxzzrUoeN4RuNk517Pg+RgA59wdB/mcysAfnXP99/PzYcAwgHr16rWdO3du8X4TiYX8/Hxq1KgRugyJqWrvvUfTCROotXUru9u1Y/t11/Ft/fo//jwq52e3bt1edc5lFz1eMYHPPBx4v9DznUCH/b3YzPoBPYGfAJP29zrn3DRgGkB2drbLyclJoERJN3l5eeickKAuvBCmTKH26NF0HDoU7rwTrrwSMjIif36W20Vb59yTzrkrnHPnO+fyyut7RUSSKiMDRo3ySzg7dfL/fdJJsH176MoOKpHA/wA4otDzRgXHEqYtDkUk8ho39jdrzZzp79Zt1YqfPfYYfPdd6Mr2K5HAXw8ca2ZHFczLDwCWJKMobXEoIinBDC6+2Af+mWdy9PTpvhnbxo2hK9un4i7LnAOsAZqa2U4zG+Kc2wuMAlYA24D5zrmtyShKI3wRSSn168OCBWy5+Wb48ENo187vuBWxZmzFCnzn3EDnXAPnXCXnXCPn3IyC48udc8c5545xzo1PVlEa4YtIKvr05JP9aH/wYL+heuvW8NJLocv6kVoriIgkU+3afpOVFSv8CL9rV7/j1hdfhK5MgS8iUiZ69PAreUaNgsmToUUL/0cgoEgGvubwRSQt1KjhO2+uXg2HHAK5uX7Hrd27g5QTycDXHL6IpJXOnX0ztrFjYfZs357hiSfKvYxIBr6ISNqpWtXvqrVhAzRsCP37+x23/va3cishkoGvKR0RSVutW8PatX4Vz1NP+dH+zJnl0owtkoGvKR0RSWuVKsHo0bBpk7+Ye+mlfsetd98t06+NZOCLiMRC06bw4oswaRKsWePDf+JE+P77Mvk6Bb6ISEgZGX4f3S1b/E5bV1/t1+6/9VbyvyrpnygiIiV35JHw9NPwyCPwzjuwd2/SvyKSga+LtiISS2Zw0UV+Lj8zM+kfH8nA10VbEYm1KlXK5GMjGfgiIpJ8CnwRkZhQ4IuIxEQkA18XbUVEki+Sga+LtiIiyRfJwBcRkeRT4IuIxIS5cujQVlpm9gnw11K+vRYQ4iJAWXxvop9Z2veX5H3Ffe3BXnewnx8KfFrMmqIuxDkaxfOztJ8R4vw82Guicn4e6Zw77D+OOufS8gFMS5fvTfQzS/v+kryvuK892OuK8fMNIf5dy+IR4hyN4vlZ2s8IcX4e7DVRPz/TeUpnaRp9b6KfWdr3l+R9xX3twV4X6t8thBC/axTPz9J+Rojzs6TfGymRntIRKcrMNjjnskPXIbIvUT8/03mEL+lpWugCRA4g0uenRvgiIjGhEb6ISEwo8EVEYkKBLyISEwp8SWlmdrSZzTCzBaFrESnKzPqa2R/MbJ6Z9QhdjwJfIsfMHjKzXWa2pcjxXDPbbmY7zGw0gHPuHefckDCVShyV8Pxc5JwbCgwHzg9Rb2EKfImimUBu4QNmVgGYDPQCMoGBZpb8TT9FDm4mJT8/xxX8PCgFvkSOc24VsLvI4fbAjoIR/R5gLtCn3IuT2CvJ+WneXcDTzrnXyrvWohT4kioOB94v9HwncLiZ1TGzqUAbMxsTpjSRfZ+fwFXAqUB/MxseorDCKoYuQCQRzrnP8POjIpHjnHsAeCB0HT/QCF9SxQfAEYWeNyo4JhIFKXF+KvAlVawHjjWzo8ysMjAAWBK4JpEfpMT5qcCXyDGzOcAaoKmZ7TSzIc65vcAoYAWwDZjvnNsask6Jp1Q+P9U8TUQkJjTCFxGJCQW+iEhMKPBFRGJCgS8iEhMKfBGRmFDgi4jEhAJfRCQmFPgiIjGhwBcRiYn/A6139ENFB8/XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, y, color = 'r')\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "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
}