{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  概要:補間と近似
2  多項式補間(polynomial interpolation)
3  Lagrange(ラグランジュ) の内挿公式
3.0.1  python code
4  Newton(ニュートン) の差分商公式
4.1  Newton補間と多項式補間の一致の検証
5  数値積分 (Numerical integration)
5.1  中点則 (midpoint rule)
5.2  台形則 (trapezoidal rule)
5.3  Simpson(シンプソン)則
6  数値積分のコード
6.0.0.1  Midpoint rule(中点法)
6.0.0.2  Trapezoidal rule(台形公式)
6.0.0.3  Simpson's rule(シンプソンの公式)
7  課題
7.1  補間法
7.2  対数関数のニュートンの差分商補間(2014期末試験，25点)
7.2.1  差分商補間の表中の開いている箇所[ XXX ]を埋めよ．
7.2.2  ニュートンの二次多項式
7.2.3  ニュートンの三次多項式の値を求めよ．
7.3  数値積分(I)
8  解答例[7.3]
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "補間(interpolation)と数値積分(Integral)\n", "
\n", "
\n", "
\n", "file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/interpolationintegral\n", "
\n", "https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python\n", "
\n", "cc by Shigeto R. Nishitani 2017 \n", "
\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": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAEICAYAAAB2/gEGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFfNJREFUeJzt3X1sXXd9x/HPh6RQ3BR10M6iTWyzwRgV42HOeOpY46hM\noRQYCErDpRMTzBNapU6DFVCYAE0BDU0IELDJ46EMXEIlHlbaoqoodks3GCRQUEKoljEnDS0rpRTq\n3lGU9rs/zil1PDu+1/d7z7nHfr+kK3yOz/3dT46bLx+f+xBHhAAAANCbR9UdAAAAYC2gVAEAACSg\nVAEAACSgVAEAACSgVAEAACSgVAEAACSgVKFWtmdtv3GV9x2xPW97Q3YuAGvbauaH7Rfavq2fudBs\nlKo1yvac7Qs6OG7VpaZqi/9MEXE0IjZFxIN15gIw2MrZ8b9liZq3PS/p+ML5sdQstB22n/zwdkR8\nLSKeWnF8NAilCj3hKhGAhnhpWaIevt1RdyCsPZSqNc72623fYvsfbP/M9n/bfnH5vd2SXijpw+Vv\nbx8u9/+u7Rtt32P7NtsXL1jvStv/aPt62/dLmij3/VN5n/ts32R7dMF9XmD7W7Z/Xv7vC5bJ+tu2\n99r+qe27bU/bPqP83qcljUj6cpn1Cttj5W+SG8tjzrZ9TZn7sO0/X7D2u2xfbftfyowHbW9NP+EA\nGmHh/FhqFtq+uTz0u+W+19jeZvvYgjXmbL/F9vfK+fY526cu+P4Vtu+0fYftNy6+8oW1h1K1PjxX\n0m2SzpT0Pkkft+2I2CXpa5IuK39zu8z2aZJulHSVpN+UdImkj9o+d8F6r5W0W9Lpkm4p97Uk/V35\nGLdKmpYk24+XdJ2kD0l6gqT3S7rO9hOWyGlJ75V0tqSnSdoi6V2SFBGXSjqqR37bfN8S998j6Vh5\n/1dJeo/t7Qu+/7LymDMkXSPpwyc/bQDWg6VmYUT8UfntZ5b7PrfM3S+WtEPSkyQ9Q9LrJcn2Dkl/\nLekCSU+WtK1/fwIMCkrV+nAkIv65fO3ApyQ9UdLwMsdeJGkuIj4ZEccj4juSPi/p1QuO+deI+LeI\neCgiflnuuy4ibo6IByTtkvR821skvUTSf0bEp8v1PivpB5JeuviBI+JwRNwYEQ9ExE9UFLDzO/kD\nlo91nqS3RsQvI+JWSR+T9KcLDrslIq4vz8OnJT2zk7UBrAlfsn1veftS4rofiog7IuIeSV+W9Kxy\n/8WSPhkRByOirfIXRKxtG+sOgEr8+OEvIqJtW5I2LXPsqKTn2r53wb6NKkrIw25f4n6/3hcR87bv\nUXHF6GxJRxYde0TSOYsXsD0s6YMqLsOfrqL0/2yZnIudLemeiLhv0eMsfIrvxwu+bks61fbGiDje\n4WMAaK4/iYivPrxheyxp3cVz5ezy67Ml7VvwvaXmJtYYrlQhFm3fLummiDhjwW1TRLzpJPeRiqfq\nJEm2N0l6vKQ7ytvoomNHJP1oiTXeU679exHxOEmvU/GU4Mke92F3SHq87dM7eBwAWOxk82U17pS0\necH2luUOxNpBqcL/SPqtBdvXSvod25faPqW8/YHtp62wzoW2/9D2o1W8tuobEXG7pOvL9V5bviD0\nNZLOLR9nsdMlzUv6ue1zJP3NCll/rXysf5f0Xtun2n6GpDdI+swKuQFAWnq+LDtzOnC1pD+z/TTb\nQ5L+tpdwaAZKFT4o6VXlOwM/VD599scqXqB+h4pL238v6TErrHOVpHdKukfSuIqrTIqIn6p4ndab\nJf1U0hWSLoqIu5dY492Sfl/Sz1W8uP0Li77/XknvKF8T8ZYl7r9T0liZ+4uS3rnwcj8AnMQJs7Dc\n9y5JnypnzsXL3/X/i4ivqHiDzoykw5K+UX7rgaS8GECOyL7iifXG9pWSjkXEO+rOAgCDqLzaf0DS\nY3gd59rFlSoAAPrA9itsP8b2b6i44v9lCtXaRqkCAKA//kLSXZL+S9KDkt508sPRdDz9BwAAkIAr\nVQAAAAlq+fDPM888M8bGxjo69v7779dpp53W30B90tTs5K7Wesm9f//+uyPirD5GqkQ380taPz/f\nQUHuaq2X3B3Pr4io/DY+Ph6dmpmZ6fjYQdPU7OSu1nrJLWlf1DBvsm/dzK+I9fPzHRTkrtZ6yd3p\n/OLpPwAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgARppcr2BtvfsX1t1poAUAXmF4AMmVeqLpd0\nKHE9AKgK8wtAz1JKle3Nkl4i6WMZ6wFAVZhfALJkfaL6ByRdIen05Q6wPSlpUpKGh4c1Ozvb0cLz\n8/MdHztopqam6o6wKk095+SuVlNzL6Fv80tq7nliflWL3NXqW+5OPiH0ZDdJF0n6aPn1NknXrnSf\n9fKJ6sXpbZ6mnnNyV2stfKJ6v+dXRHN/vsyvapG7WoP8iernSXqZ7TlJeyRtt/2ZhHUBoN+YXwDS\n9FyqIuLtEbE5IsYkXSJpb0S8rudkANBnzC8AmficKgAAgASppSoiZiPiosw1gRVNT0tjYzp/+3Zp\nbKzYBrrE/ALQq6x3/wH1mJ6WJieldluWpCNHim1JarXqTAYAWGd4+g/NtmuX1G6fuK/dLvYDAFAh\nShWa7ejR7vYDANAnlCo028hId/sBAOgTShWabfduaWjoxH1DQ8V+AAAqRKlCs7Va0tSUNDqqsKXR\n0WKbF6kDACpGqULztVrS3Jxu2rtXmpujUAEAakGpAgAASECpAgAASECpAgAASECpAgAASECpAgAA\nSECpAgAASECpAgAASECpAgAASECpAgAASECpAgAASECpAgAASECpAgAASECpAgAASECpAgAASECp\nAgAASECpAgAASECpwiOmp6WxMZ2/fbs0NlZso3843wCaivm1pI11B8CAmJ6WJieldluWpCNHim1J\narXqTLY2cb4BNBXza1lcqUJh1y6p3T5xX7td7Ec+zjeApmJ+LYtShcLRo93tR2843wCaivm1LEoV\nCiMj3e1HbzjfAJqK+bUsShUKu3dLQ0Mn7hsaKvYjH+cbQFMxv5ZFqUKh1ZKmpqTRUYUtjY4W2+v8\nRYd9w/kG0FTMr2VRqvCIVkuam9NNe/dKc3P8Bek3zjeApmJ+LannUmX7VNvftP1d2wdtvzsjGABU\ngRkGIEvG51Q9IGl7RMzbPkXSLba/EhHfSFgbAPqNGQYgRc+lKiJC0ny5eUp5i17XBYAqMMMAZHEx\nT3pcxN4gab+kJ0v6SES8dYljJiVNStLw8PD4nj17Olp7fn5emzZt6jljHSYmJjQzM1N3jK419ZyT\nu1rd5p6YmNgfEVv7GGnVVpphq51fUnN/vsyvapG7Wn2bXxGRdpN0hqQZSU8/2XHj4+PRqZmZmY6P\nHTTF6W2epp5zcler29yS9kXivOnHrZMZ1s38imjuz5f5VS1yV6tf8yv13X8RcW85kHZkrgsAVWCG\nAehFxrv/zrJ9Rvn1YyW9SNIPel0XAKrADAOQJePdf0+U9KnyNQmPknR1RFybsC4AVIEZBiBFxrv/\nvifp2QlZAKByzDAAWfhEdQAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASU\nKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAA\ngASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASUKgAAgASU\nKgAAgASUKgAAgASUKgAAgASUKgAAgAQ9lyrbW2zP2P6+7YO2L88IBgBVYIYByLIxYY3jkt4cEd+2\nfbqk/bZvjIjvJ6wNAP3GDAOQoucrVRFxZ0R8u/z6PkmHJJ3T67oAUAVmGIAsjoi8xewxSTdLenpE\n/GLR9yYlTUrS8PDw+J49ezpac35+Xps2bUrLWKWJiQnNzMzUHaNrTT3n5K5Wt7knJib2R8TWPkbq\n2XIzbLXzS2ruz5f5VS1yV6tv8ysiUm6SNknaL+mVKx07Pj4enZqZmen42EFTnN7maeo5J3e1us0t\naV8kzZt+3DqdYd3Mr4jm/nyZX9Uid7X6Nb9S3v1n+xRJn5c0HRFfyFgTAKrCDAOQIePdf5b0cUmH\nIuL9vUcCgOowwwBkybhSdZ6kSyVtt31rebswYV0AqAIzDECKnj9SISJukeSELABQOWYYajU9Le3a\npfOPHpVGRqTdu6VWq+5UWKWMz6kCAADdmp6WJieldrto9UeOFNsSxaqh+GdqAACow65dUrt94r52\nu9iPRqJUAQBQh6NHu9uPgUepAgCgDiMj3e3HwKNUAQBQh927paGhE/cNDRX70UiUKgAA6tBqSVNT\n0uiowpZGR4ttXqTeWJQqAADq0mpJc3O6ae9eaW6OQtVwlCoAAIAElCoAAIAElCoAAIAElCoAAIAE\nlCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoA\nAIAElCoAAIAElCoAAIAElKp+mJ6Wxsb0oCSNjRXbAID+Kefu+du3M3dRm411B1hzpqelyUmp3S4a\n65EjxbYktVp1JgOAtWnB3LXE3EVtuFKVbdcuqd0+cV+7XewHAORj7mJAUKqyHT3a3X4AQG+YuxgQ\nlKpsIyPd7QcA9Ia5iwFBqcq2e7c0NHTivqGhYj8AIB9zFwOCUpWt1ZKmpqTRUT0kSaOjxTYvlgSA\n/lgwd8Nm7qI2vPuvH1otqdXSBlsxN1d3GgBY+8q5e9PsrLZt21Z3GqxTKVeqbH/C9l22D2SsBwBV\nYX4ByJL19N+VknYkrQUAVbpSzC8ACVJKVUTcLOmejLUAoErMLwBZKntNle1JSZOSNDw8rNnZ2Y7u\nNz8/3/Gxg6iJ2Zt6zsldrabmXo3Vzi+p2eepibmber7JXa2+5Y6IlJukMUkHOjl2fHw8OjUzM9Px\nsYOmOL3N09RzTu5qdZtb0r5ImjfZt37Nr4jm/nyZX9Uid7X6Nb/4SAUAAIAElCoAAIAEWR+p8FlJ\nX5f0VNvHbL8hY10A6DfmF4AsKS9Uj4idGesAQNWYXwCy8PQfAABAAkoVAABAAkoVAABAAkoVAABA\nAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoV\nAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAABAAkoVAOAR\n09PS2JgelKSxsWIbQEc21h0AADAgpqelyUmp3S5+4z5ypNiWpFarzmRAI3ClCgBQ2LVLardP3Ndu\nF/sBrIhSBQAoHD3a3X4AJ6BUAQAKIyPd7QdwAkoVAKCwe7c0NHTivqGhYj+AFVGqAACFVkuampJG\nR/WQJI2OFtu8SB3oCO/+AwA8otWSWi1tsBVzc3WnARol5UqV7R22b7N92PbbMtYEgKowwwBk6LlU\n2d4g6SOSXizpXEk7bZ/b67oAUAVmGIAsGVeqniPpcET8MCJ+JWmPpJcnrAsAVWCGAUiR8ZqqcyTd\nvmD7mKTnLj7I9qSkSUkaHh7W7OxsR4tPTU1pYmKi95Q1sV13BKAvdu7cWXeELCvOsNXOL6nZM4z5\nhbWqX/OrsheqR8SUpClJ2rp1a2zbtq3j+1511VV9StVfthURdcfo2uzsrLr5+QwKclerqblXo5f5\nJTVzhjG/qkXuavUrd8bTfz+StGXB9uZyHwA0ATMMQIqMUvUtSU+x/STbj5Z0iaRrEtYFgCowwwCk\n6Pnpv4g4bvsySTdI2iDpExFxsOdkAFABZhiALCmvqYqI6yVdn7EWAFSNGQYgA/9MDQAAQAJKFQAA\nQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJK\nFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAA\nQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJKFQAAQAJK\nFQAAQIKeSpXtV9s+aPsh21uzQgFAFZhhADL1eqXqgKRXSro5IQsAVI0ZBiDNxl7uHBGHJMl2ThoA\nqBAzDECmnkpVN2xPSpqUpOHhYc3OznZ0v/n5+Y6PHTQ7d+5sZPamnnNyV6upuVdjtfNLau55Yn5V\ni9zV6lvuiDjpTdJXVVwiX3x7+YJjZiVtXWmth2/j4+PRqZmZmY6PHTRNzU7uaq2X3JL2RYczIvOW\nPcO6mV+rOU+DgtzVIne1+jW/VrxSFREX5FU4AKgWMwxAVfhIBQAAgAS9fqTCK2wfk/R8SdfZviEn\nFgD0HzMMQKZe3/33RUlfTMoCAJVihgHIxNN/AAAACShVAAAACShVAAAACShVAAAACVx8plXFD2r/\nRNKRDg8/U9LdfYzTT03NTu5qrZfcoxFxVr/CVKXL+SWtn5/voCB3tdZL7o7mVy2lqhu290VEI//1\n+KZmJ3e1yL22NfU8kbta5K5Wv3Lz9B8AAEACShUAAECCJpSqqboD9KCp2cldLXKvbU09T+SuFrmr\n1ZfcA/+aKgAAgCZowpUqAACAgUepAgAASNCIUmX71bYP2n7I9sC/ddP2Dtu32T5s+2115+mU7U/Y\nvsv2gbqzdMr2Ftsztr9f/jdyed2ZOmH7VNvftP3dMve7687UDdsbbH/H9rV1Z2mCJs0w5le1mGH1\n6NcMa0SpknRA0isl3Vx3kJXY3iDpI5JeLOlcSTttn1tvqo5dKWlH3SG6dFzSmyPiXEnPk/SXDTnf\nD0jaHhHPlPQsSTtsP6/mTN24XNKhukM0SCNmGPOrFsywevRlhjWiVEXEoYi4re4cHXqOpMMR8cOI\n+JWkPZJeXnOmjkTEzZLuqTtHNyLizoj4dvn1fSr+kpxTb6qVRWG+3DylvDXiXSO2N0t6iaSP1Z2l\nKRo0w5hfFWOGVa+fM6wRpaphzpF0+4LtY2rAX5C1wPaYpGdL+o96k3SmvPx8q6S7JN0YEY3ILekD\nkq6Q9FDdQZCO+VUjZlhl+jbDBqZU2f6q7QNL3BrxWxLqZXuTpM9L+quI+EXdeToREQ9GxLMkbZb0\nHNtPrzvTSmxfJOmuiNhfd5ZBwwxDL5hh1ej3DNvYj0VXIyIuqDtDkh9J2rJge3O5D31i+xQVw2g6\nIr5Qd55uRcS9tmdUvB5k0F9ke56kl9m+UNKpkh5n+zMR8bqac9Vujcww5lcNmGGV6usMG5grVWvI\ntyQ9xfaTbD9a0iWSrqk505pl25I+LulQRLy/7jydsn2W7TPKrx8r6UWSflBvqpVFxNsjYnNEjKn4\nb3svhWpNYX5VjBlWrX7PsEaUKtuvsH1M0vMlXWf7hrozLScijku6TNINKl5weHVEHKw3VWdsf1bS\n1yU91fYx22+oO1MHzpN0qaTttm8tbxfWHaoDT5Q0Y/t7Kv6P7MaI4OMJ1qimzDDmVy2YYWsI/0wN\nAABAgkZcqQIAABh0lCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAElCoAAIAE/wdShHKGJrUA\nVgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "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", "\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", "となる．$a_i$と$y_i$をそれぞれベクトルとみなすと\n", "\n", " \n", "\n", "Aとbからa_iを導出：\n", "\n", " \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", " \n", "\n", "$L_0(x)$：\n", "\n", " \n", "\n", " \n", "\n", "$L_1(x)$：\n", "\n", " \n", "\n", "\n", " \n", "\n", "$L_2(x)$：\n", "\n", " \n", "\n", " \n", "\n", "$F(x)$：\n", "\n", " \n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### python code\n", "pythonではLagrange補間はinterpolate.lagrangeで用意されている．" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 3 2\n", "-1 x + 3 x - 1 x + 1\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9xvHvgjDHKggEVAjqRRBFxESkIpogReAiICDW\nG5VWaRzqUG2t1Dgr9arVWqqPSutYUaSKoIiCQoCCQgEVZ7wOTBULCA4Bw7juH79YFYGc5Oxz9tn7\nvJ/nOU+Gs88+v8VOXnbWXnst571HRETio07YBYiISLAU7CIiMaNgFxGJGQW7iEjMKNhFRGJGwS4i\nEjMKdhGRmFGwi4jEjIJdRCRmcsJ40+bNm/t27drV6rUbN26kSZMmwRYUErUl88SlHaC2ZKpk2rJ4\n8eJ13vsW1W0XSrC3a9eORYsW1eq1s2bNoqioKNiCQqK2ZJ64tAPUlkyVTFucc8sT2U5dMSIiMaNg\nFxGJGQW7iEjMKNhFRGJGwS4iEjMKdhGRmFGwi4jETCjj2EWy2eefw7JlsHatPT77DCorYcsWe9Sv\nDw0bQqNGsO++kJcHrVpBfr59T6Q6SQe7c64hMAdoULW/J7331ya7X5Go+/JLWLgQliyxx1tvwUcf\nWbDXVtu20KEDdOkCxxwD3bvDAQcEV7PEQxBn7JuBXt77CudcPWCuc+557/38APYtEhlffAEzZsDM\nmTB3Lrz5JuzYYc+1agWdO1sQH3ggtGtnZ+ItWthZeaNG0KAB5OTAtm12Br9pk53Nf/oprF4NH34I\nS5fCe+/BmDF2dg+2v7597dGrF+ROHgdlZZywYoX9TzB6NJSUhPbvIumXdLB77z1QUfVlvaqHT3a/\nIlHwf/8HEyfCs8/C/PmwfTvk5lqAX3MN/PjHcOSR0LJl4vusV88ee+1l4d+p0w+32bzZ/gqYP9/+\nM3nkEbjnHmhUfxsDt9fn9O1d6MtqGixfDqWl9iKFe9YIpI/dOVcXWAz8F3C3935BEPsVyUTLlsGj\nj8KECXZWDnDUUTBqFJx0koV6vXqpraFBA+jWzR4XX2xn73PnwlODHmNCRT+e4FSas5ZzuJ/zNt1L\nu7IyBXsWcXbCHdDOnNsHeBq4yHv/1k7PlQKlAHl5eQXjx4+v1XtUVFSQm5ubbKkZQW3JPLtrR2Vl\nHcrLW/LCC6144419AOjc+XOOP34tPXuuIy9vc7pL3aUTevVim6/LDE7kL/yCyQxiB3UYyDP0vacN\nHTt+FXaJtRKXny9Iri3FxcWLvfeF1W7ovQ/0AVwD/GZP2xQUFPjaKi8vr/VrM43aknl2bsc773h/\n0UXe77239+D9IYd4f9NN3n/8cSjlVS8/3wqteqxkf1/Gjb5pnQ0evO/Xz/tXXgm7yJqLy8+X98m1\nBVjkE8jhpMexO+daVJ2p45xrBPwEeC/Z/YqExXuYPt0uRnbqBPfdB/37w+zZduGyrMwufmak0aOh\nceP/fHkA/+Kmxjez7N5p3HyzjdL58Y9h+HAboSPxFMQNSq2BcufcG8BC4EXv/ZQA9iuSVtu2wYwZ\nLenSxfrKlyyBG2+EVavgscfg+OPBubCrrEZJCYwdC/n5eOds8PvYsfzoF6cxapRdH7j+enjuOTj0\nUPjtb2HjxrCLlqAlHeze+ze8912990d47w/33t8QRGEi6bJ1K9x/P3TsCDfd1Int2+GhhywEr7rK\nhiRGSkkJLFvG7JkzrRHfuWjapImN1nn/fTj9dLjtNjj8cJg2LbxyJXiaUkCy1rZt8PDDFugjR8I+\n+8ANN7zFm2/CiBE28iSu9t/f/vOaM8fa2bcvnHlmcjdPSeZQsEvW8d7Gnh9+OPzsZ7D33jBlivU/\n9+y5jjpZ9FvRsye8/jpcfTWMH293tP7jH2FXJcnKoh9hERvr3aMHDB1q/eVPPQWLF8N//3cE+s9T\npGFDuOEGmDfPxt8XFcG119rNVhJNCnbJCsuW2UiQnj3t87Fj7eaiIUOyN9B31q0bvPaadcnccAP0\n62dTGkj0KNgl1jZutG6Gjh2tu+Xaa20agF/8wuZlke/bay/re//rX21459FH2+ggiRYFu8TSN/3o\nhx4KN91kXS9Ll8J119nIENmzc86xvvYtW2zc+6RJYVckNaFgl9j58EO7oWjoUGja1AJq3Dho0ybs\nyqKlWze7/tCli3VZ/fnPYVckiVKwS2xs2QK//72Ndpk3D/70Jwum444Lu7Loysuz2SMHDrTJxn7z\nm2+nIpbMpV5GiYVXXrF+87ffhmHDLNT32y/squKhcWMbPXTJJXD77bBund3QVbdu2JXJ7ijYJdIq\nKuzu0DFjbCWhZ5+FAQPCrip+6ta1rpiWLe0CdGUl/O1vqZ+eWGpHwS6RNXOmXeRbtgx++Uu4+WYb\n1SGp4ZxNR9Cokc0xs3mz3dQU5zt0o0p97BI5FRUW5CeeaGeMc+bAXXcp1NPl8svt7H3SJDjtNJtr\nRzKLgl0iZfZsOOIIWwbu0kvtdviePcOuKvtceKGF++TJNi2D7lLNLOqKkUiorLS+9DvugIMOsoBX\noIfrwgvtr6ff/c4usI4dq7t4M4WCXTLe66/DGWfYiJfzzrOpZmOySlrkjRpl4T56NDRrBrfcEnZF\nAgp2yWA7dtjwurIyaN4cpk61+Usks9x4I6xfD7feCm3b2vUPCZeCXTLSypVw1lkwa5bd9Th2LOy7\nb9hVya44Z8NNV62ym5jatLEbmiQ8ungqGWfiRLuNfeFCuxHmyScV6pkuJwcefxyOOgp++lP45z/D\nrii7KdglY2zaBKWlNsfLwQfbFLJnn60LclHRpInNoJmXB4MHwyefhF1R9lKwS0Z4800oLLTpYq+4\nwuZ6ad8+7KqkpvLy4Jln4MsvrQutsjLsirKTgl1C5T3ce6/NJLhhA0yfDv/7v1C/ftiVSW117gyP\nPAILFtgoJu/Drij7KNglNJ9/bqsanX8+nHCCLejQu3fYVUkQhgyxOWUefljT/YZBwS6hWLjQLrRN\nmmTD5KZOtQmmJD6uuQYGDYJf/xrmzw+7muyiYJe08h7uvNMWlN6+3eZ5ufxyqKOfxNipUwcefNBm\n3TztNBvrLumhXydJmw0b7E/0Sy+1FY5ee82WXZP4atoU/v53+PRTGDFCi3Ski4Jd0uKbrpcpU+CP\nf4Snn7Zb0CX+CgvtDuIpU+APfwi7muygYJeU8t6m1O3Rw87W5s6FX/1KY9OzzS9/afcnlJXZcoWS\nWgp2SZkvv7S+1YsugpNOsq6XY44JuyoJg3M2LUReHpSU2M1okjoKdkmJN96wP8EnTrRRL5Mnq+sl\n2zVrZsMfly61RbEldRTsErgHH7Qz84oKKC/XqBf51okn2vDHe+6xPndJjaR/3ZxzbZxz5c65d5xz\nbzvnLgmiMImer7+2NUjPPhuOPda6XrQYhuxs9GhbBWvkSPjss7CriacgzqO2Ab/23ncCugO/dM51\nCmC/EgXjxkG7duxXXEr3vd/lgQdspaPp060/VWRnDRrYlAOffWYX0iV4SQe793619/7Vqs+/At4F\n9k92vxIB48ZBaSlPLS+gkEWs2tqSqQ1O4caO46hbN+ziJJN16QJXXgmPPqoumVQItOfTOdcO6Aos\nCHK/kpm2Xnktl226kWE8xaG8y2t0pd/mSTamTaQaZWU2Ydi559q8QRIc5wOaes05lwvMBkZ77yfu\n4vlSoBQgLy+vYPz48bV6n4qKCnJjsuBllNuydm0D7hy+g5fpwUWM4Q/8hvpsBcA7x+yZM0OusHai\nfEx2FoW2LF2aywUXFNCnz6dcccXS3W4XhbYkKpm2FBcXL/beF1a7ofc+6QdQD5gGXJbI9gUFBb62\nysvLa/3aTBPVtkyf7n3z5t7nuq/8E5zqvd2H9O0jPz/sEmstqsdkV6LSllGj7MdmT+VGpS2JSKYt\nwCKfQMYGMSrGAfcD73rv70h2f5K5tm+H666zm41atYJFt8xkeOPnvr9R48Y27EEkQVdfDQceaNM3\nb94cdjXxEEQfew/gTKCXc+71qkf/APYrGWTtWpu46/rr4cwzbRGFDpcPtNsJ8/PxzkF+vn1dUhJ2\nuRIhjRvD3XfDe+9pLpmg5CS7A+/9XEAzf8TY3Lm2QPG6dfCXv9hY9f/M9VJSAiUlzJ41i6KiojDL\nlAjr1w9OPRVuusl+1g4+OOyKok33A8pueW9nUEVF0KiRLZYwcqQm8JLUuPNOqFcPLrhAy+klS8Eu\nu7R+va1+c/nltuL8okVw5JFhVyVxtt9+dsY+fbrNLSS1p2CXH1iwwOZOf+EF+NOfbKGEvfcOuyrJ\nBuefD4cdBpddBpWVYVcTXQp2+Q/vbRGM446z7pZ58+Dii9X1IulTr56dTHz8MdyhMXa1pmAXwObt\nGDjQzpQGDIBXX4Wjjw67KslGJ54Ip5xio2ZXrQq7mmhSsAvz5kHXrta3OWaMzaHetGnYVUk2u/12\nu2/iiivCriSaFOxZbPt2Oys64QT7E/jll221I3W9SNgOPNAW43jsMbvmIzWjYM9Sq1fbHaRXXQXD\nh9vc6QUFYVcl8q0rroCWLW1kloY/1oyCPQtNmWILHbzyCjzwgM2++6MfhV2VyPfttZdNYfGPf8C8\nefuGXU6kKNizSGWljXI5+WQ44ABbLf7nP1fXi2SukSOhQwcYO/Zgtm4Nu5roULBniTfesFEuf/6z\nrVozfz507Bh2VSJ7Vq8e3HILrFzZmL/+NexqokPBHnM7dtjY9KOPtom8pk61rxs0CLsykcQMHAhH\nHPE5111nC6RL9RTsMbZypV0gvewy6NsX3nzTJlsSiRLn4NxzP2TNGrt5SaqnYI8h720tyc6d7QLp\nfffBpEnQokXYlYnUTqdOX3HyyXDbbbBhQ9jVZD4Fe8ysWWPTn555Jhx+OCxZAqWlukAq0XfjjfDF\nF5qzPREK9hh58kmbQOnZZ+Hmm2H2bM1rLfHRpQucdpp1x6xZE3Y1mU3BHgNr1tjiBKeeaosYvfoq\njBoFdeuGXZlIsK6/Hr7+2k5cZPcU7BHmvd1c1KkTPP20zWX9yit21i4SRx06wIgRcM89miBsTxTs\nEbVsmc3CeMYZ0L69TQlQVmbjfkXi7JprbJ6jW28Nu5LMpWCPmG3b7OLRYYdZH/qdd9qapJ06hV2Z\nSHq0awdnnWXrpq9eHXY1mUnBHiEvvwyFhTYp0oknwjvvwCWXqC9dss+VV9pJzm23hV1JZlKwR8Ca\nNTanS48etiDGU0/ZmpBt24ZdmUg4Dj4YSkrg3nvh3/8Ou5rMo2DPYFu22O3/hxxiF0lHjYL33oMh\nQzQuXeTKK2HzZluUQ75PwZ6BvLepdTt3tukAune3SbxuvhmaNAm7OpHM0KGDDfO9+26bB0m+pWDP\nMAsXQq9eNrWuc/Dcc/DCC5qJUWRXyspg0yabtVS+pWDPEO+/b2cf3brBW2/Z2qNvvgn9+4ddmUjm\n6tQJBg2Cu+7SzI/fpWAP2ccf24XRQw+1qQCuugo+/NDWHtWYdJHqjRplE4NpvvZvKdhD8v77cOut\nHTjkEHj8cRu2+NFHNtGRlqkTSVz37rYg++2324ADUbCn3aJFNpFRx44wY0ZLzj/fztDvuAPy8sKu\nTiSarrjCphh4/PGwK8kMCvY02L7d5kM//nhbyej55+0H8fHH5zNmDOy/f9gVikRb3762QPstt9iq\nYdkukGB3zj3gnFvjnHsriP3t0rhx0K4dJ/TqZfcUjxuXsrcKyurVNjHXQQfBKafAihX25+LKlTZ0\nsVkzrc4rEgTn7GTp3XdtJFm2C+qM/SGgb0D7+qFx42y1iOXLcd7D8uX2dQaGe2WlzYs+aJDdGXr1\n1XaD0ZNPwgcf2Lj0vfcOu0qR+Bk+3H7n7rgj7ErCF0iwe+/nAOuD2NcuVQ1WXUA3nmUAW6hng1fL\nylL2ljWxebOdJZx9NrRubfOiL1wIl15qF0lffBGGDoWcnLArFYmvnBwbTTZrls12ms2c9z6YHTnX\nDpjivT98N8+XAqUAeXl5BePHj0943yf06oXznhE8xCOMoCnrGcpTDGcCftpV1K8fTBtqYu3aBvzz\nn01ZuLAZixY1Y+PGHJo02caxx66jT59/07Xrhmon56qoqCA3Nzc9BadYXNoSl3ZAdraloiKH4cO7\nc9xx67jyyvfSUFnNJXNciouLF3vvC6vd0HsfyANoB7yVyLYFBQW+RvLzvQe/hRz/HP38GTzim/CV\nB++bNPH+5JO9HzPG+4ULvd+8uWa7TsTWrd6/8473Dzzg/TnneN+hg/d247/3++3n/ciR3k+dWvP3\nLi8vD77YkMSlLXFph/fZ25aLL/Y+J8f7VatSV08ykjkuwCKfQMZGo3Ng9GgoLaXepk3053n68zyb\nGu3LS+c9ybTNRTz/vN3cA9CwIRx5pN3w06GD9W+3bg2tWkHLltCo0Q8n0Nqxw25wWLvWHitW2Jjy\nDz+Et9+2O0ErK23bZs3g2GNh5Ei7En/YYZqQSySTXHKJTTFw993w+9+HXU04ohHsJSX2sawMv2IF\nrm1bGo8ezcCSIgZi584rVsCCBfZYvBimToUHH9z17ho2hAYNbD7nLVtg624Gp7RubbcsX3CBLaRb\nWGjjz+tokKhIxjroIBg8GO67zy7DZePEeYEEu3PucaAIaO6cWwVc672/P4h9/0dJCZSUMHvWLIqK\ninZ6f1vEOT/frox/44svbCTKp5/anM1r1thCuJWV9sjJgfr1LeSbNoUWLaB5c7uynp9vZ/ciEj2X\nXWbrAP/tb3DeeWFXk36BBLv3/vQg9hO0vfeGgoKwqxCRdOvRA7p2tcnBzj03+7pL1akgIrHjHFx4\noV0jmz077GrST8EuIrF0+uk22OGuu8KuJP0U7CISS40a2ei1SZNsGo9somAXkdg6/3wbznzffWFX\nkl4KdhGJrXbtbJnJsWNt6o9soWAXkVi78EK78XDChLArSR8Fu4jEWu/edgf6vfeGXUn6KNhFJNac\ns1m+X37ZpgfJBgp2EYm9ESPsLvNsuYiqYBeR2GveHIYNsykGNm0Ku5rUU7CLSFY491ybP+qJJ8Ku\nJPUU7CKSFXr2tOm8s6E7RsEuIlnhm4uoCxbAkiVhV5NaCnYRyRpnnWXrMYwdG3YlqaVgF5Gs0awZ\nDBkCjz1mazPElYJdRLLKOefA55/b5GBxpWAXkaxSVAQHHgj3B7vGW0ZRsItIVqlTB37+c5gxAz7+\nOOxqUkPBLiJZZ8QIGyXz0ENhV5IaCnYRyTpt28JPfgIPPgjbt4ddTfAU7CKSlc45x1ZWmjEj7EqC\np2AXkaw0aJANf3zggbArCZ6CXUSyUoMGtuD15Mk2h0ycKNhFJGuddRZUVsKTT4ZdSbAU7CKStY4+\nGjp0gEceCbuSYCnYRSRrOWdn7XPmxGtMu4JdRLLaGWfYx0cfDbeOICnYRSSrtW0LxcXWHeN92NUE\nQ8EuIlnvrLPggw9g/vywKwmGgl1Est7QodCoUXwuogYS7M65vs65pc65D5xzo4LYp4hIuuy1Fwwe\nDBMmwJYtYVeTvKSD3TlXF7gb6Ad0Ak53znVKdr8iIulUUgLr18OLL4ZdSfKCOGPvBnzgvf/Ie78F\nGA8MCmC/IiJp06cP7Luvra4Udc4neRnYOTcM6Ou9H1n19ZnAMd77C3farhQoBcjLyysYP358rd6v\noqKC3NzcpGrOFGpL5olLO0BtqY0//rE906e3YuLEeTRqtCMl75FMW4qLixd77wur2y6nVnuvBe/9\nWGAsQGFhoS8qKqrVfmbNmkVtX5tp1JbME5d2gNpSG3XrwjPPwIYNx9OvX2reIx1tCaIr5l9Am+98\nfUDV90REIqVHD2jTJvrdMUEE+0KgvXPuQOdcfeCnwDMB7FdEJK3q1LEZH6dNg3Xrwq6m9pIOdu/9\nNuBCYBrwLjDBe/92svsVEQnD//wPbNsW7RkfAxnH7r2f6r0/xHt/sPd+dBD7FBEJwxFHQKdO0e6O\n0Z2nIiLf4RycdhrMnQuffBJ2NbWjYBcR2cnw4TYhWFS7YxTsIiI76djRumQmTAi7ktpRsIuI7MLw\n4TBvHqxaFXYlNadgFxHZheHD7ePf/x5uHbWhYBcR2YX27aFr12h2xyjYRUR2Y/hwW3xj+fKwK6kZ\nBbuIyG580x0TtdExCnYRkd046CAoLIQnngi7kppRsIuI7MGwYbBwIaxYEXYliVOwi4jswZAh9vHp\np8OtoyYU7CIie9C+PXTuDE89FXYliVOwi4hUY+hQmzvm00/DriQxCnYRkWoMHWpzx0yaFHYliVGw\ni4hU47DDrEtm4sSwK0mMgl1EpBrO2Vl7eTmsXx92NdVTsIuIJGDoUFtZ6ZkILPypYBcRSUBBAbRt\nG43RMQp2EZEEOGdj2l98ESoqwq5mzxTsIiIJGjwYNm+GadPCrmTPFOwiIgnq0QOaNcv8YY8KdhGR\nBOXkwIAB8NxzsHVr2NXsnoJdRKQGBg+GDRvsTtRMpWAXEamBPn2gYcPM7o5RsIuI1ECTJtC7N0ye\nbNMMZCIFu4hIDQ0ebMvlLVkSdiW7pmAXEamhAQNsXPvkyWFXsmsKdhGRGsrLg2OPzdx+dgW7iEgt\nDBoEr78OK1eGXckPJRXszrlTnXNvO+d2OOcKgypKRCTTDRhgH597Ltw6diXZM/a3gCHAnABqERGJ\njI4d4aCD4Nlnw67kh5IKdu/9u977pUEVIyISFc7BySfDjBmwcWPY1Xyf+thFRGppwACbFGzmzLAr\n+T7nqxlh75x7CWi1i6fKvPeTq7aZBfzGe79oD/spBUoB8vLyCsaPH1+rgisqKsjNza3VazON2pJ5\n4tIOUFvSYetWx+DBPejVaw2//vX7Cb0mmbYUFxcv9t5Xfz3Te5/0A5gFFCa6fUFBga+t8vLyWr82\n06gtmScu7fBebUmXYcO8328/73fsSGz7ZNoCLPIJZKy6YkREkjBgAHzyCbz2WtiVfCvZ4Y6nOOdW\nAT8GnnPOZfj08yIiwerf3y6kTpkSdiXfSnZUzNPe+wO89w2893ne+5OCKkxEJApatIDu3WMU7CIi\nYt0xCxfCp5+GXYlRsIuIJKl/f/v4wgvh1vENBbuISJK6dIHWreH558OuxCjYRUSS5Bz06wfTp8O2\nbWFXo2AXEQlEv37w+ecwf37YlSjYRUQC0bs31K2bGd0xCnYRkQDssw/06AFTp4ZdiYJdRCQw/frZ\n4huffBJuHQp2EZGAZMqwRwW7iEhAOneG/fcPv59dwS4iEhDnoG9fG/a4dWt4dSjYRUQC1K8ffPll\nuMMeFewiIgE68UQb9jh9eng1KNhFRAK0zz5wzDEwLcRJzBXsIiIBO+kkWLQI1q0L5/0V7CIiAevT\nB7yHGTPCeX8Fu4hIwI4+Gpo2Da87RsEuIhKwunVt7php0+zMPd0U7CIiKdCnj00t8M476X9vBbuI\nSAqcVLUCdBjdMQp2EZEUaNMGDj1UwS4iEisnnQRz5sDXX6f3fRXsIiIp0qcPVFbC3LnpfV8Fu4hI\nihx/PNSrBy+9lN73VbCLiKRIkyZw7LEKdhGRWOndG157Lb3TCyjYRURSqHdvu0mpvDx976lgFxFJ\nocJC+NGP0tsdo2AXEUmhnBwoLlawi4jESu/e8NFH9kiHpILdOXebc+4959wbzrmnnXP7BFWYiEhc\n9O5tH9N11p7sGfuLwOHe+yOA94HfJV+SiEi8dOgABzTbyEuXTuGEXr2gXTsYNy5l75eTzIu9999d\n1W8+MCy5ckRE4sc9No7eX+zgme398YBbvhxKS+3JkpLA3y/IPvazgecD3J+ISDyUldF7+wusZ19e\n50j73qZNUFaWkrdzvppZ4J1zLwGtdvFUmfd+ctU2ZUAhMMTvZofOuVKgFCAvL69g/PjxtSq4oqKC\n3NzcWr0206gtmScu7QC1JZOc0KsX//Yt+QV/4Vqup5DFAHjnmD1zZsL7KS4uXuy9L6x2Q+99Ug/g\nZ8ArQONEX1NQUOBrq7y8vNavzTRqS+aJSzu8V1sySn6+93af0vcf+fk12g2wyCeQscmOiukL/BYY\n6L3flMy+RERia/RoaNz4+99r3Ni+nwLJ9rHfBewFvOice905d28ANYmIxEtJCYwdC/n5eOcgP9++\nTsGFU0h+VMx/BVWIiEislZRASQmzZ82iqKgopW+lO09FRGJGwS4iEjMKdhGRmFGwi4jEjIJdRCRm\nqr3zNCVv6txaYHktX94cSOMiUymltmSeuLQD1JZMlUxb8r33LarbKJRgT4ZzbpFP5JbaCFBbMk9c\n2gFqS6ZKR1vUFSMiEjMKdhGRmIlisI8Nu4AAqS2ZJy7tALUlU6W8LZHrYxcRkT2L4hm7iIjsQcYG\nu3Our3NuqXPuA+fcqF0875xzY6qef8M5d1QYdSYigbYUOee+qJoh83Xn3DVh1Fkd59wDzrk1zrm3\ndvN8JI5JAu2IxPEAcM61cc6VO+fecc697Zy7ZBfbROW4JNKWjD82zrmGzrl/OueWVLXj+l1sk9pj\nksik7el+AHWBD4GDgPrAEqDTTtv0x5bic0B3YEHYdSfRliJgSti1JtCW44GjgLd283xUjkl17YjE\n8aiqtTVwVNXne2GLykf1dyWRtmT8san6d86t+rwesADons5jkqln7N2AD7z3H3nvtwDjgUE7bTMI\neMSb+cA+zrnW6S40AYm0JRK893OA9XvYJBLHJIF2RIb3frX3/tWqz78C3gX232mzqByXRNqS8ar+\nnSuqvqxX9dj5YmZKj0mmBvv+wMrvfL2KHx7gRLbJBInWeWzVn2TPO+cOS09pgYvKMUlE5I6Hc64d\n0BU7Q/yuyB2XPbQFInBsnHN1nXOvA2uAF733aT0mSS20IYF5FWjrva9wzvUHJgHtQ64pm0XueDjn\ncoGngF95778Mu55kVNOWSBwb7/124Ejn3D7A0865w733u7ymkwqZesb+L6DNd74+oOp7Nd0mE1Rb\np/f+y2/+dPPeTwXqOeeap6/EwETlmOxR1I6Hc64eFoTjvPcTd7FJZI5LdW2J2rHx3n8OlAN9d3oq\npcckU4N9IdDeOXegc64+8FPgmZ22eQY4q+rqcnfgC+/96nQXmoBq2+Kca+Wcc1Wfd8OOy2dprzR5\nUTkmexTBZHi+AAAA20lEQVSl41FV5/3Au977O3azWSSOSyJticKxcc61qDpTxznXCPgJ8N5Om6X0\nmGRkV4z3fptz7kJgGjaq5AHv/dvOufOqnr8XmIpdWf4A2AT8PKx69yTBtgwDznfObQO+Bn7qqy6d\nZxLn3OPYqITmzrlVwLXYhaFIHZME2hGJ41GlB3Am8GZVny7AlUBbiNZxIbG2ROHYtAYeds7Vxf7j\nmeC9n5LO/NKdpyIiMZOpXTEiIlJLCnYRkZhRsIuIxIyCXUQkZhTsIiIxo2AXEYkZBbuISMwo2EVE\nYub/Af1UZhyp0QSqAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "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": 1, "metadata": { "collapsed": true }, "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": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9xvHvgjDHKggEVAjqRRBFxESkIpogReAiICDW\nG5VWaRzqUG2t1Dgr9arVWqqPSutYUaSKoIiCQoCCQgEVZ7wOTBULCA4Bw7juH79YFYGc5Oxz9tn7\nvJ/nOU+Gs88+v8VOXnbWXnst571HRETio07YBYiISLAU7CIiMaNgFxGJGQW7iEjMKNhFRGJGwS4i\nEjMKdhGRmFGwi4jEjIJdRCRmcsJ40+bNm/t27drV6rUbN26kSZMmwRYUErUl88SlHaC2ZKpk2rJ4\n8eJ13vsW1W0XSrC3a9eORYsW1eq1s2bNoqioKNiCQqK2ZJ64tAPUlkyVTFucc8sT2U5dMSIiMaNg\nFxGJGQW7iEjMKNhFRGJGwS4iEjMKdhGRmFGwi4jETCjj2EWy2eefw7JlsHatPT77DCorYcsWe9Sv\nDw0bQqNGsO++kJcHrVpBfr59T6Q6SQe7c64hMAdoULW/J7331ya7X5Go+/JLWLgQliyxx1tvwUcf\nWbDXVtu20KEDdOkCxxwD3bvDAQcEV7PEQxBn7JuBXt77CudcPWCuc+557/38APYtEhlffAEzZsDM\nmTB3Lrz5JuzYYc+1agWdO1sQH3ggtGtnZ+ItWthZeaNG0KAB5OTAtm12Br9pk53Nf/oprF4NH34I\nS5fCe+/BmDF2dg+2v7597dGrF+ROHgdlZZywYoX9TzB6NJSUhPbvIumXdLB77z1QUfVlvaqHT3a/\nIlHwf/8HEyfCs8/C/PmwfTvk5lqAX3MN/PjHcOSR0LJl4vusV88ee+1l4d+p0w+32bzZ/gqYP9/+\nM3nkEbjnHmhUfxsDt9fn9O1d6MtqGixfDqWl9iKFe9YIpI/dOVcXWAz8F3C3935BEPsVyUTLlsGj\nj8KECXZWDnDUUTBqFJx0koV6vXqpraFBA+jWzR4XX2xn73PnwlODHmNCRT+e4FSas5ZzuJ/zNt1L\nu7IyBXsWcXbCHdDOnNsHeBq4yHv/1k7PlQKlAHl5eQXjx4+v1XtUVFSQm5ubbKkZQW3JPLtrR2Vl\nHcrLW/LCC6144419AOjc+XOOP34tPXuuIy9vc7pL3aUTevVim6/LDE7kL/yCyQxiB3UYyDP0vacN\nHTt+FXaJtRKXny9Iri3FxcWLvfeF1W7ovQ/0AVwD/GZP2xQUFPjaKi8vr/VrM43aknl2bsc773h/\n0UXe77239+D9IYd4f9NN3n/8cSjlVS8/3wqteqxkf1/Gjb5pnQ0evO/Xz/tXXgm7yJqLy8+X98m1\nBVjkE8jhpMexO+daVJ2p45xrBPwEeC/Z/YqExXuYPt0uRnbqBPfdB/37w+zZduGyrMwufmak0aOh\nceP/fHkA/+Kmxjez7N5p3HyzjdL58Y9h+HAboSPxFMQNSq2BcufcG8BC4EXv/ZQA9iuSVtu2wYwZ\nLenSxfrKlyyBG2+EVavgscfg+OPBubCrrEZJCYwdC/n5eOds8PvYsfzoF6cxapRdH7j+enjuOTj0\nUPjtb2HjxrCLlqAlHeze+ze8912990d47w/33t8QRGEi6bJ1K9x/P3TsCDfd1Int2+GhhywEr7rK\nhiRGSkkJLFvG7JkzrRHfuWjapImN1nn/fTj9dLjtNjj8cJg2LbxyJXiaUkCy1rZt8PDDFugjR8I+\n+8ANN7zFm2/CiBE28iSu9t/f/vOaM8fa2bcvnHlmcjdPSeZQsEvW8d7Gnh9+OPzsZ7D33jBlivU/\n9+y5jjpZ9FvRsye8/jpcfTWMH293tP7jH2FXJcnKoh9hERvr3aMHDB1q/eVPPQWLF8N//3cE+s9T\npGFDuOEGmDfPxt8XFcG119rNVhJNCnbJCsuW2UiQnj3t87Fj7eaiIUOyN9B31q0bvPaadcnccAP0\n62dTGkj0KNgl1jZutG6Gjh2tu+Xaa20agF/8wuZlke/bay/re//rX21459FH2+ggiRYFu8TSN/3o\nhx4KN91kXS9Ll8J119nIENmzc86xvvYtW2zc+6RJYVckNaFgl9j58EO7oWjoUGja1AJq3Dho0ybs\nyqKlWze7/tCli3VZ/fnPYVckiVKwS2xs2QK//72Ndpk3D/70Jwum444Lu7Loysuz2SMHDrTJxn7z\nm2+nIpbMpV5GiYVXXrF+87ffhmHDLNT32y/squKhcWMbPXTJJXD77bBund3QVbdu2JXJ7ijYJdIq\nKuzu0DFjbCWhZ5+FAQPCrip+6ta1rpiWLe0CdGUl/O1vqZ+eWGpHwS6RNXOmXeRbtgx++Uu4+WYb\n1SGp4ZxNR9Cokc0xs3mz3dQU5zt0o0p97BI5FRUW5CeeaGeMc+bAXXcp1NPl8svt7H3SJDjtNJtr\nRzKLgl0iZfZsOOIIWwbu0kvtdviePcOuKvtceKGF++TJNi2D7lLNLOqKkUiorLS+9DvugIMOsoBX\noIfrwgvtr6ff/c4usI4dq7t4M4WCXTLe66/DGWfYiJfzzrOpZmOySlrkjRpl4T56NDRrBrfcEnZF\nAgp2yWA7dtjwurIyaN4cpk61+Usks9x4I6xfD7feCm3b2vUPCZeCXTLSypVw1lkwa5bd9Th2LOy7\nb9hVya44Z8NNV62ym5jatLEbmiQ8ungqGWfiRLuNfeFCuxHmyScV6pkuJwcefxyOOgp++lP45z/D\nrii7KdglY2zaBKWlNsfLwQfbFLJnn60LclHRpInNoJmXB4MHwyefhF1R9lKwS0Z4800oLLTpYq+4\nwuZ6ad8+7KqkpvLy4Jln4MsvrQutsjLsirKTgl1C5T3ce6/NJLhhA0yfDv/7v1C/ftiVSW117gyP\nPAILFtgoJu/Drij7KNglNJ9/bqsanX8+nHCCLejQu3fYVUkQhgyxOWUefljT/YZBwS6hWLjQLrRN\nmmTD5KZOtQmmJD6uuQYGDYJf/xrmzw+7muyiYJe08h7uvNMWlN6+3eZ5ufxyqKOfxNipUwcefNBm\n3TztNBvrLumhXydJmw0b7E/0Sy+1FY5ee82WXZP4atoU/v53+PRTGDFCi3Ski4Jd0uKbrpcpU+CP\nf4Snn7Zb0CX+CgvtDuIpU+APfwi7muygYJeU8t6m1O3Rw87W5s6FX/1KY9OzzS9/afcnlJXZcoWS\nWgp2SZkvv7S+1YsugpNOsq6XY44JuyoJg3M2LUReHpSU2M1okjoKdkmJN96wP8EnTrRRL5Mnq+sl\n2zVrZsMfly61RbEldRTsErgHH7Qz84oKKC/XqBf51okn2vDHe+6xPndJjaR/3ZxzbZxz5c65d5xz\nbzvnLgmiMImer7+2NUjPPhuOPda6XrQYhuxs9GhbBWvkSPjss7CriacgzqO2Ab/23ncCugO/dM51\nCmC/EgXjxkG7duxXXEr3vd/lgQdspaPp060/VWRnDRrYlAOffWYX0iV4SQe793619/7Vqs+/At4F\n9k92vxIB48ZBaSlPLS+gkEWs2tqSqQ1O4caO46hbN+ziJJN16QJXXgmPPqoumVQItOfTOdcO6Aos\nCHK/kpm2Xnktl226kWE8xaG8y2t0pd/mSTamTaQaZWU2Ydi559q8QRIc5wOaes05lwvMBkZ77yfu\n4vlSoBQgLy+vYPz48bV6n4qKCnJjsuBllNuydm0D7hy+g5fpwUWM4Q/8hvpsBcA7x+yZM0OusHai\nfEx2FoW2LF2aywUXFNCnz6dcccXS3W4XhbYkKpm2FBcXL/beF1a7ofc+6QdQD5gGXJbI9gUFBb62\nysvLa/3aTBPVtkyf7n3z5t7nuq/8E5zqvd2H9O0jPz/sEmstqsdkV6LSllGj7MdmT+VGpS2JSKYt\nwCKfQMYGMSrGAfcD73rv70h2f5K5tm+H666zm41atYJFt8xkeOPnvr9R48Y27EEkQVdfDQceaNM3\nb94cdjXxEEQfew/gTKCXc+71qkf/APYrGWTtWpu46/rr4cwzbRGFDpcPtNsJ8/PxzkF+vn1dUhJ2\nuRIhjRvD3XfDe+9pLpmg5CS7A+/9XEAzf8TY3Lm2QPG6dfCXv9hY9f/M9VJSAiUlzJ41i6KiojDL\nlAjr1w9OPRVuusl+1g4+OOyKok33A8pueW9nUEVF0KiRLZYwcqQm8JLUuPNOqFcPLrhAy+klS8Eu\nu7R+va1+c/nltuL8okVw5JFhVyVxtt9+dsY+fbrNLSS1p2CXH1iwwOZOf+EF+NOfbKGEvfcOuyrJ\nBuefD4cdBpddBpWVYVcTXQp2+Q/vbRGM446z7pZ58+Dii9X1IulTr56dTHz8MdyhMXa1pmAXwObt\nGDjQzpQGDIBXX4Wjjw67KslGJ54Ip5xio2ZXrQq7mmhSsAvz5kHXrta3OWaMzaHetGnYVUk2u/12\nu2/iiivCriSaFOxZbPt2Oys64QT7E/jll221I3W9SNgOPNAW43jsMbvmIzWjYM9Sq1fbHaRXXQXD\nh9vc6QUFYVcl8q0rroCWLW1kloY/1oyCPQtNmWILHbzyCjzwgM2++6MfhV2VyPfttZdNYfGPf8C8\nefuGXU6kKNizSGWljXI5+WQ44ABbLf7nP1fXi2SukSOhQwcYO/Zgtm4Nu5roULBniTfesFEuf/6z\nrVozfz507Bh2VSJ7Vq8e3HILrFzZmL/+NexqokPBHnM7dtjY9KOPtom8pk61rxs0CLsykcQMHAhH\nHPE5111nC6RL9RTsMbZypV0gvewy6NsX3nzTJlsSiRLn4NxzP2TNGrt5SaqnYI8h720tyc6d7QLp\nfffBpEnQokXYlYnUTqdOX3HyyXDbbbBhQ9jVZD4Fe8ysWWPTn555Jhx+OCxZAqWlukAq0XfjjfDF\nF5qzPREK9hh58kmbQOnZZ+Hmm2H2bM1rLfHRpQucdpp1x6xZE3Y1mU3BHgNr1tjiBKeeaosYvfoq\njBoFdeuGXZlIsK6/Hr7+2k5cZPcU7BHmvd1c1KkTPP20zWX9yit21i4SRx06wIgRcM89miBsTxTs\nEbVsmc3CeMYZ0L69TQlQVmbjfkXi7JprbJ6jW28Nu5LMpWCPmG3b7OLRYYdZH/qdd9qapJ06hV2Z\nSHq0awdnnWXrpq9eHXY1mUnBHiEvvwyFhTYp0oknwjvvwCWXqC9dss+VV9pJzm23hV1JZlKwR8Ca\nNTanS48etiDGU0/ZmpBt24ZdmUg4Dj4YSkrg3nvh3/8Ou5rMo2DPYFu22O3/hxxiF0lHjYL33oMh\nQzQuXeTKK2HzZluUQ75PwZ6BvLepdTt3tukAune3SbxuvhmaNAm7OpHM0KGDDfO9+26bB0m+pWDP\nMAsXQq9eNrWuc/Dcc/DCC5qJUWRXyspg0yabtVS+pWDPEO+/b2cf3brBW2/Z2qNvvgn9+4ddmUjm\n6tQJBg2Cu+7SzI/fpWAP2ccf24XRQw+1qQCuugo+/NDWHtWYdJHqjRplE4NpvvZvKdhD8v77cOut\nHTjkEHj8cRu2+NFHNtGRlqkTSVz37rYg++2324ADUbCn3aJFNpFRx44wY0ZLzj/fztDvuAPy8sKu\nTiSarrjCphh4/PGwK8kMCvY02L7d5kM//nhbyej55+0H8fHH5zNmDOy/f9gVikRb3762QPstt9iq\nYdkukGB3zj3gnFvjnHsriP3t0rhx0K4dJ/TqZfcUjxuXsrcKyurVNjHXQQfBKafAihX25+LKlTZ0\nsVkzrc4rEgTn7GTp3XdtJFm2C+qM/SGgb0D7+qFx42y1iOXLcd7D8uX2dQaGe2WlzYs+aJDdGXr1\n1XaD0ZNPwgcf2Lj0vfcOu0qR+Bk+3H7n7rgj7ErCF0iwe+/nAOuD2NcuVQ1WXUA3nmUAW6hng1fL\nylL2ljWxebOdJZx9NrRubfOiL1wIl15qF0lffBGGDoWcnLArFYmvnBwbTTZrls12ms2c9z6YHTnX\nDpjivT98N8+XAqUAeXl5BePHj0943yf06oXznhE8xCOMoCnrGcpTDGcCftpV1K8fTBtqYu3aBvzz\nn01ZuLAZixY1Y+PGHJo02caxx66jT59/07Xrhmon56qoqCA3Nzc9BadYXNoSl3ZAdraloiKH4cO7\nc9xx67jyyvfSUFnNJXNciouLF3vvC6vd0HsfyANoB7yVyLYFBQW+RvLzvQe/hRz/HP38GTzim/CV\nB++bNPH+5JO9HzPG+4ULvd+8uWa7TsTWrd6/8473Dzzg/TnneN+hg/d247/3++3n/ciR3k+dWvP3\nLi8vD77YkMSlLXFph/fZ25aLL/Y+J8f7VatSV08ykjkuwCKfQMZGo3Ng9GgoLaXepk3053n68zyb\nGu3LS+c9ybTNRTz/vN3cA9CwIRx5pN3w06GD9W+3bg2tWkHLltCo0Q8n0Nqxw25wWLvWHitW2Jjy\nDz+Et9+2O0ErK23bZs3g2GNh5Ei7En/YYZqQSySTXHKJTTFw993w+9+HXU04ohHsJSX2sawMv2IF\nrm1bGo8ezcCSIgZi584rVsCCBfZYvBimToUHH9z17ho2hAYNbD7nLVtg624Gp7RubbcsX3CBLaRb\nWGjjz+tokKhIxjroIBg8GO67zy7DZePEeYEEu3PucaAIaO6cWwVc672/P4h9/0dJCZSUMHvWLIqK\ninZ6f1vEOT/frox/44svbCTKp5/anM1r1thCuJWV9sjJgfr1LeSbNoUWLaB5c7uynp9vZ/ciEj2X\nXWbrAP/tb3DeeWFXk36BBLv3/vQg9hO0vfeGgoKwqxCRdOvRA7p2tcnBzj03+7pL1akgIrHjHFx4\noV0jmz077GrST8EuIrF0+uk22OGuu8KuJP0U7CISS40a2ei1SZNsGo9somAXkdg6/3wbznzffWFX\nkl4KdhGJrXbtbJnJsWNt6o9soWAXkVi78EK78XDChLArSR8Fu4jEWu/edgf6vfeGXUn6KNhFJNac\ns1m+X37ZpgfJBgp2EYm9ESPsLvNsuYiqYBeR2GveHIYNsykGNm0Ku5rUU7CLSFY491ybP+qJJ8Ku\nJPUU7CKSFXr2tOm8s6E7RsEuIlnhm4uoCxbAkiVhV5NaCnYRyRpnnWXrMYwdG3YlqaVgF5Gs0awZ\nDBkCjz1mazPElYJdRLLKOefA55/b5GBxpWAXkaxSVAQHHgj3B7vGW0ZRsItIVqlTB37+c5gxAz7+\nOOxqUkPBLiJZZ8QIGyXz0ENhV5IaCnYRyTpt28JPfgIPPgjbt4ddTfAU7CKSlc45x1ZWmjEj7EqC\np2AXkaw0aJANf3zggbArCZ6CXUSyUoMGtuD15Mk2h0ycKNhFJGuddRZUVsKTT4ZdSbAU7CKStY4+\nGjp0gEceCbuSYCnYRSRrOWdn7XPmxGtMu4JdRLLaGWfYx0cfDbeOICnYRSSrtW0LxcXWHeN92NUE\nQ8EuIlnvrLPggw9g/vywKwmGgl1Est7QodCoUXwuogYS7M65vs65pc65D5xzo4LYp4hIuuy1Fwwe\nDBMmwJYtYVeTvKSD3TlXF7gb6Ad0Ak53znVKdr8iIulUUgLr18OLL4ZdSfKCOGPvBnzgvf/Ie78F\nGA8MCmC/IiJp06cP7Luvra4Udc4neRnYOTcM6Ou9H1n19ZnAMd77C3farhQoBcjLyysYP358rd6v\noqKC3NzcpGrOFGpL5olLO0BtqY0//rE906e3YuLEeTRqtCMl75FMW4qLixd77wur2y6nVnuvBe/9\nWGAsQGFhoS8qKqrVfmbNmkVtX5tp1JbME5d2gNpSG3XrwjPPwIYNx9OvX2reIx1tCaIr5l9Am+98\nfUDV90REIqVHD2jTJvrdMUEE+0KgvXPuQOdcfeCnwDMB7FdEJK3q1LEZH6dNg3Xrwq6m9pIOdu/9\nNuBCYBrwLjDBe/92svsVEQnD//wPbNsW7RkfAxnH7r2f6r0/xHt/sPd+dBD7FBEJwxFHQKdO0e6O\n0Z2nIiLf4RycdhrMnQuffBJ2NbWjYBcR2cnw4TYhWFS7YxTsIiI76djRumQmTAi7ktpRsIuI7MLw\n4TBvHqxaFXYlNadgFxHZheHD7ePf/x5uHbWhYBcR2YX27aFr12h2xyjYRUR2Y/hwW3xj+fKwK6kZ\nBbuIyG580x0TtdExCnYRkd046CAoLIQnngi7kppRsIuI7MGwYbBwIaxYEXYliVOwi4jswZAh9vHp\np8OtoyYU7CIie9C+PXTuDE89FXYliVOwi4hUY+hQmzvm00/DriQxCnYRkWoMHWpzx0yaFHYliVGw\ni4hU47DDrEtm4sSwK0mMgl1EpBrO2Vl7eTmsXx92NdVTsIuIJGDoUFtZ6ZkILPypYBcRSUBBAbRt\nG43RMQp2EZEEOGdj2l98ESoqwq5mzxTsIiIJGjwYNm+GadPCrmTPFOwiIgnq0QOaNcv8YY8KdhGR\nBOXkwIAB8NxzsHVr2NXsnoJdRKQGBg+GDRvsTtRMpWAXEamBPn2gYcPM7o5RsIuI1ECTJtC7N0ye\nbNMMZCIFu4hIDQ0ebMvlLVkSdiW7pmAXEamhAQNsXPvkyWFXsmsKdhGRGsrLg2OPzdx+dgW7iEgt\nDBoEr78OK1eGXckPJRXszrlTnXNvO+d2OOcKgypKRCTTDRhgH597Ltw6diXZM/a3gCHAnABqERGJ\njI4d4aCD4Nlnw67kh5IKdu/9u977pUEVIyISFc7BySfDjBmwcWPY1Xyf+thFRGppwACbFGzmzLAr\n+T7nqxlh75x7CWi1i6fKvPeTq7aZBfzGe79oD/spBUoB8vLyCsaPH1+rgisqKsjNza3VazON2pJ5\n4tIOUFvSYetWx+DBPejVaw2//vX7Cb0mmbYUFxcv9t5Xfz3Te5/0A5gFFCa6fUFBga+t8vLyWr82\n06gtmScu7fBebUmXYcO8328/73fsSGz7ZNoCLPIJZKy6YkREkjBgAHzyCbz2WtiVfCvZ4Y6nOOdW\nAT8GnnPOZfj08yIiwerf3y6kTpkSdiXfSnZUzNPe+wO89w2893ne+5OCKkxEJApatIDu3WMU7CIi\nYt0xCxfCp5+GXYlRsIuIJKl/f/v4wgvh1vENBbuISJK6dIHWreH558OuxCjYRUSS5Bz06wfTp8O2\nbWFXo2AXEQlEv37w+ecwf37YlSjYRUQC0bs31K2bGd0xCnYRkQDssw/06AFTp4ZdiYJdRCQw/frZ\n4huffBJuHQp2EZGAZMqwRwW7iEhAOneG/fcPv59dwS4iEhDnoG9fG/a4dWt4dSjYRUQC1K8ffPll\nuMMeFewiIgE68UQb9jh9eng1KNhFRAK0zz5wzDEwLcRJzBXsIiIBO+kkWLQI1q0L5/0V7CIiAevT\nB7yHGTPCeX8Fu4hIwI4+Gpo2Da87RsEuIhKwunVt7php0+zMPd0U7CIiKdCnj00t8M476X9vBbuI\nSAqcVLUCdBjdMQp2EZEUaNMGDj1UwS4iEisnnQRz5sDXX6f3fRXsIiIp0qcPVFbC3LnpfV8Fu4hI\nihx/PNSrBy+9lN73VbCLiKRIkyZw7LEKdhGRWOndG157Lb3TCyjYRURSqHdvu0mpvDx976lgFxFJ\nocJC+NGP0tsdo2AXEUmhnBwoLlawi4jESu/e8NFH9kiHpILdOXebc+4959wbzrmnnXP7BFWYiEhc\n9O5tH9N11p7sGfuLwOHe+yOA94HfJV+SiEi8dOgABzTbyEuXTuGEXr2gXTsYNy5l75eTzIu9999d\n1W8+MCy5ckRE4sc9No7eX+zgme398YBbvhxKS+3JkpLA3y/IPvazgecD3J+ISDyUldF7+wusZ19e\n50j73qZNUFaWkrdzvppZ4J1zLwGtdvFUmfd+ctU2ZUAhMMTvZofOuVKgFCAvL69g/PjxtSq4oqKC\n3NzcWr0206gtmScu7QC1JZOc0KsX//Yt+QV/4Vqup5DFAHjnmD1zZsL7KS4uXuy9L6x2Q+99Ug/g\nZ8ArQONEX1NQUOBrq7y8vNavzTRqS+aJSzu8V1sySn6+93af0vcf+fk12g2wyCeQscmOiukL/BYY\n6L3flMy+RERia/RoaNz4+99r3Ni+nwLJ9rHfBewFvOice905d28ANYmIxEtJCYwdC/n5eOcgP9++\nTsGFU0h+VMx/BVWIiEislZRASQmzZ82iqKgopW+lO09FRGJGwS4iEjMKdhGRmFGwi4jEjIJdRCRm\nqr3zNCVv6txaYHktX94cSOMiUymltmSeuLQD1JZMlUxb8r33LarbKJRgT4ZzbpFP5JbaCFBbMk9c\n2gFqS6ZKR1vUFSMiEjMKdhGRmIlisI8Nu4AAqS2ZJy7tALUlU6W8LZHrYxcRkT2L4hm7iIjsQcYG\nu3Our3NuqXPuA+fcqF0875xzY6qef8M5d1QYdSYigbYUOee+qJoh83Xn3DVh1Fkd59wDzrk1zrm3\ndvN8JI5JAu2IxPEAcM61cc6VO+fecc697Zy7ZBfbROW4JNKWjD82zrmGzrl/OueWVLXj+l1sk9pj\nksik7el+AHWBD4GDgPrAEqDTTtv0x5bic0B3YEHYdSfRliJgSti1JtCW44GjgLd283xUjkl17YjE\n8aiqtTVwVNXne2GLykf1dyWRtmT8san6d86t+rwesADons5jkqln7N2AD7z3H3nvtwDjgUE7bTMI\neMSb+cA+zrnW6S40AYm0JRK893OA9XvYJBLHJIF2RIb3frX3/tWqz78C3gX232mzqByXRNqS8ar+\nnSuqvqxX9dj5YmZKj0mmBvv+wMrvfL2KHx7gRLbJBInWeWzVn2TPO+cOS09pgYvKMUlE5I6Hc64d\n0BU7Q/yuyB2XPbQFInBsnHN1nXOvA2uAF733aT0mSS20IYF5FWjrva9wzvUHJgHtQ64pm0XueDjn\ncoGngF95778Mu55kVNOWSBwb7/124Ejn3D7A0865w733u7ymkwqZesb+L6DNd74+oOp7Nd0mE1Rb\np/f+y2/+dPPeTwXqOeeap6/EwETlmOxR1I6Hc64eFoTjvPcTd7FJZI5LdW2J2rHx3n8OlAN9d3oq\npcckU4N9IdDeOXegc64+8FPgmZ22eQY4q+rqcnfgC+/96nQXmoBq2+Kca+Wcc1Wfd8OOy2dprzR5\nUTkmexTBZHi+AAAA20lEQVSl41FV5/3Au977O3azWSSOSyJticKxcc61qDpTxznXCPgJ8N5Om6X0\nmGRkV4z3fptz7kJgGjaq5AHv/dvOufOqnr8XmIpdWf4A2AT8PKx69yTBtgwDznfObQO+Bn7qqy6d\nZxLn3OPYqITmzrlVwLXYhaFIHZME2hGJ41GlB3Am8GZVny7AlUBbiNZxIbG2ROHYtAYeds7Vxf7j\nmeC9n5LO/NKdpyIiMZOpXTEiIlJLCnYRkZhRsIuIxIyCXUQkZhTsIiIxo2AXEYkZBbuISMwo2EVE\nYub/Af1UZhyp0QSqAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "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": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADcBJREFUeJzt3WGoXGedx/Hfb5NKbaaLL1pnJam5vhAhlG3LXEqXLrv3\nhFWiFmUXFhqir4R5s0plFVEKgkLZN4vsGxf2osUFUwfZNaxEsUQ61yJI6701dZOmhdJNNUHIFnG3\n04AS/e+LmUK7pMm9z3mSc+af7weG3Jl7znP+/3vSX06fee4ZR4QAAHn8UdcFAADqItgBIBmCHQCS\nIdgBIBmCHQCSIdgBIBmCHQCSIdgBIBmCHQCS2d3FQW+77bZYWVkp2ve1117Tnj176hbUEXrpnyx9\nSPTSV2162draeiUibr/adp0E+8rKijY3N4v23djY0NraWt2COkIv/ZOlD4le+qpNL7Zf3s52TMUA\nQDIEOwAkQ7ADQDIEOwAkQ7ADQDJVVsXYPivpVUm/l3QpIlZrjAsA2Lmayx2biHil4ngAgAJMxQBA\nMq7xmae2/0vS/2g+FfMvEbF+mW3GksaSNBwOR5PJpOhYs9lMg8GgRbX9sb6+rvF43HUZVWQ5L1n6\nkOilr9r00jTN1ramuiOi9UPS3sWf75T0rKS/uNL2o9EoSk2n0+J9+2b+488hy3nJ0kcEvfRVm14k\nbcY2MrnKVExEnF/8eUHSMUn31hgXALBzrYPd9h7bt77+taQPSDrVdlwAQJkaq2KGko7Zfn28xyLi\nBxXGBQAUaB3sEfGSpLsq1AIAqIDljgCQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ\n7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQDMEOAMkQ7ACQ\nDMEOAMkQ7ACQDMEOAMlUC3bbu2z/zPbxWmMCAHau5hX7Q5LOVBwPAFCgSrDb3ifpw5K+VmM8AEA5\nR0T7Qex/k/QPkm6V9NmIeOAy24wljSVpOByOJpNJ0bFms5kGg0GLavujaRpNp9Ouy6giy3nJ0odE\nL33VppemabYiYvWqG0ZEq4ekByT98+LrNUnHr7bPaDSKUtPptHjfvpn/+HPIcl6y9BFBL33VphdJ\nm7GNXK4xFXO/pI/YPitpIumg7W9WGBcAUKB1sEfEFyJiX0SsSHpQ0hMR8bHWlQEAirCOHcjm6FFp\nZUV/efCgtLIyf44byu6ag0XEhqSNmmMC2IGjR6XxWLp4UZakl1+eP5ekI0e6rAzXEVfsQCYPPyxd\nvPjm1y5enL+OGwbBDmTyi1/s7HWkRLADmbz73Tt7HSkR7EAmjzwi3XLLm1+75Zb567hhEOxAJkeO\nSOvr0v79Clvav3/+nDdObygEO5DNkSPS2bP60RNPSGfPEuo3IIIdAJIh2AEgGYIdAJIh2AEgGYId\nAJIh2AGJG2chlao3AQOWEjfOQjJcsQPcOAvJEOwAN85CMgQ7wI2zkAzBDnDjLCRDsAPcOAvJEOyA\nxI2z+oplqEVY7gign1iGWowrdgD9xDLUYgQ7gH5iGWoxgh1AP7EMtVjrYLd9s+2nbT9r+7TtL9Uo\nDMANjmWoxWpcsf9W0sGIuEvS3ZIO2b6vwrgAbmQsQy3WOthjbrZ4etPiEW3HBQCWoZapMsdue5ft\nk5IuSDoREU/VGBcAsHOOqHdxbfsdko5J+lREnPp/3xtLGkvScDgcTSaTomPMZjMNBoO2pfZC0zSa\nTqddl1FFlvOSpQ+JXvqqTS9N02xFxOrVtqsa7JJk+4uSLkbEP77VNqurq7G5uVk0/sbGhtbW1gqr\n6xfbqv3z70qW85KlD4le+qpNL7a3Few1VsXcvrhSl+23S3q/pOfbjgsAKFPjlgLvkvSvtndp/g/F\ntyPieIVxAQAFaqyK+XlE3BMRfxoRd0bEl2sUhiXBTZqA3uEmYCjHTZqAXuKWAijHTZqAXiLYUY6b\nNAG9RLCjHDdpAnqJYEc5btIE9BLBjnLcpAnoJYId7XCTJqB3CHYASIZgB4BkCHYASIZgB4BkCHYA\nSIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASKZ1\nsNu+w/bU9nO2T9t+qEZhAIAyuyuMcUnSZyLiGdu3StqyfSIinqswNgBgh1pfsUfEryLimcXXr0o6\nI2lv23EBAGWqzrHbXpF0j6Snao4LANg+R0SdgeyBpB9JeiQivnOZ748ljSVpOByOJpNJ0XFms5kG\ng0GbUnujaRpNp9Ouy6giy3nJ0odEL33VppemabYiYvWqG0ZE64ekmyQ9Lunvt7P9aDSKUtPptHjf\nvpn/+HPIcl6y9BFBL33VphdJm7GNjK2xKsaSvi7pTER8pe14AIB2asyx3y/p45IO2j65eHyowrgA\ngAKtlztGxI8luUItAIAK+M1TAEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiG\nYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeAZAh2AEiGYAeA\nZAh2AEiGYAeAZKoEu+1HbV+wfarGeACAcrWu2L8h6VClsQAALVQJ9oh4UtKva4wFAGiHOXYASMYR\nUWcge0XS8Yi48y2+P5Y0lqThcDiaTCZFx5nNZhoMBoVV9kvTNJpOp12XUUWW85KlD4le+qpNL03T\nbEXE6lU3jIgqD0krkk5tZ9vRaBSlptNp8b59M//x55DlvGTpI4Je+qpNL5I2YxsZy1QMACRTa7nj\ntyT9RNL7bJ+z/Yka4wIAdm53jUEi4nCNcQAA7TEVAwDJEOwAkAzBDgDJEOwAkAzBDgDJEOwAkAzB\nDgDJEOwAkAzBDgDJEOwAkAzBDgDJEOwAkAzBDgDJEOwAkAzBDgDJEOwAkAzBDgDJEOwAkAzBDgDJ\nEOwAkAzBDgDJEOwAkAzBDgDJEOwAkEyVYLd9yPYLtl+0/fkaYwIAyrQOdtu7JH1V0gclHZB02PaB\ntuMCAMrUuGK/V9KLEfFSRPxO0kTSRyuMCwAoUCPY90r65Ruen1u8BgDowO7rdSDbY0ljSRoOh9rY\n2CgaZ319XU3TVKysW7a7LgHAdXT48OFrfowawX5e0h1veL5v8dqbRMS6pHVJWl1djbW1teIDPvbY\nY8X79oltRUTXZVSxsbGhNue0L7L0IdFLX12PXmpMxfxU0nttv8f22yQ9KOm7FcYFABRofcUeEZds\nf1LS45J2SXo0Ik63rgwAUKTKHHtEfF/S92uMBQBoh988BYBkCHYASIZgB4BkCHYASIZgB4BkCHYA\nSIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASIZg\nB4BkCHYASIZgB4BkCHYASIZgB4BkCHYASKZVsNv+W9unbf/B9mqtogAA5dpesZ+S9DeSnqxQCwCg\ngt1tdo6IM5Jku041AIDWmGMHgGQcEVfewP6hpD+5zLcejoj/WGyzIemzEbF5hXHGksaSNBwOR5PJ\npKjg2WymwWBQtG/frK+vazwed11GFVnOS5Y+JHrpqza9NE2zFRFXfz8zIlo/JG1IWt3u9qPRKEpN\np9PiffuGXvonSx8R9NJXbXqRtBnbyFimYgAgmbbLHf/a9jlJfybpe7Yfr1MWAKBU21UxxyQdq1QL\nAKACpmIAIBmCHQCSIdgBIBmCHQCSIdgBIJmr/ubpNTmo/d+SXi7c/TZJr1Qsp0v00j9Z+pDopa/a\n9LI/Im6/2kadBHsbtjdjO79SuwTopX+y9CHRS19dj16YigGAZAh2AEhmGYN9vesCKqKX/snSh0Qv\nfXXNe1m6OXYAwJUt4xU7AOAKljLYl/1DtG0fsv2C7Rdtf77retqw/ajtC7ZPdV1LG7bvsD21/dzi\n79ZDXddUyvbNtp+2/eyily91XVMbtnfZ/pnt413X0obts7b/0/ZJ22/5oUQ1LGWwa4k/RNv2Lklf\nlfRBSQckHbZ9oNuqWvmGpENdF1HBJUmfiYgDku6T9HdLfF5+K+lgRNwl6W5Jh2zf13FNbTwk6UzX\nRVTSRMTdLHe8jIg4ExEvdF1HoXslvRgRL0XE7yRNJH2045qKRcSTkn7ddR1tRcSvIuKZxdevah4k\ne7utqsziw3Zmi6c3LR5L+Waa7X2SPizpa13XskyWMtiX3F5Jv3zD83Na0gDJyvaKpHskPdVtJeUW\n0xcnJV2QdCIilrWXf5L0OUl/6LqQCkLSD21vLT4D+ppp9UEb19J2PkQbqM32QNK/S/p0RPxv1/WU\niojfS7rb9jskHbN9Z0Qs1fsgth+QdCEitmyvdV1PBX8eEedtv1PSCdvPL/6Pt7reBntE/FXXNVwj\n5yXd8Ybn+xavoWO2b9I81I9GxHe6rqeGiPiN7anm74MsVbBLul/SR2x/SNLNkv7Y9jcj4mMd11Uk\nIs4v/rxg+5jm07LXJNiZirn+firpvbbfY/ttkh6U9N2Oa7rh2bakr0s6ExFf6bqeNmzfvrhSl+23\nS3q/pOe7rWrnIuILEbEvIlY0/+/kiWUNddt7bN/6+teSPqBr+A/tUgb7Mn+IdkRckvRJSY9r/gbd\ntyPidLdVlbP9LUk/kfQ+2+dsf6LrmgrdL+njkg4ulqOdXFwpLqN3SZra/rnmFxInImKplwomMJT0\nY9vPSnpa0vci4gfX6mD85ikAJLOUV+wAgLdGsANAMgQ7ACRDsANAMgQ7ACRDsANAMgQ7ACRDsANA\nMv8H1yLPSwdFRJcAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "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", " \n", "\n", " \n", "\n", " \n", "\n", "Simpson則の導出（数式変形）：\n", "\n", " \n", "\n", " \n", "\n", " \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": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9//H3F4iKRESKBoTQKEQUUIEgIiASFRXqdG8d\noL0OqFAcqlZR63BtxXJrW7UVZ1Gu1f40aBUHLo4YVKwT80wZtA6lzgUDVgW+vz/WicY0cE5OTrLP\n2efzep79cIZ1Tr7L/fjJztp7r2XujoiIxEuzqAsQEZHMU7iLiMSQwl1EJIYU7iIiMaRwFxGJIYW7\niEgMKdxFRGJI4S4iEkMKdxGRGGoR1Q9u166dl5SUpPXZDRs20KpVq8wWFBH1JTvFpS9x6QeoL9Xm\nzJnzsbvvmqxdZOFeUlLC7Nmz0/rszJkzGTJkSGYLioj6kp3i0pe49APUl2pm9rdU2mlYRkQkhhTu\nIiIxpHAXEYkhhbuISAwp3EVEYijlcDez5mY2z8ym1fGemdlEM1tlZgvNrE9myxQRkfqoz5H7BcCy\nrbw3DChNbGOA2xtYl4iINEBK17mbWSfgB8AE4KI6mhwH3Odhzb7XzKyNmXVw97WZKzVhyRJKJk+G\nefNgl12gXTvYd1/o3BnMMv7jRERykaWyhqqZ/Rn4NbATMM7dj671/jTgOneflXg+A7jM3WfXajeG\ncGRPUVFRWUVFRb0L3rWyku7XXovVqvvr1q35fK+9+LRfPz46+GC+bN++3t8dhaqqKgoLC6MuIyPU\nl+wTl36A+lKtvLx8jrv3TdrQ3be5AUcDtyUeDwGm1dFmGjCoxvMZQN9tfW9ZWZmnq3LGDPdPP3Vf\nvdr9L39xv/1299Gj3Xv2dIew9e3rfued7hs3pv1zmkJlZWXUJWSM+pJ94tIPd/WlGjDbk+S2u6c0\nLDMQONbMhgM7AK3N7E/u/l812rwPFNd43inxWuNo1iwMyeyyC+y5Jxx00LfvrVoFjz4KDzwAP/kJ\nXHUVnHMOnH8+tG3baCWJiGSTpCdU3f1yd+/k7iXACOCFWsEO8ARwauKqmf7AOm+M8fZUdO0Kl14a\nxuRfeAH69YNrroHSUrj5Zvj660jKEhFpSmlf525mY81sbOLpdGANsAqYBJyTgdoaxgzKy2HaNFiw\nAHr3Dkfv++0XQl9EJMbqFe7uPtMTJ1Pd/Q53vyPx2N39XHfv4u77eq0TqZHbbz947jl4/HHYtAkO\nOywE/caNUVcmItIo8ucOVTM49thwFH/++WGIpndvePPNqCsTEcm4/An3ajvuCDfdFIZm/vUvGDQI\nJk2KuioRkYzKv3CvVl4eTroOGQJjxsDo0fDll1FXJSKSEfkb7hAujZw+Ha64Au6+Gw49FD75JOqq\nREQaLL/DHaB5c5gwAR56CObMgYED4e23o65KRKRBFO7VTjwRnn0WPvgg3BQ1b17UFYmIpE3hXtPg\nwTBrFhQUhDH511+PuiIRkbQo3Gvr0SMEfLt2MHQovPJK1BWJiNSbwr0unTvDiy9C+/Zw5JHw0ktR\nVyQiUi8K963p2DEEfHExDB+uIRoRySkK923p0CHc7NS+PQwbBgsXRl2RiEhKFO7JdOgAzz8f7mw9\n4ghYuTLqikREklK4p6KkJEw8tnlzOMm6NprZjEVEUqVwT9U++8DTT8PHH8PRR8Pnn0ddkYjIVinc\n66OsLNzJumABnHSSFv4QkaylcK+v4cPhjjvCUfzZZ4cVW0VEskwqa6hKbWedBe+8A9deC926wSWX\nRF2RiMh3KNzTdc01sGIFXHYZ7L03HHNM1BWJiHwj6bCMme1gZm+Y2QIzW2Jm19TRZoiZrTOz+Ynt\n6sYpN4uYwf/+bxiH/9GPYNGiqCsSEflGKmPuXwKHuvv+QC/gKDPrX0e7l929V2Ibn9Eqs9WOO4Z1\nWVu3DkfuH38cdUUiIkAK4Z5Y/Loq8bQgseksYrXddw8B/49/wMknhwW4RUQiltLVMmbW3MzmAx8C\nz7l7XROtDDCzhWb2lJn1yGiV2a5vX7jrrjBVwWWXRV2NiAjm9biUz8zaAFOBn7r74hqvtwa2uHuV\nmQ0HbnL30jo+PwYYA1BUVFRWUVGRVtFVVVUUFham9dnG1HXiRDpNncrSK6/kw8MPT+kz2dqXdKgv\n2Scu/QD1pVp5efkcd++btKG712sDrgbGJWnzNtBuW23Kyso8XZWVlWl/tlF99ZX7wQe7t2zpPn9+\nSh/J2r6kQX3JPnHph7v6Ug2Y7SlkdSpXy+yaOGLHzFoCQ4Hltdq0NzNLPO5HGO7Jv5WmCwrg4Yeh\nTRs44QRYty7qikQkT6Uy5t4BqDSzhcCbhDH3aWY21szGJtqcACw2swXARGBE4jdM/ikqgilT4K23\n4IwzdAeriEQi6U1M7r4Q6F3H63fUeHwLcEtmS8thBx8Mv/kNjBsHN94IF18cdUUikmc0t0xjuegi\n+M//DFfPaB1WEWliCvfGYgaTJ4e54EeM0A1OItKkFO6NaeedwxTBH34Ip50GW7ZEXZGI5AmFe2Pr\n0yeMu0+fDjfcEHU1IpInFO5N4ZxzwqWRl18Or74adTUikgcU7k3BDO6+Gzp3hpEj4bPPoq5IRGJO\n4d5Udt4ZHnwQ3n8fRo/W9e8i0qgU7k3pwANhwgR45JEw0ZiISCNRuDe1cePgiCPgwgth8eLk7UVE\n0qBwb2rNmsF994VhmhEjaPbll1FXJCIxpHCPQlER/PGPsGQJXW6/PepqRCSGFO5ROfJIuOgiOj7+\neFjJSUQkgxTuUfqf/+Hz0tIwe+T770ddjYjEiMI9Sttvz9KrroJ//QtOPVXTE4hIxijcI/ZF585w\n001h/VVNTyAiGaJwzwZnnhmmB77ySpg7N+pqRCQGFO7ZwAwmTYLddgvTE2zYEHVFIpLjFO7Zom1b\nuP9+WLkSfvazqKsRkRyXygLZO5jZG2a2wMyWmNk1dbQxM5toZqvMbKGZ9WmccmOuvBwuvTQcxU+d\nGnU1IpLDUjly/xI41N33B3oBR5lZ/1pthgGliW0MoDtz0jV+fJgD/qyz4O9/j7oaEclRScPdg6rE\n04LEVntKw+OA+xJtXwPamFmHzJaaJ7bbDh54IFweqdWbRCRNKY25m1lzM5sPfAg85+6v12rSEXi3\nxvP3Eq9JOrp1gz/8AZ5/Hn7/+6irEZEcZF6PecXNrA0wFfipuy+u8fo04Dp3n5V4PgO4zN1n1/r8\nGMKwDUVFRWUVFRVpFV1VVUVhYWFan802W+2LOz2uvprvvf46c2+7jaquXZu+uHrKi/2SY+LSD1Bf\nqpWXl89x975JG7p7vTbgamBcrdfuBEbWeL4C6LCt7ykrK/N0VVZWpv3ZbLPNvnz0kXuHDu777OO+\nYUOT1ZSuvNkvOSQu/XBXX6oBsz2FrE7lapldE0fsmFlLYCiwvFazJ4BTE1fN9AfWufvaVH8TyVa0\naxemB162LFxFIyKSolTG3DsAlWa2EHiTMOY+zczGmtnYRJvpwBpgFTAJOKdRqs1Hhx8OF18Mt94K\n06ZFXY2I5IgWyRq4+0Kgdx2v31HjsQPnZrY0+caECTBjBowaBYsWQfv2UVckIllOd6jmgu23D5dH\nbtgAp5+uyyNFJCmFe67YZx+48UZ45hmYODHqakQkyyncc8lPfgLHHguXXQYLFkRdjYhkMYV7LjGD\nu+8Ok4z96EewcWPUFYlIllK455pddw2XRy5dCuPGRV2NiGQphXsuGjo0BPvtt8Njj0VdjYhkIYV7\nrpowIcweeeaZWlxbRP6Nwj1X1Zw98pRTYPPmqCsSkSyicM9l3brBzTdDZSX89rdRVyMiWUThnutG\njYKTT4b//m947bWoqxGRLKFwz3VmcMcdUFwcFtdety7qikQkCyjc46BNmzD+/u674UaneszRLyLx\npHCPi4MOCuuvTpkC99wTdTUiEjGFe5z8/OdhiuDzz4clS6KuRkQipHCPk2bN4P77YaedwklWTU8g\nkrcU7nHTvn0I+CVL4MILo65GRCKicI+jI44IQzSTJsGDD0ZdjYhEQOEeV+PHw8CBMGYM/PWvUVcj\nIk1M4R5XBQXhqH377eGkk8I0BSKSN5KGu5kVm1mlmS01syVmdkEdbYaY2Tozm5/Yrm6ccqVeiovD\n9MALFmj8XSTPJF0gG9gEXOzuc81sJ2COmT3n7ktrtXvZ3Y/OfInSIMOHw6WXhrlnBg8Oi3yISOwl\nPXJ397XuPjfx+HNgGdCxsQuTDPrVr2DQoDD+vmxZ1NWISBOo15i7mZUAvYHX63h7gJktNLOnzKxH\nBmqTTCkogIoK2HFHOOEE2LAh6opEpJGZpzgPiZkVAi8CE9z90VrvtQa2uHuVmQ0HbnL30jq+Ywww\nBqCoqKisoqIiraKrqqooLCxM67PZpin7ssucOex3ySV8cPjhLL/88jDpWAZpv2SfuPQD1Jdq5eXl\nc9y9b9KG7p50AwqAZ4CLUmz/NtBuW23Kyso8XZWVlWl/Nts0eV+uucYd3G+/PeNfrf2SfeLSD3f1\npRow21PI4VSuljHgHmCZu9+4lTbtE+0ws36E4Z5PUvktJE3sqqtg2LAw/8zrdY2uiUgcpHK1zEDg\nFGCRmc1PvHYF0BnA3e8ATgDONrNNwBfAiMRvGMk2zZrBn/4EZWVh/H3uXNh116irEpEMSxru7j4L\n2ObgrLvfAtySqaKkkbVtC488AgMGhAU+nnkGmjePuioRySDdoZqv+vSBW2+FGTPgyiujrkZEMkzh\nns/OPDOs3PSb38DDD0ddjYhkkMI93910U1jFadQoWLw46mpEJEMU7vlu++3hz38OC3wcfzx89lnU\nFYlIBijcBXbfPZxgfecdGDECNm2KuiIRaSCFuwQDBoQTrM8+Gxb6EJGclsp17pIvRo8O0wPfcAPs\nvz+cckrUFYlImnTkLt/1+9/DkCEh6HUHq0jOUrjLdxUUhMsid989nGB9992oKxKRNCjc5d+1awdP\nPhmmBj72WE0RLJKDFO5Stx49YMoUWLgwjL1v2RJ1RSJSDwp32bphw8LJ1alT4Yoroq5GROpBV8vI\ntl1wAaxYEaYo6NoVzjor6opEJAUKd9k2M7j5Znj7bRg7Fr7/fRg6NOqqRCQJDctIci1ahPH37t3D\nHPCag0Yk6yncJTWtW8P//R+0agXDh8P770ddkYhsg8JdUldcDNOnh8nFhg+H9eujrkhEtkLhLvXT\nq1eYZGzpUvjhD+Grr6KuSETqkMoC2cVmVmlmS81siZldUEcbM7OJZrbKzBaaWZ/GKVeywhFHwN13\nw/PPwxln6Bp4kSyUytUym4CL3X2ume0EzDGz59x9aY02w4DSxHYgcHviX4mr004L4+5XXgm77Rau\nhxeRrJHKAtlrgbWJx5+b2TKgI1Az3I8D7nN3B14zszZm1iHxWYmryy+HDz4Ik421bw/9+kVdkYgk\n1GvM3cxKgN5A7ekCOwI1Z5h6L/GaxJlZCPaRI+Gyy2j/1FNRVyQiCSnfxGRmhcAjwIXuntZlEmY2\nBhgDUFRUxMyZM9P5GqqqqtL+bLaJQ19s1Cj2XbmSbtdfz+Idd+TjQw6JuqQGi8N+gfj0A9SXenP3\npBtQADwDXLSV9+8ERtZ4vgLosK3vLCsr83RVVlam/dlsE5u+VFX5P3v0cC8ocH/66airabC47Je4\n9MNdfakGzPYUcjuVq2UMuAdY5u43bqXZE8Cpiatm+gPrXOPt+aVVKxZdd12YTfI//gNeeSXqikTy\nWipj7gOBU4BDzWx+YhtuZmPNbGyizXRgDbAKmASc0zjlSjbbVFgIzzwTbnYaPhxmz466JJG8lcrV\nMrMAS9LGgXMzVZTksN12gxkzYPDgcD18ZWVYj1VEmpTuUJXM69QJXngBCgvh8MPD3awi0qQU7tI4\nSkrCEXxBARx6KCxbFnVFInlF4S6Np7Q0HMGbQXm5Al6kCSncpXHtvXcYdwcFvEgTUrhL49t7b6i+\nYWPIEC32IdIEFO7SNKoDvkWLEPDz5kVdkUisKdyl6ey9N7z0UljN6dBD4Y03oq5IJLYU7tK0unQJ\nAd+2bbhM8sUXo65IJJYU7tL0vv/9EPDFxXDUUWFtVhHJKIW7RKNjx3DU3rMnHH88PPhg1BWJxIrC\nXaLTrl240WngQPjxj+GWW6KuSCQ2FO4Srdat4amn4Ljj4Kc/hauugjBttIg0gMJdoteyJTz8MIwe\nDRMmhH83bYq6KpGclvJKTCKNqkULuPPOsBbrtdfC2rUwZUqYfExE6k1H7pI9zGD8+BDyzzwDhxwS\nQl5E6k3hLtlnzBh44glYsQL699d0BSJpULhLdho+PFwq+fXXMGBAOOkqIilTuEv2KisLUxR06QJH\nHw033xx1RSI5Q+Eu2a1TJ3j55RDu558PY8fCV19FXZVI1ksa7mY22cw+NLM6Bz7NbIiZrauxePbV\nmS9T8lphITz6KPz85+Fk6+GHw0cfRV2VSFZL5cj9XuCoJG1edvdeiW18w8sSqaV5c/j1r+GBB+DN\nN6FvX00bLLINScPd3V8CPm2CWkSSGzkSZs2CLVvCidb77ou6IpGslKkx9wFmttDMnjKzHhn6TpG6\nlZXBnDlw0EFw2mlw3nkahxepxTyFeTzMrASY5u4963ivNbDF3avMbDhwk7uXbuV7xgBjAIqKisoq\nKirSKrqqqorCmNy5qL6kzzZvZo9Jk+g8ZQrr99mHJb/4BV8WFWXku+OyX+LSD1BfqpWXl89x975J\nG7p70g0oARan2PZtoF2ydmVlZZ6uysrKtD+bbdSXDHj4YfeddnJv29Z9+vSMfGVc9ktc+uGuvlQD\nZnsKWdzgYRkza29mlnjcjzDU80lDv1ckZSecEIZpOnUKNz9dfnm4+Ukkj6VyKeSDwKtANzN7z8zO\nNLOxZjY20eQEYLGZLQAmAiMSv11Emk5pKbz2WphR8rrrwrw0b78ddVUikUk6K6S7j0zy/i2AVlmQ\n6LVsCXfdBYcdFuan6dULJk2CE0+MujKRJqc7VCV+Tj45XAPfrRucdBKMGgWffx51VSJNSuEu8bTn\nnuF6+KuuCtfC9+oFr74adVUiTUbhLvFVUBAW/njxRdi8GQYNgiuu0DXxkhcU7hJ/gwbBwoVw+ulh\nCoMDDoAFC6KuSqRRKdwlP7RuDffcExYB+eCDMDfNL3+po3iJLYW75JdjjoElS8JJ12uuCUfxc+dG\nXZVIxincJf9873vwpz/B44+HqYP79YNLLoGNG6OuTCRjFO6Sv449FpYuhTPOgOuvh5494bnnoq5K\nJCMU7pLf2rQJNz7NnBmurjniCBg5ku0+0QwaktsU7iIQpitYsCCcZJ06lX6nnRbWbN20KerKRNKi\ncBeptsMO8ItfwKJFrN9nn7Bma1lZWMNVJMco3EVqKy1l4W9/C3/+M3z2GQweDD/+Mbz3XtSViaRM\n4S5SFzP44Q9h+fIwhcEjj4S5asaP11U1khMU7iLbsuOOYQqD5cvh6KPDsE23bmG+mi1boq5OZKsU\n7iKpKCmBKVPgpZegQ4ewdmtZGTz/fNSVidRJ4S5SHwcfHBYFeeCBMB4/dGi4fFJ3uUqWUbiL1Fez\nZjByZBiqueGGsMRfWRmMGAErVkRdnQigcBdJ3w47wEUXwZo14aTrtGnQvXuYfXLNmqirkzyXyhqq\nk83sQzNbvJX3zcwmmtkqM1toZn0yX6ZIFtt553DSdc0auPDCMDbfrRucdZZCXiKTypH7vcBR23h/\nGFCa2MYAtze8LJEctNtuYZhm9WoYOzZMTrbXXmHumlWroq5O8kzScHf3l4BPt9HkOOA+D14D2phZ\nh0wVKJJzdt89TF2wejWce244+dqtW7gRanGdfwCLZFwmxtw7Au/WeP5e4jWR/NaxI9x0E7z9Nlx8\ncVgoZN99w2yUr7wSdXUSc+buyRuZlQDT3L1nHe9NA65z91mJ5zOAy9x9dh1txxCGbigqKiqrqKhI\nq+iqqioKCwvT+my2UV+yU2P0pcX69XScOpVOjz5Kwfr1rOvRg3dHjODjgw6C5s0z+rOqaZ9kp4b0\npby8fI67903a0N2TbkAJsHgr790JjKzxfAXQIdl3lpWVeboqKyvT/my2UV+yU6P2parK/eab3UtK\n3MG9a1f3W24Jr2eY9kl2akhfgNmeQm5nYljmCeDUxFUz/YF17r42A98rEk+tWsF558HKlfDQQ2Fl\nqPPOg06dYNy4MIwj0kCpXAr5IPAq0M3M3jOzM81srJmNTTSZDqwBVgGTgHMarVqROGnRAk48EV59\nNYzBH3EE/OEP0KULHH98WBUqhWFTkbq0SNbA3Ucmed+BczNWkUi+MYMBA8L27rtw221w991hjde9\n9oKzz4ZTT4W2baOuVHKI7lAVySbFxfDrX4eQv//+EOg/+1m48ubUU8MRvo7mJQUKd5FstMMO8F//\nFYZs5s+HUaPgscdg0KAwxcH118OHH0ZdpWQxhbtIttt//zBU8/e/w+TJ4Wj+kkvC0fzxx4fhm6+/\njrpKyTIKd5FcUVgYjuBfeQWWLAnz2Lz2Wgj4jh3Dmq+zZ2vYRgCFu0hu6t4dfve7sK7rk0/CkCFw\n111wwAHhvWuv1Xw2eU7hLpLLWrQIy/899BD84x8h4IuK4OqrobSUPmefDTfeqMW985DCXSQu2rSB\n0aNh5kx45x343e+wzZvDvDbFxWEVqYkT4f33o65UmoDCXSSOioth3Djm3HUX/PWvYZjmn/+ECy4I\nd8IOHBimJ9Z887GlcBeJu9LSsFLUokVhacBf/Qo2bgxTHXTpAr16wS9/CfPm6WRsjCjcRfJJt25w\n5ZUhyNesCUfvhYUwfjz06QMlJWEO+qefhn/9K+pqpQEU7iL5ao89whqws2aFk7GTJ0Pv3nDvvTBs\nGLRrB8cdB3feGe6YlZySdG4ZEckDu+0WrqEfNQq++CKclH3ySZg+PSwyAtCzJxx5JBx1VLhTdocd\nIi1Ztk3hLiLf1bJlOHIfNiyMwS9bFkL+6afD8oE33BDaDB4MQ4eGrWdPaKaBgGyicBeRrTMLN0V1\n7x5OwG7YEI7qn302TEk8blxot+uucOihcNhhUF4eTtSaRVp6vlO4i0jqWrWCH/wgbBDG4mfM+Hab\nMiW8XlwcQv6QQ8K2554K+yamcBeR9BUXw+mnh80dVqyAykp44QV46im4777QrmPHcBPV4MHh3+7d\nNYzTyBTuIpIZZrD33mE7++xvx+tffBFeeilsFRWhbZs2YXGSgQPDvwccEP4qkIxRuItI46g5Xl8d\n9m+9BS+/HGa2nDUrnKgFaN48TG180EHQvz8ceCB07aqhnAZQuItI0zALY+977gmnnRZe+/TTMG3x\nX/4Stj/+EW69Nby3yy7Qr1/YDjiA7TRnfb2kFO5mdhRwE9AcuNvdr6v1/hDgceCtxEuPuvv4DNYp\nInHUti0MHx42gM2bYenSEPhvvBG2CRNgyxYGQBi779sXysq+3YqKouxB1koa7mbWHLgVGAq8B7xp\nZk+4+9JaTV9296MboUYRyRfNm8O++4Zt9Ojw2oYNMG8eqyoq6PrZZzBnTlh9qtruu4c7a3v3DvPk\n9O4dplHI8xO2qRy59wNWufsaADOrAI4Daoe7iEjmtWoFgwbx3qZNdB0yJLy2fn2YH2fePJg7N2xP\nPx2O/AF22gn22y+M4++3X9h69gyv5wnzJLPAmdkJwFHuflbi+SnAge5+Xo02Q4BHCUf27wPj3H1J\nHd81BhgDUFRUVFZRfea8nqqqqigsLEzrs9lGfclOcelLXPoByfvS7MsvafXWWxSuXEnhmjW0Wr2a\nwjVraLFhwzdtvmjfng177smGkhI2lJSwcY892Ni5M1u2264puvCNhuyX8vLyOe7eN1m7TJ1QnQt0\ndvcqMxsOPAaU1m7k7ncBdwH07dvXh1T/Fq6nmTNnku5ns436kp3i0pe49APS7Is7/O1vYbrjhQtp\nuWgRLRcvpt1DD8GmTaFNs2bhjtoePcKVPfvsE7Zu3cKMmY2gKfZLKuH+PlBc43mnxGvfcPf1NR5P\nN7PbzKydu3+cmTJFRNJgFsbfS0rgmGO+ff2rr8IiJkuWhG3p0rBNm/Zt6EO4Sav62v1u3cK2115h\nwZMsH9NPJdzfBErNbA9CqI8AflSzgZm1Bz5wdzezfoSphD/JdLEiIhmx3XZhDL5nz+++/tVXYWHx\n5cvDDVjLl4e7bu+9Fz7//Nt2LVuG6/BLS0PYVz/u2hU6dMiK6/OThru7bzKz84BnCJdCTnb3JWY2\nNvH+HcAJwNlmtgn4AhjhyQbzRUSyzXbbfXvjVU3usHZtONpfsSJsK1eGo/4nn4Sa1+DvuGO4lr9r\n1zDc06VLeN6lC3TuHH5GE0hpzN3dpwPTa712R43HtwC3ZLY0EZEsYRYuudx9d6g9Vr5pU1iQfOXK\nsK1eHbYVK8L8Ol9++W3bZs2gUyc6/eAH//49GaY7VEVEGqJFi2/vvD3yyO++t2VLOOJfvTpMvbBm\nDbz1Fl9973uNX1aj/wQRkXzVrFm4q7ZjxzAjZsKHM2fSfRsfy8iPbuTvFxGRCCjcRURiSOEuIhJD\nCncRkRhSuIuIxJDCXUQkhhTuIiIxpHAXEYmhpPO5N9oPNvsI+FuaH28HxGXGSfUlO8WlL3HpB6gv\n1b7v7rsmaxRZuDeEmc1OZbL6XKC+ZKe49CUu/QD1pb40LCMiEkMKdxGRGMrVcL8r6gIySH3JTnHp\nS1z6AepLveTkmLuIiGxbrh65i4jINmR1uJvZUWa2wsxWmdnP63jfzGxi4v2FZtYnijpTkUJfhpjZ\nOjObn9iujqLOZMxsspl9aGaLt/J+Lu2TZH3JlX1SbGaVZrbUzJaY2QV1tMmJ/ZJiX3Jlv+xgZm+Y\n2YJEX66po03j7Rd3z8qNsF7ramBPYDtgAdC9VpvhwFOAAf2B16OuuwF9GQJMi7rWFPoyGOgDLN7K\n+zmxT1LsS67skw5An8TjnYC/5vD/K6n0JVf2iwGFiccFwOtA/6baL9l85N4PWOXua9z9K6ACOK5W\nm+OA+zx4DWhjZh2autAUpNKXnODuLwGfbqNJruyTVPqSE9x9rbvPTTz+HFgGdKzVLCf2S4p9yQmJ\n/9ZViacBa1g5AAAB5klEQVQFia32Sc5G2y/ZHO4dgXdrPH+Pf9/JqbTJBqnWOSDxp9lTZtajaUrL\nuFzZJ6nKqX1iZiVAb8JRYk05t1+20RfIkf1iZs3NbD7wIfCcuzfZftEaqtljLtDZ3avMbDjwGFAa\ncU35Lqf2iZkVAo8AF7r7+qjraYgkfcmZ/eLum4FeZtYGmGpmPd29znM8mZbNR+7vA8U1nndKvFbf\nNtkgaZ3uvr76Tzh3nw4UmFm7pisxY3JlnySVS/vEzAoIYfj/3P3ROprkzH5J1pdc2i/V3P2fQCVw\nVK23Gm2/ZHO4vwmUmtkeZrYdMAJ4olabJ4BTE2ec+wPr3H1tUxeagqR9MbP2ZmaJx/0I++aTJq+0\n4XJlnySVK/skUeM9wDJ3v3ErzXJiv6TSlxzaL7smjtgxs5bAUGB5rWaNtl+ydljG3TeZ2XnAM4Sr\nTSa7+xIzG5t4/w5gOuFs8ypgIzAqqnq3JcW+nACcbWabgC+AEZ44nZ5NzOxBwtUK7czsPeAXhBNF\nObVPIKW+5MQ+AQYCpwCLEuO7AFcAnSHn9ksqfcmV/dIB+KOZNSf8AnrI3ac1VYbpDlURkRjK5mEZ\nERFJk8JdRCSGFO4iIjGkcBcRiSGFu4hIDCncRURiSOEuIhJDCncRkRj6/7sXSfNXP+hvAAAAAElF\nTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "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": 26, "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": 5, "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": 40, "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", "を通る多項式を(a)ラグランジュ補間, (b)逆行列で求めよ．" ] }, { "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": 2, "metadata": { "collapsed": true }, "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": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHtFJREFUeJzt3XucVXW9//HXZxgQ5SagjAoIyAxyU4ThZt4glQYB0bIf\n3jU9YppaWsc0y8z0yDniHZKM0NSSNM0LktpPQTuGCIgphiii3JJIRWS8AfI5f3z3jj0jw2yYvWet\ntff7+Xish6y19+z5LHHm7fpezd0REZHiUxJ1ASIiEg0FgIhIkVIAiIgUKQWAiEiRUgCIiBQpBYCI\nSJFSAIiIFCkFgIhIkVIAiIgUqdKoC9iePfbYw7t27Vrv+z7++GNatGiR/4LyKOn3oPqjl/R7UP25\ns2DBgvfcfc/63hfLADCzMcCY8vJy5s+fX+/7Z8+ezbBhw/JeVz4l/R5Uf/SSfg+qP3fMbHk274tl\nE5C7P+bu49u0aRN1KSIiBSuWASAiIvmnABARKVIKABGRIqUAEBEpUgoAEZEipQAQESlShRkADz8M\nN9wQdRUiIrFWmAHwxBNw3XVRVyEiEmuxDAAzG2Nmd6xfv37nPqCiAt5/PxwiIrJNsQyABs8E7tEj\n/PPNN3NXlIhIgYllADSYAkBEpF6FGQDdukFJCbzxRtSViIjEVmEGQLNmIQQUACIidSrMAIDQDKQm\nIBGROhVuAFRUhCcA96grERGJpcINgB494OOPYc2aqCsREYmlwg2AiorwT/UDiIhsU+EGQHooqAJA\nRGSbCjcAOneGXXZRR7CISB0KNwCaNIHu3fUEICJSh8INAAjNQAoAEZFtatQAMLPjzOxXZvZ7MxuR\n92/Yowe89RZ88UXev5WISNJkHQBmNs3M1prZolrXq8xsiZktNbPLtvcZ7v6wu58DfBsYt3Ml74CK\nCti4EVasyPu3EhFJmh15ArgLqMq8YGZNgMnASKA3cJKZ9TazA8xsRq2jQ8aX/jj1dfmlReFEROpU\nmu0b3f05M+ta6/JgYKm7LwMws+nAWHe/Dhhd+zPMzIAJwJ/c/aWdLTprmXMBRuS/xUlEJEmyDoA6\ndARWZpyvAoZs5/0XAkcBbcys3N2n1H6DmY0HxgOUlZUxe/bseouorq7e9vvcOXTXXVnzzDMs7du3\n3s+JUp33kBCqP3pJvwfVHwF3z/oAugKLMs5PAKZmnJ8GTNqRz9zeUVlZ6dmYNWtW3S8OGOBeVZXV\n50Rpu/eQAKo/ekm/B9WfO8B8z+J3bENHAa0GOmecd0pda5AGbwmZKb0onIiI1NDQAJgHVJhZNzNr\nBpwIPNrQoryhW0Jm6tED3nknjAYSEZF/25FhoPcBc4D9zWyVmZ3t7puBC4AngcXA/e7+Wn5K3Uk9\nesCWLbBsWdSViIjEyo6MAjqpjuszgZk5q4jQBASMKS8vb/iHZY4E6tmz4Z8nIlIgYrkURE6bgNIB\noLkAIiI1xDIAcqpdO2jfXh3BIiK1xDIAcjoKCLQonIjINsQyAHLaBATaIF5EZBtiGQA5V1EBq1eH\nPYJFRAQolgBILwq3dGm0dYiIxEgsAyAvfQCgfgARkQyxDICc9wGk5xMoAERE/i2WAZBzLVpAx44K\nABGRDMURAAD9+sGcOVFXISISG7EMgJz3AQCMHBmGguopQEQEiGkA5LwPAGDUqPDPxx/P3WeKiCRY\nLAMgL7p1g169FAAiIinFEwAQngKeew42bIi6EhGRyBVfAGzaBH/+c9SViIhELpYBkJdOYIBDDoE2\nbdQMJCJCTAMgL53AAE2bwogRMHNm2CVMRKSIxTIA8mrUKFizBhYujLoSEZFIFV8AjBwJZmoGEpGi\nV3wB0KEDDBqkABCRold8AQChGWjePFi7NupKREQiU7wB4A5/+lPUlYiIRCaWAZC3YaBp/fvD3nur\nGUhEilosAyBvw0DTSkrgmGPgySfDxDARkSIUywBoFKNGwUcfwfPPR12JiEgkijcAjjoqTAxTM5CI\nFKniDYBWreCIIxQAIlK0ijcAIDQDLV4Mb78ddSUiIo1OAQB6ChCRolTcAVBREQ4FgIgUoeIOAAhP\nAbNmwccfR12JiEijUgCMGgWffw7PPBN1JSIijSqWAZD3mcCZDj8cWrZUM5CIFJ1YBkDeZwJnatYM\njj46BIB7/r+fiEhMxDIAGt2oUbBqFbzyStSViIg0GgUAhHWBQM1AIlJUFAAQVgYdMEABICJFRQGQ\nNno0vPACvP9+1JWIiDQKBUDaqFGwZQs88UTUlYiINAoFQNrAgWG/YDUDiUiRUACklZTAyJHhCWDz\n5qirERHJOwVAplGjYN06+Otfo65ERCTvFACZvvY1aNsWrr5ak8JEpOApADK1bg3XXANPPw0PPRR1\nNSIieaUAqG38eDjwQLjkEvjkk6irERHJm0YLADPrZWZTzOwPZnZeY33fHVZaCrfeCitWwPXXR12N\niEjeZBUAZjbNzNaa2aJa16vMbImZLTWzy7b3Ge6+2N2/Dfw/4JCdL7kRHHEEjBsHEybA8uVRVyMi\nkhfZPgHcBVRlXjCzJsBkYCTQGzjJzHqb2QFmNqPW0SH1NccCjwMzc3YH+XL99WAGP/hB1JWIiOSF\neZajXcysKzDD3fumzg8GrnL3r6XOLwdw9+uy+KzH3X1UHa+NB8YDlJWVVU6fPr3e2qqrq2nZsmVW\n97EjutxzD92mTePlG27gwwEDcv75mfJ1D41F9Ucv6feg+nNn+PDhC9x9YL1vdPesDqArsCjj/ARg\nasb5acCk7Xz9MOBW4JfAd7L5npWVlZ6NWbNmZfW+Hfbpp+7durn36eO+aVN+vkdK3u6hkaj+6CX9\nHlR/7gDzPYvfsY3WCezus939Inc/190nN9b3bZDmzeGmm+C11+D226OuRkQkpxoSAKuBzhnnnVLX\nGqxRt4Ssz7HHwogRcOWV8K9/RV2NiEjONCQA5gEVZtbNzJoBJwKP5qIob8wtIetjBjffDNXVcMUV\nUVcjIpIz2Q4DvQ+YA+xvZqvM7Gx33wxcADwJLAbud/fX8ldqhHr1gosugqlTYcGCqKsREcmJ0mze\n5O4n1XF9JnkY0mlmY4Ax5eXluf7onXfllXDvvXDhhfD88+HJQEQkwWK5FESsmoDS2rQJE8PmzIHf\n/jbqakREGiyWARBbZ5wBgwbBpZfChg1RVyMi0iCxDIBYjQLKVFICt90G774L114bdTUiIg0SywCI\nZRNQ2pAhcOaZcOON8OabUVcjIrLTYhkAsXfddWGS2MUXR12JiMhOi2UAxLYJKG2vveCnPw0byGsT\neRFJqFgGQKybgNIuvBB69oTvfQ8+/zzqakREdlgsAyARmjWDW26BpUvDTGERkYRRADTEiBEwdiz8\n/Ofwj39EXY2IyA5RADTUjTfC5s3wwx9GXYmIyA6JZQDEvhM40377hV3D7r03LBEhIpIQsQyARHQC\nZ7r8cujYMSwY98UXUVcjIpKVWAZA4rRoARMnwksvwbRpUVcjIpIVBUCujBsHhx0GP/oRrFsXdTUi\nIvVSAOSKGdx6K3zwAVx1VdTViIjUK5YBkKhO4EwHHQTnnguTJ8OiRVFXIyKyXbEMgMR1Amf6+c/D\n3gEXXQTuUVcjIlKnWAZAorVvD9dcA7NmwYMPRl2NiEidFAD5MH489OsH3/8+fPJJ1NWIiGyTAiAf\nmjQJHcIrVoRRQWoKEpEYUgDky+GHw3nnhQXjvvUt2Lgx6opERGoojbqAgjZ5MpSVhWGh77wDDz0E\n7dpFXZWICBDTJ4DEDgOtzSxsHHPvvTBnDgwdqm0kRSQ2YhkAiR4Gui2nnAJPPx0miQ0dCn/5S9QV\niYjEMwAK0qGHwty5sOeecOSRcM89UVckIkVOAdCYuncPTUGHHgqnnw5XXqkRQiISGQVAY2vbFp54\nAs46K8waPvlk+OyzqKsSkSKkAIhCs2YwdSpMmADTp8NXv0rTDz+MuioRKTIKgKiYhW0kH3gAFi5k\nwPnnw+LFUVclIkVEARC1E06AZ5+lyWefwcEHh9FCIiKNQAEQB4MHs+AXv4DOnaGqCn71q6grEpEi\noACIic/32itsKn/UUWExuUsvhS1boi5LRApYLAOgYGYC76jWreGxx8IaQtdfH5qHtJqoiORJLAOg\n4GYC74jS0rCG0E03wcMPwxFHwLvvRl2ViBSgWAZA0TOD730PHnkkjAwaMgReeSXqqkSkwCgA4mzM\nmLBu0JYtcMghMHNm1BWJSAFRAMRd//5hDaGKihAIkyZFXZGIFAgFQBJ07AjPPQejR8OFF4YN57/4\nIuqqRCThFABJ0bJl2FDmkkvgtttg7FjYsCHqqkQkwRQASdKkCdxwA9x+e1hQ7rDDYOXKqKsSkYRS\nACTRt78dOoTffjuMEFqwIOqKRCSBFABJNWJEmDncrFl4Enj44agrEpGEUQAkWd++YYTQgQfC178O\nEydqgxkRyZoCIOnKymDWrLBsxH/+Z2ge2rQp6qpEJAEaNQDMrIWZzTez0Y35fQverruGjWV+9CO4\n4w445hjQBjMiUo+sAsDMppnZWjNbVOt6lZktMbOlZnZZFh/1Q+D+nSlU6lFSAtdeC3feCc8+C1/5\nSugkFhGpQ7ZPAHcBVZkXzKwJMBkYCfQGTjKz3mZ2gJnNqHV0MLOjgb8Da3NYv9R25pnw1FOwZk0Y\nITRnTtQViUhMZRUA7v4c8EGty4OBpe6+zN03AtOBse7+qruPrnWsBYYBQ4GTgXPMTP0P+TJsWPjF\n37o1DB8Ov/991BWJSAyZZzlqxMy6AjPcvW/q/ASgyt3/I3V+GjDE3S+o53POBN5z9xl1vD4eGA9Q\nVlZWOX369Hprq66upmXLllndR1zl4x6arl9PnyuvZPdXXuHts85i+amnhpVG8yDpfwdJrx+Sfw+q\nP3eGDx++wN0H1vtGd8/qALoCizLOTwCmZpyfBkzK9vOyOSorKz0bs2bNyup9cZa3e/jsM/dTT3UH\n99NPD+d5kPS/g6TX7578e1D9uQPM9yx+x5Y2IGRWA50zzjulrkmc7LIL3H039OgBV14J77wT1hRq\n3z7qykQkYg1ph58HVJhZNzNrBpwIPJqLoop2S8h8MYOf/AR+9zt44QUYOhTeeCPqqkQkYtkOA70P\nmAPsb2arzOxsd98MXAA8CSwG7nf313JRlBfzlpD5dNJJ8MwzYY7A0KFw882gkBUpWtmOAjrJ3fd2\n96bu3sndf526PtPde7h7d3e/NldF6Qkgjw45JCwf0acPXHxx2GvgO98JW0+KSFGJ5VBMPQHk2X77\nha0m588PS0hMnQq9e8PRR8Ojj2qzGZEiEcsAkEZSWQl33RX2FLjmmvAUMHZs2H5y4kRYty7qCkUk\njxQAAh06wBVXhKUj7r8fOnUKC8t16gTnnguLFtX/GSKSOLEMAPUBRKRpU/jmN8P+wwsXhk7ju++G\nAw4IM4ofegg2b466ShHJkVgGgPoAYuCgg0LfwKpVMGECLFsG3/gGdO8O//3f8P77UVcoIg0UywCQ\nGGnfHn74Q3jrrfAE0L07XHZZaB46+2x4+eWoKxSRnRTLAFATUAyVlsLxx4d5BK++CmecEfYg6N8f\nDjuMPWfN0kY0IgkTywBQE1DM9e0LU6aE5qGJE2H1avpcfTV06xb2JFirFb9FkiCWASAJ0bYtfP/7\n8OabvHrttWEuwY9/DJ07hyeE+fOjrlBEtkMBIA3XpAnvf+UrYSOaxYvhnHPgwQdh0CA4+OCwBtHG\njVFXKSK1KAAkt3r2hEmTYPXqsNbQe+/BKadAly7ws5+FncpEJBZiGQDqBC4AbdrAd78LS5bAzJmh\ns/iqq2DffUMgvPACZLkZkYjkRywDQJ3ABaSkBEaODCGwZAmcdx489lhoGho8GO65Bz7/POoqRYpS\nLANAClSPHnDLLaF5aNIkqK6G008PTwU/+Um4LiKNRgEgja9Vq7AE9d//HjqOhwwJw0e7doVx4+B/\n/1fNQyKNQAEg0THbugT10qVw0UXw5JNw2GEwYADceSd8+mnUVYoUrFgGgDqBi9B++8ENN4RmoClT\nwqzis84KcwouvxxWrIi6QpGCE8sAUCdwEWvRIixB/eqrYdmJww+H//mfMMv4G9+A2bPVPCSSI7EM\nABHMti5BvWwZ/OAH4Zf/8OHQrx/86lfwySdRVymSaAoAib8uXcIS1CtXhiWqS0pg/PitG9e8807U\nFYokkgJAkmO33cIS1AsXhk1rjjwSbrop9B8cdxw8/bSah0R2gAJAkscsjBR64IGwjeXll8Pzz8NR\nR4WVSm+/PcwxEJHtUgBIsnXuHOYQrFwZNrjfdVc4//zQPHTxxWF4qYhsUywDQMNAZYc1bx6WoJ43\nD/7617D8xKRJYfbxqFFhfsGWLVFXKRIrsQwADQOVnWYW1hm67z5YvjwsMbFgAVRVQa9ecNtt8NFH\nUVcpEguxDACRnNhnn7AE9fLlcO+9YQObiy4KzUMXXhgWpxMpYgoAKXy77LJ1Ceq5c2HsWPjlL8Pe\nBVVV8Pjjah6SoqQAkOKSXoJ65Uq4+mp45RUYPZohp50WhpR++GHUFYo0GgWAFKeystA/sHw5TJ/O\nxrZt4ZJLQvPQ+eeHlUpFCpwCQIpb06YwbhwLJ00KncXf/CZMmwZ9+oR5BY88Al98EXWVInmhABBJ\nSy9BvXIl/Nd/hU7i446D8nK4/nr44IOoKxTJKQWASG177hlmF7/9dphtvO++cOmloXno+ONhwoSw\nUqmGk0rClUZdwLaY2RhgTHl5edSlSDErLYUTTgjH3/4Wlph45hl4+OHwulmYWzBkSOhcHjIkLEXR\ntGm0dYtkKZYB4O6PAY8NHDjwnKhrEQHCEtRTpoQ/f/ABvPhiOObODZvc33lneG3XXUNTUmYodOkS\nwkIkZmIZACKx1q5dmD9QVRXO3UNzUToQ5s6FyZPhxhvD6x06bA2DwYPDsfvu0dUvkqIAEGkos7Ak\n9X77wYknhmubNoU5BulQePFFmDFj69fsv3/NUOjXD5o1i6Z+KVoKAJF8aNoUKivDcd554dr69WGx\nunQoPPVUmJQG4Zd///41m466d1fTkeSVAkCksbRpE+YWHHVUOHcPQ07TTwhz54Ydz269Nbzert2X\nm4722CO6+qXgKABEomIWhpjuu2+YgAaweTO89lrN/oQnn9y609l++4VASIdC//5hKWyRnaAAEImT\n0tLQH9CvH5yTGgS3YUOYpZwOhb/8JSx3DaGpqV8/GDyYstatYa+9wh4IJZriI/VTAIjEXatWMGxY\nONL+8Y+aTUd3302v6uowSa1NGxg0qGZ/QllZVNVLjCkARJJon33CrOTjjw/nX3zBi3ffzWDY2nQ0\nYcLWdYy6dKnZn1BZCbvtFlX1EhMKAJFC0KQJn3TrFp4SvvWtcO2TT+Cll2oORX3ggX+/n759a/Yn\n9OoVrkvRUACIFKrddoNDDw1H2j//WXMW8+9/D3fcEV5r2TI0HWU+KXTsGE3t0igUACLFpKwMxowJ\nB4Sd0N58s+ZTwo03holsEAIgMxAGDgx9ElIQGi0AzGwY8HPgNWC6u89urO8tInUoKQmzkvffH047\nLVz77DN4+eWaQ1H/+MfwmlnYKyEzFPr2DaOXJHGy+lszs2nAaGCtu/fNuF4F3AI0Aaa6+4TtfIwD\n1UBzYNVOVywi+dW8OQwdGo6099+v2XT0yCNh4xwIC+BVVtbsT9h3X81iToBsY/suYBJwd/qCmTUB\nJgNHE36hzzOzRwlhcF2trz8L+Iu7P2tmZcCNwCkNK11EGk379jByZDggTExbtqzmUNRJk+CGG8Lr\nZWVbnxKGDAlNR1oAL3ayCgB3f87Muta6PBhY6u7LAMxsOjDW3a8jPC3UZR2wy46XKiKxYRbWKure\nHU4+OVzbuPHLC+A99tjWr+nZs2bT0YEHagG8iJmnp5jX98YQADPSTUBmdgJQ5e7/kTo/DRji7hfU\n8fVfB74G7A7cXlcfgJmNB8YDlJWVVU6fPr3e2qqrq2nZsmVW9xFXSb8H1R+9ON5DaXU1rV5/nVav\nv07rxYtpvXgxzdatA2BL06Zs6NGDDT178lGvXqzZd19KyssT23QUp3//w4cPX+DuA+t7X6MFwM4Y\nOHCgz58/v973zZ49m2GZsyQTKOn3oPqjl4h7cIcVK2p2MC9YAJ9+Gl5v375m09GgQeFaAsTp37+Z\nZRUADem6Xw10zjjvlLrWYNoSUqRAmYVZyV261FwAb9EiltxzD/t/+GEIhyee2LoAXnl5zVA46CDY\nRa3IudCQAJgHVJhZN8Iv/hOBk3NRlLaEFCkipaVw0EG8++GH7J/+P+gNG2D+/K1PCs8+C7/7XXit\nadMQApn9CRUVWgBvJ2Q7DPQ+YBiwh5mtAn7q7r82swuAJwkjf6a5+2t5q1REikerVjB8eDjSVq+u\n2XT0m9+ErTchjDBK75mQDoUOHaKpPUGyHQV0Uh3XZwIzc1oRagISkW3o2PFLC+CxeHHNULjuuq0L\n4HXtWrPpqH9/LYBXSyyn76kJSETqlV7Qrm9fOOuscO3jj2sugDd3Ltx//9b3H3hgzVDo2bOom45i\nGQAiIjulRQs47LBwpK1ZU3MW8/Tp8MtfhtdatfryAnj77BNN7RGIZQCoCUhEcmavveDYY8MBYQG8\nN96o+ZQwcWIYjQTQqVPNzXQqK8NKqQUolgGgJiARyZuSktD007MnnH56uPbZZ7BwYc1ZzA8+uPX9\nmQvgDRkCvXsXxAJ4yb8DEZGGat4cDj44HGnvvVez6eiPf4Rf/zq8tttuYX2jzFDIclJtnCgARES2\nZY894JhjwgHhF/xbb9VsOrr11rAGEnBwu3ah7yEdCgMHhv2ZYyyWAaA+ABGJHbMwK7m8vOYCeH/7\nG7z4IuseeYS9Xn89LJWdfn/PnjX7Ew44IExki4lYBoD6AEQkEZo1C6OIBg3i9T592GvYMFi3DubN\n2/qk8PjjcNdd4f3Nm8OAATWbjrp2jWwBvFgGgIhIYrVtCyNGhANC09Hy5TX3TpgyBW6+Oby+555f\nnsXctm2jlKoAEBHJJ7Pwf/ldu8K4ceHapk2waFHN/oSZM7d2JFdUhKUuMjul8yCWAaA+ABEpaE2b\nhqUp+veHc88N1z76qOYCeHvvnfcyYhkA6gMQkaLTujV89avhaCTFuwiGiEiRUwCIiBQpBYCISJFS\nAIiIFKlYBoCZjTGzO9avXx91KSIiBSuWAeDuj7n7+DYxX0dDRCTJYhkAIiKSfwoAEZEiZR7jNazN\n7F/A8izeugfwXp7Lybek34Pqj17S70H1504Xd9+zvjfFOgCyZWbz3X1g1HU0RNLvQfVHL+n3oPob\nn5qARESKlAJARKRIFUoA3BF1ATmQ9HtQ/dFL+j2o/kZWEH0AIiKy4wrlCUBERHZQ4gPAzKrMbImZ\nLTWzy6Kupz5m1tnMZpnZ383sNTP7bup6OzP7s5m9mfpn4+wJt5PMrImZLTSzGanzpNW/u5n9wcxe\nN7PFZnZwku7BzC5O/fezyMzuM7Pmca7fzKaZ2VozW5Rxrc56zezy1M/0EjP7WjRV11THPVyf+m/o\nFTP7o5ntnvFa7O6htkQHgJk1ASYDI4HewElm1jvaquq1Gfi+u/cGhgLfSdV8GfC0u1cAT6fO4+y7\nwOKM86TVfwvwhLv3BPoR7iUR92BmHYGLgIHu3hdoApxIvOu/C6iqdW2b9aZ+Hk4E+qS+5hepn/Wo\n3cWX7+HPQF93PxB4A7gcYn0PNSQ6AIDBwFJ3X+buG4HpwNiIa9oud3/X3V9K/XkD4RdPR0Ldv0m9\n7TfAcdFUWD8z6wSMAqZmXE5S/W2Aw4FfA7j7Rnf/kATdA2E3v13NrBTYDfgHMa7f3Z8DPqh1ua56\nxwLT3f1zd38bWEr4WY/Utu7B3Z9y982p0xeATqk/x/Ieakt6AHQEVmacr0pdSwQz6wr0B+YCZe7+\nbuqlNUBZRGVl42bgUmBLxrUk1d8N+BdwZ6oZa6qZtSAh9+Duq4GJwArgXWC9uz9FQurPUFe9Sf25\nPgv4U+rPibiHpAdAYplZS+BB4Hvu/lHmax6GZsVyeJaZjQbWuvuCut4T5/pTSoEBwO3u3h/4mFrN\nJXG+h1Rb+VhCkO0DtDCzUzPfE+f6tyVp9dZmZlcQmnd/G3UtOyLpAbAa6Jxx3il1LdbMrCnhl/9v\n3f2h1OV/mtneqdf3BtZGVV89DgGONbN3CE1uXzWze0lO/RD+b2yVu89Nnf+BEAhJuYejgLfd/V/u\nvgl4CPgKyak/ra56E/VzbWZnAqOBU3zruPpE3EPSA2AeUGFm3cysGaHT5dGIa9ouMzNC2/Nid78x\n46VHgTNSfz4DeKSxa8uGu1/u7p3cvSvh3/cz7n4qCakfwN3XACvNbP/UpSOBv5Oce1gBDDWz3VL/\nPR1J6EtKSv1pddX7KHCime1iZt2ACuDFCOqrl5lVEZpDj3X3TzJeSsY9uHuiD+AYQu/7W8AVUdeT\nRb2HEh51XwFeTh3HAO0JIyHeBP4/0C7qWrO4l2HAjNSfE1U/cBAwP/X38DDQNkn3APwMeB1YBNwD\n7BLn+oH7CP0VmwhPYGdvr17gitTP9BJgZNT1b+celhLa+tM/y1PifA+1D80EFhEpUklvAhIRkZ2k\nABARKVIKABGRIqUAEBEpUgoAEZEipQAQESlSCgARkSKlABARKVL/B8fu99fEM7jIAAAAAElFTkSu\nQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "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": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHVxJREFUeJzt3XmUVNW1x/HvbkY1PBJRBsWIiiItzSAtiIC2EzSKYBAV\nJDihOIBDTIw4JRpnonEAlIegiPpAxAkVxSG2EEWkFRGQh/KMiWgUFZPYmojE8/7YmJAOQw9VfW7V\n/X3W6rWsy62q3Yu7/HHvOWcfCyEgIiLpUxC7ABERiUMBICKSUgoAEZGUUgCIiKSUAkBEJKUUACIi\nKaUAEBFJKQWAiEhKKQBERFJKASAiklL1YxewJTvssENo06ZN7DIkQb788ku222672GWIbFJSrs/X\nXnvt0xDCjls7L5EBYGZHAUe1bduW8vLy2OVIgpSVlVFSUhK7DJFNSsr1aWZ/qMp5iXwEFEJ4PIQw\nsmnTprFLERHJW4kMABERyT4FgIhISikARERSSgEgIpJSCgARkZTKzwBYvhwWLIhdhYhIouVnAFx2\nGfTsCeefDxUVsasREUmkRAaAmR1lZpP+8pe/1OwDpk2Ds8+GW2+FoiJ49tnMFigikgcSGQC1XgjW\npAmMHw/z5kHDhtCnD5x6Knz+eWYLFRHJYYkMgIzp3RuWLIGLL/a7gsJCePjh2FWJiCRCfgcAQOPG\ncO21sGgRtGoFxxwDgwfDRx/FrkxEJKr8D4DvdOkCCxfC9dfDE09A+/YwdSqEELsyEZEo0hMAAA0a\nwEUX+WOhoiI45RTo2xd+//vYlYmI1Ll0BcB32rWDsjK4/XZfL9Chg88Y+sc/YlcmIlJn0hkAAAUF\ncNZZvmjsoIN8zUCvXvDWW7ErExGpE+kNgO/88Ifw5JNw333wzjs+VnDVVbBuXezKRESySgEAYAbD\nhvm//gcNgl/8AoqLfeaQiEieUgBsrHlzmD4dHnsMPvsM9t8fLrwQvvoqdmUiIhmnANiUAQP8buD0\n0+HGG6FjR3jhhdhViYhklAJgc5o2hYkT//U//kMOgTPOgJr2JxIRSZg6DQAzO9rM7jSzB8ysT11+\nd42VlMCbb/qjoMmTvZ3E7NmxqxIRqbUqB4CZ3WVma8xsWaXjpWa20sxWmdmYLX1GCOHREMLpwJnA\n8TUrOYJtt4WxY30lcbNmMHAgDBkCa9bErkxEpMaqcwcwFSjd+ICZ1QMmAP2AQmComRWaWZGZPVHp\np/lGb71sw/tyS3ExlJf7NNFHHvF2Evfdp3YSIpKTqhwAIYR5wNpKh7sBq0II74YQ1gEzgIEhhKUh\nhP6VftaYuwF4KoTweuZ+jTrUsKFvOLN4sa8oHj4cjjwS/vjH2JWJiFRL/Vq+f2fg/Y1erwa6b+H8\nc4DDgKZm1jaEMLHyCWY2EhgJ0KJFC8rKympZYhZddRU7P/oou0+eTNh7b94dOZIPBwzwVcaSFRUV\nFcm+JiTVcu36tFCNxxdm1gZ4IoTQYcPrwUBpCOG0Da+HA91DCKMzUVxxcXEoLy/PxEdl13vvwciR\nvvNYr14+WNyuXeyq8lJZWRklJSWxyxDZpKRcn2b2WgiheGvn1fafqh8Au2z0uvWGY7VS6y0h61qb\nNjB3rreXXr4cOnXyttPffBO7MhGRzaptACwC9jSz3cysITAEqPUcyVpvCRmDGZx0ki8gO+oo34Ws\nWzcfKxARSaDqTAOdDiwA2pnZajMbEUJYD4wG5gIrgJkhhOXZKTVHtGwJDz4IDz3ku47tt5+Hwd/+\nFrsyEZF/U51ZQENDCK1CCA1CCK1DCFM2HJ8TQtgrhLBHCOGaTBSVc4+ANmXQIL8bOOkkfxzUuTPM\nnx+7KhGRf0rkdJWcfAS0KT/4AUyZ4oPD69bBgQfCqFHw17/GrkxEJJkBkHcOOwyWLfNNZ+64w3cg\nmzMndlUiknKJDIC8eARU2Xbbwc03w8svQ5Mmvnhs+HD49NPYlYlISiUyAPLmEdCm7L8/vP66bzoz\nY4Y3l3vgAbWTEJE6l8gAyHuNGsGVV3oQtGnjjeWOPho+qPUSChGRKlMAxFRUBAsWwE03+UBxYSFM\nmgTffhu7MhFJgUQGQF6OAWxOvXpwwQWwdKl3Gz3jDDj0UFi1KnZlIpLnEhkAeT0GsDl77AHPPed9\nhBYv9ruDG2+E9etjVyYieSqRAZBaZjBihC8g69vXdyHr0cN3JBMRyTAFQBLttJNvODNzpu8z0LWr\nzxr6+uvYlYlIHklkAKRqDGBzzODYY/1u4IQTfBeyffeFV16JXZmI5IlEBkAqxwA2p1kzuOceXzn8\nxRdwwAHwk5/Al1/GrkxEclwiA0A2oV8/32vg7LPhllt8kPj552NXJSI5TAGQS5o0gfHjYd48qF/f\newyddhr8+c+xKxORHKQAyEW9e8OSJTBmjO9CVlgIjz0WuyoRyTEKgFy1zTZw3XWwcCE0b+6tJI4/\nHj7+OHZlIpIjEhkAmgVUDV27wqJFcM018Oijfjdw331qLiciW5XIANAsoGpq0AAuuQTeeAPatfM2\n00ce6WsIREQ2I5EBIDXUvr1vO3nbbT5QvM8+vgGNmsuJyCYoAPJNvXpwzjm+A1mPHj5ttKQE3n47\ndmUikjAKgHzVpg3MnQt33+2dRjt1grFj1VxORP5JAZDPzODkk72dRL9+cNFF0L27TyEVkdRLZABo\nFlCGtWoFDz8Ms2b5rmPFxXD55WouJ5JyiQwAzQLKkmOO8buBYcPg6quhSxffpF5EUimRASBZtP32\nvnr46ae9oVyvXnDeeVBREbsyEaljCoC06tvXZwqNGuXTRouKfF9iEUkNBUCaNWkC48b52oFGjaBP\nHzj1VPj889iViUgdUACIPwZ64w24+GKYNs3bSTzySOyqRCTLFADiGjeGa6+FV1+Fli1h0CDfkeyj\nj2JXJiJZogCQf7fvvh4C114Ljz/udwPTpqm5nEgeUgDIf2rQwB8HvfGG9xc66SRfSPaHP8SuTEQy\nSAEgm7f33j5APG4c/O530KEDTJig5nIieSKRAaCVwAlSUACjR/uU0QMO8P8+6CBYuTJ2ZSJSS4kM\nAK0ETqA2bXzx2NSpvjl9p05w/fXwzTexKxORGkpkAEhCmfl4wFtvQf/+Pk7QvTssXhy7MhGpAQWA\nVF/Llt5YbtYs+PBD2G8/uPRS+PvfY1cmItWgAJCa+6653PDhPm20c2d46aXYVYlIFSkApHa23943\nnZk71+8AeveGc89VczmRHKAAkMzo08dnCo0eDePH+5TRZ56JXZWIbIECQDLne9/zzqLz53trib59\n4ZRTYO3a2JWJyCYoACTzevb0VcSXXAL33uvtJB56KHZVIlKJAkCyo3FjuOYaKC+HnXaCwYP9R83l\nRBJDASDZ1bkzLFzoi8aeeMLvBqZOVXM5kQRQAEj2NWgAF10ES5bAPvv4uEBpKbz3XuzKRFJNASB1\np107ePFFnyX08ss+U2jcODWXE4mkzgLAzNqb2UQzm2VmZ9XV90rCFBT4PsTLlvlOZOeeCwceCP/7\nv7ErE0mdKgWAmd1lZmvMbFml46VmttLMVpnZmC19RghhRQjhTOA4oGfNS5a8sOuu8NRTcM89vpq4\nUydfTazmciJ1pqp3AFOB0o0PmFk9YALQDygEhppZoZkVmdkTlX6ab3jPAOBJYE7GfgPJXWZw4omw\nYgUMGOD9hLp1U3M5kTpSpQAIIcwDKq/m6QasCiG8G0JYB8wABoYQloYQ+lf6WbPhc2aHEPoBwzL5\nS0iOa9ECHnzQ1wp89JE3l7v4YjWXE8my+rV4787A+xu9Xg1039zJZlYCDAIasYU7ADMbCYwEaNGi\nBWVlZbUoUXLK9ttTf9Ik9rj9dlpdfz1f3X8/Ky+8kL8UFf3zlIqKCl0Tkli5dn3WJgCqJYRQBpRV\n4bxJwCSA4uLiUFJSktW6JIGOOgqefZZtR46ky7nn+qDxdddBkyaUlZWha0KSKteuz9rMAvoA2GWj\n1603HKs1bQkpHH44LF3qs4Ruv92njM6dG7sqkbxSmwBYBOxpZruZWUNgCDA7E0VpS0gBvLncrbf6\nhvTbbgulpex93XVqLieSIVWdBjodWAC0M7PVZjYihLAeGA3MBVYAM0MIy7NXqqTWAQf4zKBLL6X5\n889D+/a+G5mI1EpVZwENDSG0CiE0CCG0DiFM2XB8TghhrxDCHiGEazJVlB4ByX9o3BiuvprXJ06E\n1q3h2GN9R7I//Sl2ZSI5K5GtIPQISDanom3bfzWXe/JJby53991qLidSA4kMAJEtql/fm8u9+SYU\nFcGpp/rmM2ouJ1ItiQwAPQKSKtlrLygrgwkTYMECnyl0223wj3/ErkwkJyQyAPQISKqsoADOPhuW\nL/cN6c87z5vLrVgRuzKRxEtkAIhU2w9/CHPmwLRp3lm0c2ffkUzN5UQ2SwEg+cMMhg/37qJHHw2X\nXeZ9hV5/PXZlIomUyADQGIDUSosW8MAD8Mgj8PHH3mF0zBj4299iVyaSKIkMAI0BSEYcfbTfDZx8\nMtxwgz8Wmj8/dlUiiZHIABDJmB/8ACZPhmefhXXrfIB41Cj44ovYlYlEpwCQdDjsMN+G8vzz4Y47\nfHP6p56KXZVIVIkMAI0BSFZstx3cfDO89JI3mjviCN+R7LPPYlcmEkUiA0BjAJJVPXp4c7nLL4fp\n072dxIMPqp2EpE4iA0Ak6xo1gl/9CsrLYZdd4LjjYNAgNZeTVFEASLp16gSvvAJjx8LTT3ur6bvu\n0t2ApIICQKR+fbjwQliyxANhxAjo0wd+//vYlYlkVSIDQIPAEsVee8ELL/gsoYUL1VxO8l4iA0CD\nwBJNQQGceaY3lzvoIG8u17u3mstJXkpkAIhEt8suvuHMfffB22+ruZzkJQWAyOaYwbBh3k7iRz/y\n5nLFxfDaa7ErE8kIBYDI1jRvDjNmwKOPwiefQPfuviOZmstJjlMAiFTVwIH/ai43dqzPGJo3L3ZV\nIjWmABCpju9/35vLPfccrF/vA8WjRsFf/xq7MpFqS2QAaBqoJN6hh8LSpfCTn/i00Q4d1FxOck4i\nA0DTQCUnbLcd/OY38PLL0KSJmstJzklkAIjklP33920nf/ELby7Xvj3MnKl2EpJ4CgCRTGjUCK68\n0qeI7rorHH+8N5f78MPYlYlslgJAJJM6doQFC+DXv/bmcoWFMGWK7gYkkRQAIplWvz787Gc+SNy5\nM5x2Ghx+OLz7buzKRP6NAkAkW9q2hd/+FiZOhFdfhaIiuOUWNZeTxFAAiGRTQQGccYYvIDv4YJ82\n2quXvxaJTAEgUhdat4bHH4f774d33oEuXeCqq2DdutiVSYopAETqihmccIK3lh40yKeN7refb0sp\nEkEiA0ArgSWv7bijrxd47DH49FNvLvfzn6u5nNS5RAaAVgJLKgwY4GMBI0b4tNGOHeHFF2NXJSmS\nyAAQSY2mTWHSJHj+efj2WygpgbPOUnM5qRMKAJEkOOQQXzdwwQUeCPvsA3PmxK5K8pwCQCQptt0W\nbrrJm8s1bQpHHgk//rGPE4hkgQJAJGm6d/fmcr/8pTeVKyyEBx5QOwnJOAWASBI1bAhXXOHN5dq0\ngSFD4Oij1VxOMkoBIJJkRUX+SOjGG+GZZ/xuYPJk3Q1IRigARJKufn346U//1Vzu9NPhsMPUXE5q\nTQEgkiu+ay733/8Nixb5NpQ336zmclJjCgCRXFJQACNH+gKyQw7xaaM9e8Ly5bErkxykABDJRd81\nl/uf/4H/+z9vLverX6m5nFRLnQaAmW1nZuVm1r8uv1ckL5nB0KF+NzB4sE8bLS72x0MiVVClADCz\nu8xsjZktq3S81MxWmtkqMxtThY+6CJhZk0JFZDN23NHvBGbPhrVrfZP6Cy+Er76KXZkkXFXvAKYC\npRsfMLN6wASgH1AIDDWzQjMrMrMnKv00N7PDgbeANRmsX0S+c9RRPhZw2mk+bbRTJygri12VJFiV\nAiCEMA9YW+lwN2BVCOHdEMI6YAYwMISwNITQv9LPGqAE2B84ATjdzDT+IJJpTZv6LKHf/tbXChx8\nMJx5Jqi1umxC/Vq8d2fg/Y1erwa6b+7kEMKlAGZ2MvBpCOHbTZ1nZiOBkQAtWrSgTP+CkY1UVFTo\nmqgKMwrGj2e3u++m9Z13su6hh1h5wQWs7dEjdmV5Ldeuz9oEQI2EEKZu5c8nAZMAiouLQ0lJSR1U\nJbmirKwMXRPVUFoKr75KoxEj6HjJJb4j2S23+LiBZFyuXZ+1eQzzAbDLRq9bbzgmIknSrZv3FLri\nCnjwQW8nMWOG2klIrQJgEbCnme1mZg2BIcDsTBSlLSFFMqxhQ58m+vrrsPvuPn104ED4QP9mS7Oq\nTgOdDiwA2pnZajMbEUJYD4wG5gIrgJkhhIwsR9SWkCJZ0qGDN5e76SZ47jm/G7jzTt0NpFRVZwEN\nDSG0CiE0CCG0DiFM2XB8TghhrxDCHiGEa7JbqohkRL163kJi6VLo2tVbSxx6qK8ollRJ5FRMPQIS\nqQN77OF7EU+a5GMERUXwm9+ouVyKJDIA9AhIpI6YeXvpt97yFtM//SkccAAsW7b190rOS2QAiEgd\n23lneOwxmD7d9xnYd1+48ko1l8tziQwAPQISicDMt55csQKOPdanjXbtCq++GrsyyZJEBoAeAYlE\ntMMOcP/93m7688+hRw/42c/UXC4PJTIARCQB+vf35nKnn+7TRouK4IUXYlclGZTIANAjIJGEaNoU\nJk70//Gb+S5kZ5yh5nJ5IpEBoEdAIglTUgJvvumPgiZP9gVkjz8euyqppUQGgIgk0Lbbwq9/Da+8\nAs2awYAB3lzuk09iVyY1pAAQkerZbz8oL/c9iGfNgvbtfUcytZPIOQoAEam+hg3h8sth8WJo2xaG\nDfM7gtWrY1cm1ZDIANAgsEiO2GcfeOkluPlm34WssNB3JPt2k/s9ScIkMgA0CCySQ+rVg/PP9+Zy\n3br5FpSHHgqrVsWuTLYikQEgIjlo993h2Wd9ltDixb5u4MYbYf362JXJZigARCRzzGDECG8u17cv\nXHihN5dbujR2ZbIJCgARybyddoJHHoEHHoD33vPmcr/8JXz9dezKZCOJDAANAovkATM47jhvLjdk\niE8b7doVFi6MXZlskMgA0CCwSB5p1gzuvReefNJbSPTo4TuSffll7MpSL5EBICJ56IgjvLncmWf6\ntNGOHX3qqESjABCRuvNf/wW33w4vvujTRw891LuN/vnPsStLJQWAiNS9Aw+EJUvg5z+Hu+7yBWWz\nZ8euKnUUACISxzbbwA03+KDwDjvAwIE+WLxmTezKUkMBICJxFRd7c7mrrvKpo4WFviOZmstlXSID\nQNNARVKmQQO47DJfQbznnvDjH/uOZO+/H7uyvJbIANA0UJGUKiyE3/0ObrkFysp8bOCOO9RcLksS\nGQAikmL16sF558GyZdC9O5x9Nhx8MLzzTuzK8o4CQESSabfd4JlnYMoUnzHUsSOMHavmchmkABCR\n5DKDU0/15nKlpXDRRbD//h4IUmsKABFJvp12gocfhpkzfWC4uNh3JFNzuVpRAIhIbjCDY4/1u4ET\nToCrr4YuXWDBgtiV5SwFgIjklmbN4J57YM4cqKiAnj19RzI1l6s2BYCI5KZ+/by53Nlnw623QocO\n8NxzsavKKQoAEcldTZrA+PEwbx40bAiHH+47kqm5XJUkMgC0ElhEqqV3b3jjDRgzxh8PFRbCo4/G\nrirxEhkAWgksItW2zTZw3XXeXK55c/jRj3xHso8/jl1ZYiUyAEREaqxrV1i0CK65Bh57zO8G7r1X\nzeU2QQEgIvmnQQO45BJ/LNSuHZx4Ihx5JPzxj7ErSxQFgIjkr/btYf58uO02HyjeZx/fkUzN5QAF\ngIjku3r14JxzvLlcjx4wahSUlMDbb8euLDoFgIikQ5s2MHcu3H03LF3qzeVuuCHVzeUUACKSHmZw\n8sneTuKII3zaaPfuqW0upwAQkfRp1cqby82aBR984M3lLrsM/v732JXVKQWAiKTXMcf43cCwYT5t\ntEsXePnl2FXVGQWAiKTb9tvD1Knw9NPw1VfQq5fvSFZREbuyrFMAiIgA9O3rM4VGjfJpo0VF8Oyz\nsavKqjoLADMrMbP5ZjbRzErq6ntFRKqsSRMYN87XDjRqBH36+I5kn38eu7KsqFIAmNldZrbGzJZV\nOl5qZivNbJWZjdnKxwSgAmgMrK5ZuSIidaBXL19FfPHFMG2at5N45JHYVWVcVe8ApgKlGx8ws3rA\nBKAfUAgMNbNCMysysycq/TQH5ocQ+gEXAVdm7lcQEcmCxo3h2mvh1VehZUsYNMh3JPvoo9iVZUyV\nAiCEMA9YW+lwN2BVCOHdEMI6YAYwMISwNITQv9LPmhDCd2uvPwcaZew3EBHJpn339RC49lp4/HG/\nG5g2LS+ay9WvxXt3Bt7f6PVqoPvmTjazQUBf4PvA+C2cNxIYCdCiRQvKyspqUaLkm4qKCl0TEkeP\nHmw7aRLtxo6l6UknsXb8eFZecAFft2z5z1Ny7fq0UMUUM7M2wBMhhA4bXg8GSkMIp214PRzoHkIY\nnaniiouLQ3l5eaY+TvJAWVkZJSUlscuQNPv2W28oN2aMryy+/no46ywoKEjM9Wlmr4UQird2Xm1m\nAX0A7LLR69YbjomI5K+CAhg92qeMHnCA//eBB8LKlbErq7baBMAiYE8z283MGgJDgNmZKEpbQopI\n4rVp44vHpk711cSdOvHD+++Hb76JXVmVVXUa6HRgAdDOzFab2YgQwnpgNDAXWAHMDCEsz0RR2hJS\nRHKCGZx0kgdA//7sPnmyN5dbvDh2ZVVS1VlAQ0MIrUIIDUIIrUMIUzYcnxNC2CuEsEcI4ZpMFaU7\nABHJKS1bwqxZLLviCvjwQ9hvP9+RLOHN5RLZCkJ3ACKSiz496CC/Gxg+3Deo79wZXnopdlmblcgA\nEBHJWdtv75vOzJ3rdwC9e/uOZF98Ebuy/6AAEBHJhj59fKbQ6NEwYQJ06OChkCCJDACNAYhIXvje\n97yz6Pz5sM02UFrqO5KtrdxYIY5EBoDGAEQkr/Ts6c3lLrkE7rvP20k89FDsqpIZACIieadxY991\nrLwcdtoJBg/2Hcn+9KdoJSUyAPQISETyVufOsHChzxJ68km/G5g6NUpzuUQGgB4BiUhea9DAewkt\nWeKDw6ec4juSvfdenZaRyAAQEUmFdu3gxRdh/HhYsMDDYNw4bzhXBxQAIiIxFRT4PsTLlvlOZOee\n62sH3n47+1+d9W8QEZGt23VXeOopuOceePddWL8+61+ZyADQILCIpJIZnHiijwUUFmb96xIZABoE\nFpFUa1Q3u+YmMgBERCT7FAAiIimlABARSalEBoAGgUVEsi+RAaBBYBGR7EtkAIiISPYpAEREUspC\nhA50VWVmnwB/qOHbmwIxBhGy8b21/cyavr8676vquVs7b2t/vgPwaRVrSroY12gSr8+afkaM63Nr\n5yTl+tw1hLDjVs8KIeTlDzApX763tp9Z0/dX531VPXdr51Xhz8tj/L1m4yfGNZrE67OmnxHj+tza\nObl2febzI6DH8+h7a/uZNX1/dd5X1XO3dl6sv7cYYvyuSbw+a/oZMa7P6n5voiX6EZBIZWZWHkIo\njl2HyKbk2vWZz3cAkp8mxS5AZAty6vrUHYCISErpDkBEJKUUACIiKaUAEBFJKQWA5DQz293MppjZ\nrNi1iFRmZkeb2Z1m9oCZ9YldT2UKAEkcM7vLzNaY2bJKx0vNbKWZrTKzMQAhhHdDCCPiVCppVM3r\n89EQwunAmcDxMerdEgWAJNFUoHTjA2ZWD5gA9AMKgaFmlv1NU0X+01Sqf31etuHPE0UBIIkTQpgH\nrK10uBuwasO/+NcBM4CBdV6cpF51rk9zNwBPhRBer+tat0YBILliZ+D9jV6vBnY2s2ZmNhHoYmYX\nxylNZNPXJ3AOcBgw2MzOjFHYltSPXYBIbYQQPsOfr4okTgjhNuC22HVsju4AJFd8AOyy0evWG46J\nJEFOXp8KAMkVi4A9zWw3M2sIDAFmR65J5Ds5eX0qACRxzGw6sABoZ2arzWxECGE9MBqYC6wAZoYQ\nlsesU9Ipn65PNYMTEUkp3QGIiKSUAkBEJKUUACIiKaUAEBFJKQWAiEhKKQBERFJKASAiklIKABGR\nlFIAiIik1P8DhUMjNffr+8cAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "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": { "collapsed": true }, "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.6.1" }, "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": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": { "height": "12px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }