{ "cells": [ { "cell_type": "markdown", "id": "ae7e2851", "metadata": {}, "source": [ "# Integration methods\n", "- Börge Göbel" ] }, { "cell_type": "markdown", "id": "efe2a841", "metadata": {}, "source": [ "![Integral](figure_05_integral.png)" ] }, { "cell_type": "code", "execution_count": 1, "id": "a527df57", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "988683ad", "metadata": {}, "source": [ "### Define and plot function" ] }, { "cell_type": "code", "execution_count": 2, "id": "f4855ad1", "metadata": {}, "outputs": [], "source": [ "def func(x):\n", " return 0.5 + 0.1 * x + 0.2 * x**2 + 0.03 * x**3" ] }, { "cell_type": "code", "execution_count": 3, "id": "d29b311a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeRElEQVR4nO3deXzV1YH38c/JThYSyL5B2EHCHhGUutUFFLFq9bFarbXqdJu2Txd9Zjozfbp3ptPO9nQ6blRabSuKiopatS6ALJoECFtYQ8i+Qvb13jN/JPhQhpAEcvP73Xu/79crryTmcu/3FeGbk3PO7/yMtRYREXGvEKcDiIjIuamoRURcTkUtIuJyKmoREZdTUYuIuFyYL540KSnJ5uTk+OKpRUQCUkFBQb21NvlsX/NJUefk5JCfn++LpxYRCUjGmNKBvqapDxERl1NRi4i4nIpaRMTlVNQiIi6nohYRcTkVtYiIy6moRURcTkUtIjIC3imuYfXmEno83hF/bhW1iMgIeHJzCWu2HiMsxIz4c6uoRUQuUG1zJ1uONHDzvAyMUVGLiLjOy7sqsRZWzc/0yfOrqEVELtDLuyrJzRzL1JRYnzy/ilpE5AIcrWulqLyJT/loNA0qahGRC7J+ZyXGwMq5GT57DRW1iMh5stayfmcFSycnkhYf5bPXUVGLiJynovImjjW0c/N8342mQUUtInLeXtpZQURoCMtz0336OipqEZHz4PFaXtlVxVUzk4kfE+7T11JRi4ichy1H6qlv7eJmH+72OEVFLSJyHtbvrCQuMoyrZ6b4/LVU1CIiw9TZ4+GNPdVcn5tGVHioz19PRS0iMkxv7quhtauXWxf4ftoDVNQiIsO2rqCcjPgolkxOHJXXU1GLiAxDbXMnmw7VccvCTEJ8cKTp2aioRUSGYf3OSrwWblmQNWqvqaIWERmGdYXlzMtO8NlJeWejohYRGaJ9lc0UV7dw28LRWUQ8RUUtIjJE6wrLCQ813OTDk/LORkUtIjIEvR4v63dWcPXMFMbFRIzqa6uoRUSGYNOheupbu7l14egtIp6iohYRGYJ1heUkRIdz1QzfXzJ+JhW1iMggmjp6eHNfDavmZRARNvq1qaIWERnEhqIqunu9jkx7wDCK2hgTaozZYYx51ZeBRETcZm1+GdNTY5mXFe/I6w9nRP11YL+vgoiIuNHBmhZ2lp3kjrxsjBmdS8bPNKSiNsZkATcCT/g2joiIu6z9qIywEMMto3RS3tkMdUT9r8DDgNd3UURE3KW718uLOyq4ZlYqibGRjuUYtKiNMSuBWmttwSCPe8gYk2+Mya+rqxuxgCIiTnmnuJaGtm7uuNiZRcRThjKivgxYZYw5BvwRuNoY8/SZD7LWPmatzbPW5iUnJ49wTBGR0fdcfhkpcZFcPs3ZThu0qK21f2OtzbLW5gB3Au9Yaz/r82QiIg6qae7k3QO13LYoi7BQZ3cyax+1iMhZvFBYgdfCHXnZTkchbDgPtta+B7znkyQiIi5hreW5/DIW54xnUlKM03E0ohYROVN+6QmO1rdxe56zi4inqKhFRM7w7EdlxESEcsOcdKejACpqEZG/0NTRw6tFlayan0FM5LBmh31GRS0icpr1Oyvo7PFy1+KJTkf5mIpaRKSftZbfbz9ObuZY5jh0ANPZqKhFRPrtKDtJcXWLq0bToKIWEfnY77cfJyYilFXzR/fmtYNRUYuIcPoiYiaxLllEPEVFLSLC6YuIE5yO8j+oqEUk6J1aRJyTGe+qRcRTVNQiEvQ+XkS8xH2jaVBRi4h8vIh40zx3LSKeoqIWkaDW1O7eRcRTVNQiEtSeKyijs8fLPUvctXf6dCpqEQlaXq/l6W2l5E0cx0UZY52OMyAVtYgErU2H6znW0M49S907mgYVtYgEsd9tPUZSbCQrct1xnOlAVNQiEpTKGtv5c3Etn1mcTUSYu6vQ3elERHzk6e2lhBjj2r3Tp1NRi0jQ6ezxsPajMq6dlUp6/Bin4wxKRS0iQefVoipOtPdwr8sXEU9RUYtI0Pnd1mNMTYll6ZREp6MMiYpaRILKrrKT7Cpv4p4lEzHGOB1nSFTUIhJUntpyjNjIMG5dmOl0lCFTUYtI0Khp7uTVokpuz8siLirc6ThDpqIWkaDx9LZSer2W+y7NcTrKsKioRSQodPZ4eGb7ca6ZlcrExBin4wyLilpEgsL6nRU0tnXz+ctynI4ybCpqEQl41lpWbz7GzLQ4lk72jy15p1NRi0jA23KkgQM1Ldy/bJLfbMk7nYpaRALe6s0lJMZEsMqlt9oajIpaRAJaSX0b7xyo5e4lE4kKD3U6znlRUYtIQHvqgxLCQgyfXeL+U/IGoqIWkYB1oq2btfnlrJqXSUpclNNxzpuKWkQC1tPbSuno8fDQ5ZOdjnJBVNQiEpA6ezys2XqMK2ckMyMtzuk4F0RFLSIB6YXCCupbu/1+NA0qahEJQF6v5YlNR5mTGe+XF7icSUUtIgHnrf01HK1v46HLJ/vlBS5nUlGLSMB5bONRssaNYUVumtNRRsSgRW2MiTLGfGiM2WWM2WuM+f5oBBMROR8FpY0UlJ7ggWWTCAsNjLFo2BAe0wVcba1tNcaEA5uNMa9ba7f5OJuIyLA9+v5REqLDuePibKejjJhBf9zYPq39n4b3v1mfphIROQ+Ha1t4a38N9y6ZSHTEUMah/mFIvxcYY0KNMTuBWuAta+32szzmIWNMvjEmv66uboRjiogM7j/fO0JUWCj3XTbJ6SgjakhFba31WGvnA1nAYmNM7lke85i1Ns9am5ecnDzCMUVEzq2ssZ31Oyu565IJjI+JcDrOiBrWTLu19iTwHrDcF2FERM7Xf71/hFBjePAT/n+By5mGsusj2RiT0P/xGOAaoNjHuUREhqymuZPn8su5bVEWafH+e/jSQIYy254OrDHGhNJX7Gutta/6NpaIyNA9sekovV4vX7piitNRfGLQorbWFgELRiGLiMiwnWjr5pntx1k1L4MJidFOx/GJwNgNLiJB6zcflNDe7eHLV011OorPqKhFxG+1dPbw1JZjXD87lemp/n2U6bmoqEXEb/12aynNnb18+crAHU2DilpE/FRLZw+PbzrKVTOSmZed4HQcn1JRi4hfWrPlGCfbe/jGNdOdjuJzKmoR8TvNnT08vqmET85MCfjRNKioRcQPPfXBMZo6gmM0DSpqEfEzzZ09PLHpKNfMSmVOVrzTcUaFilpE/MpvNh+jubOXb1wzzekoo0ZFLSJ+o6mjhyc2H+Xai1LJzQyO0TSoqEXEj6zeXEJLkI2mQUUtIn6isa2bJzeXsHx2GrMzgmc0DSpqEfETv37vMO3dvXzruuDY6XE6FbWIuF7lyQ7WbC3l1oVZTAvgMz0GoqIWEdf79z8fAkvQzU2foqIWEVc7UtfK2vwy7l4ygaxxgXne9GBU1CLiar988yBR4aF8JYDPmx6MilpEXGt3eRMbdlfxwLJJJMVGOh3HMSpqEXGtf/pTMeOiw3ng8sC7s/hwqKhFxJU2H6pn06F6vnzlVMZGhTsdx1EqahFxHY/X8qMN+8gaN4Z7lk50Oo7jVNQi4jrrCssprm7hkeUziQoPdTqO41TUIuIq7d29/OLNA8zPTmDl3HSn47iCilpEXOXxjSXUNHfxdzfOwhjjdBxXUFGLiGvUNnfy6MYjrMhNIy9nvNNxXENFLSKu8S9vH6TH4+WR5TOdjuIqKmoRcYUD1S08+1EZ9yzJIScpxuk4rqKiFhHHWWv5/it7iYsK56+vDt5LxQeiohYRx72xp5otRxr41nXTGRcT4XQc11FRi4ijOro9/GjDfmamxXHX4glOx3GlMKcDiEhwe3TjESpOdvDHh5YQFqqx49nouyIijilrbOfX7x1h5dx0lkxOdDqOa6moRcQxP3ltP8bA394wy+korqaiFhFHfHC4ntf3VPOVK6eSkTDG6TiupqIWkVHX3evley/vJXv8GB4M8rOmh0KLiSIy6h7fdJTDta2svi9Pp+MNgUbUIjKqjje08+9/PsTy2WlcPTPV6Th+QUUtIqPGWsvfr99DWIjhe6sucjqO31BRi8ioeW13Ne8frOOb180gPV4LiEM1aFEbY7KNMe8aY/YbY/YaY74+GsFEJLC0dPbw/Vf2MjtjLJ/T7bWGZSiLib3At6y1hcaYOKDAGPOWtXafj7OJSAD5xZsHqWvt4rF783QF4jAN+t2y1lZZawv7P24B9gOZvg4mIoFjx/ET/HbrMe5ZMpH52QlOx/E7w/qxZozJARYA28/ytYeMMfnGmPy6uroRiici/q6r18PDzxeROjaKb18/w+k4fmnIRW2MiQXWAd+w1jaf+XVr7WPW2jxrbV5ycvJIZhQRP/ardw5zqLaVn9wyh7FR4U7H8UtDKmpjTDh9Jf2MtfYF30YSkUCxt7KJ/3zvCLcuyOSqmSlOx/FbQ9n1YYAngf3W2l/6PpKIBIIej5eHny8iITqCf7hJe6YvxFBG1JcB9wBXG2N29r/d4ONcIuLnHtt4lL2VzfzoU7NJiNZdWy7EoNvzrLWbATMKWUQkQByubeHf3j7EjXPSWZ6b7nQcv6fNjCIyono8Xr65dhcxkaH831WznY4TEHR6noiMqP945zBF5U3812cXkhwX6XScgKARtYiMmMLjJ/jVu4e5dWGmpjxGkIpaREZEe3cv33x2J2ljozTlMcI09SEiI+LHG/ZT2tjOHx5cogtbRphG1CJywd4truWZ7cd58BOTdTdxH1BRi8gFqW3p5DvP72JGahzfvHa603ECkqY+ROS8eb2Wbz67i5bOXp55YInuf+gjKmoROW+/fv8Imw/X89Nb5zAjLc7pOAFLUx8icl7yjzXyy7cOsnJuOndenO10nICmohaRYTvZ3s3X/rCDzIQx/PTWOfSd3Sa+oqkPERkWay0PP19EXWsX6750KXHaiudzGlGLyLA8samEN/fV8MjymczNSnA6TlBQUYvIkG090sDP3ihmRW4aX1g2yek4QUNFLSJDUtXUwVd/X0hOYjQ/v32e5qVHkYpaRAbV1evhi08X0tnj4dF78oiN1PLWaNJ3W0QG9f1X9rGr7CS/vnshU1NinY4TdDSiFpFzevaj4/x++3H+6orJrJijo0udoKIWkQFtO9rA3720h2VTk/jOdTOcjhO0VNQiclalDW186ekCssdH86u7FxIWqrpwir7zIvI/NHf28IU1+Vhg9ecuJn6MLmpxkopaRP5Cr8fLV54p5Fh9G7++exE5STFORwp62vUhIh+z1vLDV/ex6VA9/3jbHJZO0U0A3EAjahH52KMbj7JmaykPfmIS/+viCU7HkX4qahEB4IXCcn72ejE3zcvgb1bMcjqOnEZFLSK8f7COh58v4tIpifzz7XMJCdHl4W7iqjnqrzxTSHp8FHk541g0cTzJcZFORxIJeEXlJ/nS0wVMS43j0XsWERmm22m5jWuKuqvXQ11rF2/tr+GJzSUA5CRGk5cznqWTE1k6JZGMhDEOpxQJLCX1bdz/1EeMj4lgzecv1tnSLuWaoo4MC2XtXy2lq9fDnopmCkob+ejYCd7eX8PzBeUATEyM5tIpiSybmsyyqUnER+svlcj5Kmts567Ht2EtrLl/MSljo5yOJAMw1toRf9K8vDybn58/Is/l9VqKq1vYerSBrUca2F7SQEtnLyEG5mUncPm0ZK6ckcy8rATNq4kMUVVTB3c8upXmjl7++NASZqWPdTpS0DPGFFhr8876NbcX9Zl6PV52lZ/k/YP1bDxYR1H5SbwWkmIjuHJGCp+cmcKyaUn6FU5kALUtndz56DbqWrp45sFLdJcWlwiooj7TyfZu3j9YxzvFtbx3oI6mjh4iQkO4dGoi112UxjUXpZASp1/pRAAa27q587GtlJ/o4HdfWMyiieOdjiT9ArqoT9fr8VJQ2jev/ae9NRxvbMcYWJCdwIrcdJbnppE9PnrUc4m4QV1LF/c8uZ2S+jZ+8/mLuXRKktOR5DRBU9Sns9ZysKaVN/dW88beavZWNgMwNyueFbnp3DgnnQmJKm0JDtVNndz1xDaqTnbyxOfyuGyqStptgrKoz3S8oZ3X91Tx2p5qdpWdBPoWI2+am86Nc9NJj9fWPwlM5Sfauevx7TS2dbP6votZPEnTHW6koj5D+Yl2NhRV8UpRJXsq+kbai3PGs2p+BjfOSWdcTITDCUVGxrH6Nu56fButXb2suX8xCyaMczqSDEBFfQ4l9W28squSl3dVcri2lbAQwxXTk7l5QSbXzkplTISu0hL/tLeyift+8xG9Hi+/+8Il5GbGOx1JzkFFPQTWWvZWNrN+ZwUv76qkprmL2MgwVuSmccvCTJZMStQ+bfEbmw/V88WnCxgbFcaa+xczLTXO6UgyCBX1MHm8lu1HG3hxRwWv7a6irdtDRnwUn1qQyW2LspiSrLswi3u9tKOCbz+3i6kpsTz1+cWkxWt7qj+4oKI2xqwGVgK11trcobygvxf16Tq6Pby1v4YXCsvZeLAOr4X52QnctiiLm+amkxCt+WxxB2stj208yk9fL2bJ5PE8dm8eY3Xhl9+40KK+HGgFfhuMRX262pZO1u+oZF1hOcXVLUSEhnDt7FQ+vSiLy6clE6qpEXFId6+Xf1i/hz9+VMbKuen84o55OgXPz1zw1IcxJgd4NdiL+pRT89nPF5SzfmcFJ9p7SB0byS0Lsrg9T1MjMroaWrv40tOFfHiskS9fOYVvXzdD6yl+aFSK2hjzEPAQwIQJExaVlpaeX1o/093r5Z3ivhP+3j1Qh8drWTghgdvzslk5N11njohP7ats5sHf5lPf2sU/fXouN8/PdDqSnCeNqEdJbUsnL+2oYG1+OYdrW4kKD+GG3HRuz8vmkknjNcqREbWhqIpvP7eL+DHhPHbvIh2u5OfOVdSuOY86EKTERfHQ5VN48BOT2Vl2krX55by6q5IXdlQwYXw0n16UxW2LssjUDRDkAnT1evjJhv2s2VrKggkJPPrZRTpLOsBpRO1jHd0e3thbxXP55Ww50oAxsGxqEnfkZXPtRalEhWvBR4bueEM7X/1DIUXlTTywbBIPL59JRJhufRoILnTXxx+AK4EkoAb4nrX2yXP9GRX12ZU1tvNcQTnrCsqpONlB/Jhwbp6fwe2LssnNHIsxmhqRgb2xp4rvPF+EAf759nlcNzvN6UgygnTBi8t4vZYtRxpYm1/GG3ur6e71MjMtjk8vyuKWBZkkxuqmvvL/NXf28INX9vF8QTnzsuL5f3ct1HG9AUhF7WJN7T28XFTJ8/ll7CpvIizEcPXMFD69KIurZqYQHqpfa4PZliP1fOe5IqqaOvjylVP52ienaaojQKmo/cTBmhaeyy/jxR2V1Ld2kRgTwc3zM7ltUSazM3SgTjDp6Pbw8z8dYPUHJUxKiuEXd8xjoU6+C2gqaj/T4/Gy8WAd6wrLeXtfLd2evqmR2xZmcfP8DK3wB7h3D9Ty9y/tofxEB59bOpFHVswkOkIbtAKditqPnWjr5pWiSl4orGBn2UlCDFw+PZlbFmRy3UVpOoY1gNQ0d/KDV/axYXcVU5Jj+PEtc1gyOdHpWDJKVNQB4nBtKy/uKOfFwgoqmzqJiQjl+tw0blmQyaVTknTWiJ/q7vXy9LZS/uWtg3R5vPz1VVN56IrJOqsjyKioA4zXa/nwWCMv7ahgw+4qWjp7SY6LZOXcdFbNy2B+doK2+vkBay1v7qvhZ68XU1LfxiemJfHDm3PJSYpxOpo4QEUdwDp7PLxTXMv6nRW8W1xHt8fLhPHRrJqXwcp56cxIjVNpu9Du8iZ+tGEf20samZoSy3dvmMWVM5L1/yqIqaiDRHNnD3/aU83Luyr54HA9XgtTkmO4cW4GK+emM113+XDc3som/u3tQ7y5r4bxMRH872un85mLswnTNsygp6IOQvWtXbyxp5pXiyrZXtKI7S/tFbnpLM9NY3aGroQcTacXdFxUGA8sm8znl+XoYH/5mIo6yNW2dPLGnmpe313N9pIGvBayx4/h+ovSuPaiVBZNHKcRnQ9Ya9l4qJ4nN5ew8WDdxwV932U5xI9RQctfUlHLxxpau3h7fw2v76lmy+EGuj1exkWHc9XMFK6ZlcqyaUka5V2gzh4PLxRWsPqDEg7XtpISF8m9Sydyz1IVtAxMRS1n1drVy8aDdby9r4Z3DtRysr2HsBDDoonjuGpmClfNSGF6aqymSIbAWsueimaezT/O+p2VtHT2MjtjLF9YNomVczN02bcMSkUtg+r1eCkoPcF7B+t470Ad+6uaAUgdG8llU5P4xLQkLpuSpKsiz1Dd1Mlru6tYm19GcXULkWEh3DAnnTsvzmbxpPH6ISdDpqKWYatu6uT9g7VsPFTPlsP1nGjvAWBaSiyXTB7PJZMSuWTS+KAs7qqmDl7fXc1ru6vILz0BwNyseG7Py2bVvAxNb8h5UVHLBfF6Lfuqmtl0qJ6tRxsoONZIW7cHgElJMSycMI4FExJYOGEc01NjA25hsru377eNjYfqeP9AHfv6f9uYmRbHjXPSWTEnnakpuqGxXBgVtYyoXo+XvZXNfFjSyPaSRnYcP0FDWzcA0RGh5GbEMztz7Mfvpyb7V3m3d/eyq6yJgtJGCkpP8GFJ3w+msBDDwonjuGJ6Mity05isu83LCFJRi09Zaylr7GBH2QkKS0+wu6KJ/VUtdPT0jbojQkOYnBzD1JRYpqfGMS0llomJMUxIjCY20tlT4RrbuimuamZ/dUv/+2b2V7Xg8fb9u5iWEsviSeO5YnoyS6ck6q7y4jMqahl1Hq+lpL6VvZXN7Ktq5nBNKwdrWyhr7PiLxyXGRDAhMZqMhDGkxkWROjaStPgokmMjGTsmnITocBKiI4iJCB3Wwlyvx0tbl4eGti7qW7tpaO2ivrWLipOdlDW2c7z/ramj5+M/kxQbycy0OOZnJ7Bo4jgWThhHfLSKWUaHilpco727l6N1bZQ2nCrLvo+rmzqpbu6kvX/u+0yhIYYx4aFEhoUQGRZCVHgoxoAFrO0b1fd4LB09Hlq7eunu9Z71ecJDDdnjoskeH82E8dFMTIxmRlocM9PGkhynW6CJc85V1DqNXEZVdEQYuZnx5Gae/Y41LZ091DR3UdfSRVNHD00d3f3ve+js8dLZ46Grt++9tWAMGGMIMX1lHhMRRnRkaN/7iFASYyNIio0kKTaSxNgIEmMidRys+B0VtbhKXFQ4cVHh2kUhchr/WYoXEQlSKmoREZdTUYuIuJyKWkTE5VTUIiIup6IWEXE5FbWIiMupqEVEXM4nl5AbY+qA0hF/4guTBNQ7HWKIlNV3/CmvP2UF/8rrxqwTrbXJZ/uCT4rajYwx+QNdR+82yuo7/pTXn7KCf+X1p6ygqQ8REddTUYuIuFwwFfVjTgcYBmX1HX/K609Zwb/y+lPW4JmjFhHxV8E0ohYR8UsqahERlwuaojbG/NAYU2SM2WmMedMYk+F0pnMxxvzcGFPcn/lFY0yC05kGYoy53Riz1xjjNca4csuTMWa5MeaAMeawMeb/OJ3nXIwxq40xtcaYPU5nGYwxJtsY864xZn//34GvO53pXIwxUcaYD40xu/rzft/pTEMRNHPUxpix1trm/o+/Blxkrf2iw7EGZIy5DnjHWttrjPlHAGvtIw7HOitjzCzACzwKfNta66obZhpjQoGDwLVAOfAR8Blr7T5Hgw3AGHM50Ar81lqb63SeczHGpAPp1tpCY0wcUAB8ysXfWwPEWGtbjTHhwGbg69babQ5HO6egGVGfKul+MfTdF9W1rLVvWmt7+z/dBmQ5medcrLX7rbUHnM5xDouBw9bao9babuCPwM0OZxqQtXYj0Oh0jqGw1lZZawv7P24B9gOZzqYamO3T2v9peP+bq7sAgqioAYwxPzbGlAF3A//gdJ5huB943ekQfiwTKDvt83JcXCb+yhiTAywAtjsc5ZyMMaHGmJ1ALfCWtdbVeSHAitoY87YxZs9Z3m4GsNZ+11qbDTwDfNXZtIPn7X/Md4Fe+jI7ZihZXexstx13/SjKnxhjYoF1wDfO+O3Vday1HmvtfPp+S11sjHH19BIE2F3IrbXXDPGhvwc2AN/zYZxBDZbXGPM5YCXwSevwYsIwvrduVA5kn/Z5FlDpUJaA0z/Xuw54xlr7gtN5hspae9IY8x6wHHD1wm1AjajPxRgz7bRPVwHFTmUZCmPMcuARYJW1tt3pPH7uI2CaMWaSMSYCuBN42eFMAaF/ce5JYL+19pdO5xmMMSb51A4qY8wY4Bpc3gUQXLs+1gEz6NudUAp80Vpb4WyqgRljDgORQEP/f9rm1l0qxphbgP8AkoGTwE5r7fWOhjqDMeYG4F+BUGC1tfbHziYamDHmD8CV9B3FWQN8z1r7pKOhBmCMWQZsAnbT928L4G+tta85l2pgxpi5wBr6/h6EAGuttT9wNtXggqaoRUT8VdBMfYiI+CsVtYiIy6moRURcTkUtIuJyKmoREZdTUYuIuJyKWkTE5f4bBpthcdnre3MAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x_list = np.linspace(-3.5, 3.5, 71)\n", "plt.plot(x_list, func(x_list))" ] }, { "cell_type": "markdown", "id": "68683ccc", "metadata": {}, "source": [ "### Analytical solution" ] }, { "cell_type": "markdown", "id": "e7884b29", "metadata": {}, "source": [ "\\\\( f(x) = \\frac{1}{2} + \\frac{1}{10}x + \\frac{1}{5}x^2 + \\frac{3}{100}x^3 \\\\)\n", "\n", "\\\\( A = \\int_{-3}^3 f(x)\\,\\mathrm{d}x = \\int_{-3}^3\\left(\\frac{1}{2} + \\frac{1}{10}x + \\frac{1}{5}x^2 + \\frac{3}{100}x^3\\right)\\,\\mathrm{d}x = \\left[\\frac{1}{2}x + \\frac{1}{20}x^2 + \\frac{1}{15}x^3 + \\frac{3}{400}x^4\\right]_{-3}^3\\\\)" ] }, { "cell_type": "code", "execution_count": 4, "id": "a6d17b48", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.6" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "( 1/2*(3) + 1/20*(3)**2 + 1/15*(3)**3 + 3/400*(3)**4 ) - ( 1/2*(-3) + 1/20*(-3)**2 + 1/15*(-3)**3 + 3/400*(-3)**4 )" ] }, { "cell_type": "markdown", "id": "5ec28ee8", "metadata": {}, "source": [ "### Create data" ] }, { "cell_type": "code", "execution_count": 5, "id": "fd0607bb", "metadata": {}, "outputs": [], "source": [ "x_points = np.linspace(-3,3,13)\n", "data = np.array([x_points, func(x_points)])" ] }, { "cell_type": "code", "execution_count": 6, "id": "73bae40c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-3. , -2.5 , -2. , -1.5 , -1. , -0.5 ,\n", " 0. , 0.5 , 1. , 1.5 , 2. , 2.5 ,\n", " 3. ],\n", " [ 1.19 , 1.03125, 0.86 , 0.69875, 0.57 , 0.49625,\n", " 0.5 , 0.60375, 0.83 , 1.20125, 1.74 , 2.46875,\n", " 3.41 ]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data" ] }, { "cell_type": "code", "execution_count": 7, "id": "55ae583c", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAARUklEQVR4nO3df4hdZX7H8c9nx+k6rJb5IxdMxsT0jzC0q60jQ1zZUsLWdjTImi4uKGW32D+ColRhd9jGBWW3LFIGpLgW04CyCtZFajqVbWTWUkX9I+5OfmiMcUpoETMjzawyiYPDNsl++8ecyU6ud+beO3PuPfc89/2CS84955lzv4fcfHLmOc95jiNCAIDy+0LRBQAA8kGgA0AiCHQASASBDgCJINABIBEEOgAk4rJ6DWxfLul1SV/M2v9LRDxS1WaHpH+T9D/Zqv0R8cPV9rthw4bYunVr8xUDQBc7dOjQryKiUmtb3UCX9GtJX4uIedu9kt60/XJEHKxq90ZE3NZoUVu3btXk5GSjzQEAkmx/sNK2uoEei3cezWdve7MXdyMBQIdpqA/ddo/to5JOS3olIt6q0ewm22/bftn2l/MsEgBQX0OBHhEXIuJ6SVdL2m772qomhyVdExF/JOnHksZr7cf2btuTtidnZ2fXXjUA4HOaGuUSEXOSXpN0S9X6sxExny0fkNRre0ONn98XEcMRMVyp1OzTBwCsUd1At12x3Z8t90m6WdL7VW2usu1seXu2349zrxYAsKJGRrlslPSM7R4tBvULEfEz2/dIUkTslXSHpHttn5e0IOnOYBpHALjE+JFpjU1MaWZuQZv6+zQ6MqhdQwO57d9F5e7w8HAwbBFAtxg/Mq09+49p4dyFi+v6env06DeuayrUbR+KiOFa27hTFADaYGxi6pIwl6SFcxc0NjGV22cQ6ADQBjNzC02tXwsCHQDaYFN/X1Pr14JAB4A2GB0ZVF9vzyXr+np7NDoymNtnNDLKBQCwTksXPls5yoVAB4A22TU0kGuAV6PLBQASQaADQCIIdABIBIEOAIkg0AEgEQQ6ACSCQAeARBDoAJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEhE3UC3fbntX9h+2/Zx2z+o0ca2H7d90vY7tm9oTbkAgJU08pDoX0v6WkTM2+6V9KbtlyPi4LI2t0ralr1ulPRk9icAoE3qnqHHovnsbW/2iqpmt0t6Nmt7UFK/7Y35lgoAWE1Dfei2e2wflXRa0isR8VZVkwFJHy57fypbBwBok4YCPSIuRMT1kq6WtN32tVVNXOvHqlfY3m170vbk7Oxs08UCAFbW1CiXiJiT9JqkW6o2nZK0edn7qyXN1Pj5fRExHBHDlUqluUoBAKtqZJRLxXZ/ttwn6WZJ71c1e0nSt7PRLl+RdCYiPsq7WADAyhoZ5bJR0jO2e7T4H8ALEfEz2/dIUkTslXRA0k5JJyV9JunuFtULAFhB3UCPiHckDdVYv3fZcki6L9/SAADN4E5RAEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQASQaADQCIIdABIBIEOAIkg0AEgEQQ6ACSCQAeARBDoAJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEhE3UC3vdn2q7ZP2D5u+4EabXbYPmP7aPZ6uDXlAgBWclkDbc5L+k5EHLZ9paRDtl+JiPeq2r0REbflXyIAoBF1z9Aj4qOIOJwtfyrphKSBVhcGAGhOI2foF9neKmlI0ls1Nt9k+21JM5K+GxHH118eALTe+JFpjU1MaWZuQZv6+zQ6MqhdQ+U7b2040G1fIelFSQ9GxNmqzYclXRMR87Z3ShqXtK3GPnZL2i1JW7ZsWWvNAJCb8SPT2rP/mBbOXZAkTc8taM/+Y5JUulBvaJSL7V4thvlzEbG/entEnI2I+Wz5gKRe2xtqtNsXEcMRMVypVNZZOgCs39jE1MUwX7Jw7oLGJqYKqmjtGhnlYklPSToREY+t0OaqrJ1sb8/2+3GehQJAK8zMLTS1vpM10uXyVUnfknTM9tFs3UOStkhSROyVdIeke22fl7Qg6c6IiPzLBYB8berv03SN8N7U31dANetTN9Aj4k1JrtPmCUlP5FUUALTL6MjgJX3oktTX26PRkcECq1qbpka5AEBqli58dtUoFwBI1a6hgVIGeDXmcgGARBDoAJAIAh0AEkGgA0AiCHQASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQASQaADQCIIdABIBIEOAIkg0AEgEXUD3fZm26/aPmH7uO0HarSx7cdtn7T9ju0bWlMuAGAllzXQ5ryk70TEYdtXSjpk+5WIeG9Zm1slbcteN0p6MvsTANAmdc/QI+KjiDicLX8q6YSkgapmt0t6NhYdlNRve2Pu1QIAVtRUH7rtrZKGJL1VtWlA0ofL3p/S50MfANBCDQe67SskvSjpwYg4W725xo9EjX3stj1pe3J2dra5SgEAq2oo0G33ajHMn4uI/TWanJK0edn7qyXNVDeKiH0RMRwRw5VKZS31AgBW0MgoF0t6StKJiHhshWYvSfp2NtrlK5LORMRHOdYJAKijkVEuX5X0LUnHbB/N1j0kaYskRcReSQck7ZR0UtJnku7OvVIAwKrqBnpEvKnafeTL24Sk+/IqCgDQPO4UBYBEEOgAkAgCHQASQaADQCIIdABIBIEOAIkg0AEgEQQ6ACSCQAeARBDoAJCIRuZyAYDCjB+Z1tjElGbmFrSpv0+jI4PaNcTjFmoh0AF0rPEj09qz/5gWzl2QJE3PLWjP/mOSRKjXQJcLgI41NjF1McyXLJy7oLGJqYIq6mwEOoCONTO30NT6bkegA+hYm/r7mlrf7Qh0AB1rdGRQfb09l6zr6+3R6MhgQRV1Ni6KAuhYSxc+GeXSGAIdQEfbNTRAgDeILhcASASBDgCJINABIBEEOgAkgkAHgEQQ6ACQCAIdABJRN9BtP237tO13V9i+w/YZ20ez18P5lwkAqKeRG4t+IukJSc+u0uaNiLgtl4oAAGtS9ww9Il6X9EkbagEArENefeg32X7b9su2v5zTPgEATchjLpfDkq6JiHnbOyWNS9pWq6Ht3ZJ2S9KWLVty+GgAwJJ1n6FHxNmImM+WD0jqtb1hhbb7ImI4IoYrlcp6PxoAsMy6A932VbadLW/P9vnxevcLAGhO3S4X289L2iFpg+1Tkh6R1CtJEbFX0h2S7rV9XtKCpDsjIlpWMQCgprqBHhF31dn+hBaHNQIACsSdogCQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJIJAB4BEEOgAkAgCHQASQaADQCIIdABIRB4PuADQpcaPTGtsYkozcwva1N+n0ZFB7RoaKLqsrlWqQOfLA3SO8SPT2rP/mBbOXZAkTc8taM/+Y5LEv8uClKbLZenLMz23oNBvvzzjR6aLLg3oSmMTUxfDfMnCuQsam5gqqCKUJtD58gCdZWZuoan1aL3SBDpfHqCzbOrva2o9Wq80gc6XB+gsoyOD6uvtuWRdX2+PRkcGC6oIpQl0vjxAZ9k1NKBHv3GdBvr7ZEkD/X169BvXcUG0QKUZ5bL0JWGUC9A5dg0N8G+wg5Qm0CW+PACwmtJ0uQAAVkegA0AiCHQASESp+tDbhSkGAJRR3TN020/bPm373RW22/bjtk/afsf2DfmX2T5MMQCgrBrpcvmJpFtW2X6rpG3Za7ekJ9dfVnGYYgBAWdUN9Ih4XdInqzS5XdKzseigpH7bG/MqsN2YYgBAWeVxUXRA0ofL3p/K1pUSUwwAKKs8At011kXNhvZu25O2J2dnZ3P46PwxxQCAssoj0E9J2rzs/dWSZmo1jIh9ETEcEcOVSiWHj84f81MAKKs8hi2+JOl+2z+VdKOkMxHxUQ77LQxTDAAoo7qBbvt5STskbbB9StIjknolKSL2Sjogaaekk5I+k3R3q4oFAKysbqBHxF11toek+3KrCACwJtz6DwCJ4Nb/gjC9AIC8EegFWJpeYOmO1KXpBSQR6gDWjEAvwGrTCxDoyAO/AXYnAr0ATC+AVuI3wO7FRdECML0AWokJ5roXgV4AphdAK/EbYPci0AvA9AJoJX4D7F70oReE6QXQKqMjg5f0oUv8BtgtCHQgMUsnCoxy6T4EOpAgfgPsTvShA0AiOENPHDeYAN2DQE8YN5gA3YUul4RxgwnQXQj0hHGDCdBdCPSEcYMJ0F0I9IQxxQDQXbgomjBuMAG6C4GeOG4w6TwMJUWrEOhAGzGUFK1EoGPdOONsHE+rQisR6FgXzjibw1BStBKjXLAu3LzUHIaSopUIdKwLZ5zNYSgpWqmhQLd9i+0p2ydt/22N7Ttsn7F9NHs9nH+p6ESccTaHp1Whler2odvukfSPkv5M0ilJv7T9UkS8V9X0jYi4rQU1ooO16+k47brw2o7PYSgpWqWRi6LbJZ2MiP+WJNs/lXS7pOpARxdqx81L7brwygVelF0jgT4g6cNl709JurFGu5tsvy1pRtJ3I+J4DvWhBFp9xtmuoX4MKUTZNRLorrEuqt4flnRNRMzb3ilpXNK2z+3I3i1ptyRt2bKluUrRtdp14ZULvCi7Ri6KnpK0edn7q7V4Fn5RRJyNiPls+YCkXtsbqncUEfsiYjgihiuVyjrKRjdp14VXLvCi7BoJ9F9K2mb792z/jqQ7Jb20vIHtq2w7W96e7ffjvItFd2rXUD+GFKLs6na5RMR52/dLmpDUI+npiDhu+55s+15Jd0i61/Z5SQuS7oyI6m4ZYE3aNWsks1Oi7FxU7g4PD8fk5GQhnw0AZWX7UEQM19rGnaIAkAgCHQASQaADQCIIdABIBIEOAIkobJSL7VlJH6zxxzdI+lWO5RSJY+lMqRxLKschcSxLromImndmFhbo62F7cqVhO2XDsXSmVI4lleOQOJZG0OUCAIkg0AEgEWUN9H1FF5AjjqUzpXIsqRyHxLHUVco+dADA55X1DB0AUKW0gW7772y/kz2U+ue2NxVd01rZHrP9fnY8/2q7v+ia1sr2N20ft/0b26UbkVDvgehlYftp26dtv1t0Letle7PtV22fyL5bDxRd01rYvtz2L2y/nR3HD3L/jLJ2udj+3Yg4my3/jaQ/iIh7Ci5rTWz/uaT/zKYq/ntJiojvFVzWmtj+fUm/kfRPWnwUYWmm1MweiP5fWvZAdEl31Xggesez/SeS5iU9GxHXFl3PetjeKGljRBy2faWkQ5J2le3vJXtmxJeyJ7v1SnpT0gMRcTCvzyjtGfpSmGe+pM8/Fq80IuLnEXE+e3tQi0+FKqWIOBERU0XXsUYXH4geEf8naemB6KUTEa9L+qToOvIQER9FxOFs+VNJJ7T4rONSiUXz2dve7JVrbpU20CXJ9o9sfyjpLyU9XHQ9OflrSS8XXUSXqvVA9NIFR8psb5U0JOmtgktZE9s9to9KOi3plYjI9Tg6OtBt/4ftd2u8bpekiPh+RGyW9Jyk+4utdnX1jiVr831J57V4PB2rkWMpqUYeiI6C2L5C0ouSHqz6Db00IuJCRFyvxd/Ct9vOtTus7iPoihQRNzfY9J8l/bukR1pYzrrUOxbbfyXpNkl/2umP72vi76Vs6j4QHcXI+pxflPRcROwvup71iog5269JukVSbheuO/oMfTW2ty17+3VJ7xdVy3rZvkXS9yR9PSI+K7qeLlb3gehov+xi4lOSTkTEY0XXs1a2K0sj2Gz3SbpZOedWmUe5vChpUIsjKj6QdE9ETBdb1drYPinpi5I+zlYdLPGInb+Q9GNJFUlzko5GxEihRTXB9k5J/6DfPhD9R8VWtDa2n5e0Q4uz+v2vpEci4qlCi1oj238s6Q1Jx7T4712SHoqIA8VV1TzbfyjpGS1+t74g6YWI+GGun1HWQAcAXKq0XS4AgEsR6ACQCAIdABJBoANAIgh0AEgEgQ4AiSDQASARBDoAJOL/ASLrZgknsY5nAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.scatter(data[0], data[1])" ] }, { "cell_type": "markdown", "id": "c577a8d8", "metadata": {}, "source": [ "## 1. Sum" ] }, { "cell_type": "markdown", "id": "a514cd55", "metadata": {}, "source": [ "\\\\( A = \\int_{a}^b f(x)\\,\\mathrm{d}x \\approx \\frac{b-a}{n-1}\\sum_{i=1}^n f(x_i)\\\\)\n", "\n", "Edges are problematic, as they are overrepresented. Furthermore, this only really works if the data is equidistant." ] }, { "cell_type": "code", "execution_count": 8, "id": "f516346d", "metadata": {}, "outputs": [], "source": [ "def integralSum(data):\n", " return np.sum(data[1]) * (data[0,-1] - data[0,0]) / ( len(data[1]) - 1 )" ] }, { "cell_type": "code", "execution_count": 9, "id": "85729d4b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "7.800000000000001" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integralSum(data)" ] }, { "cell_type": "markdown", "id": "7b9ca0b4", "metadata": {}, "source": [ "## 2. Trapezoidal method " ] }, { "cell_type": "markdown", "id": "d30bc808", "metadata": {}, "source": [ "Corresponds to integrating a linear spline through the data points. It is now possible to properly deal with non-equidistant data.\n", "\n", "\\\\( A = \\int_{a}^b f(x)\\,\\mathrm{d}x \\approx \\sum_{i=1}^{n-1} \\frac{f(x_{i+1})+f(x_{i})}{2}(x_{i+1}-x_i)\\\\)" ] }, { "cell_type": "code", "execution_count": 10, "id": "77f26c8a", "metadata": {}, "outputs": [], "source": [ "def integralTrapezoidal(data):\n", " a = 0\n", " for i in range( len(data[0]) - 1 ):\n", " a = a + ( data[1,i+1] + data[1,i] ) / 2 * ( data[0,i+1] - data[0, i] )\n", " return a" ] }, { "cell_type": "code", "execution_count": 11, "id": "2e453f2b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.650000000000001" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integralTrapezoidal(data)" ] }, { "cell_type": "markdown", "id": "01ba6e7d", "metadata": {}, "source": [ "For equidistant data this simplfies to \n", "\n", "\\\\( A \\approx \\frac{b-a}{n-1}\\left[\\frac{1}{2}f(x_1) + \\sum_{i=2}^{n-1} f(x_i) + \\frac{1}{2}f(x_n)\\right] \\\\)\n", "\n", "Therefore the edge issue is resolved." ] }, { "cell_type": "code", "execution_count": 12, "id": "d22b56ec", "metadata": {}, "outputs": [], "source": [ "def integralTrapezoidalEQ(data):\n", " return ( 1/2*data[1,0] + np.sum(data[1,1:-1]) + 1/2*data[1,-1] ) * (data[0,-1] - data[0,0]) / ( len(data[1]) - 1 )" ] }, { "cell_type": "code", "execution_count": 13, "id": "57ef7406", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.650000000000001" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integralTrapezoidalEQ(data)" ] }, { "cell_type": "markdown", "id": "02671520", "metadata": {}, "source": [ "## 3. Simpson rule & Newton-Cortes equations" ] }, { "cell_type": "markdown", "id": "92cb777a", "metadata": {}, "source": [ "This method corresponds to integrating a polynomial interpolation function through the data points. The coefficients have been optimized accordingly." ] }, { "cell_type": "markdown", "id": "c32f84f7", "metadata": {}, "source": [ "### Simpson rule" ] }, { "cell_type": "markdown", "id": "593b3861", "metadata": {}, "source": [ "Careful! This method works only if there are an odd number of data points.\n", "\n", "\\\\( A = \\int_{a}^b f(x)\\,\\mathrm{d}x \\approx \\frac{b-a}{n-1} \\left[ \\frac{1}{3}f(x_1) + \\frac{4}{3}f(x_2) + \\frac{2}{3}f(x_3) + \\frac{4}{3}f(x_4) + \\dots + \\frac{4}{3}f(x_{n-3}) + \\frac{2}{3}f(x_{n-2}) + \\frac{4}{3}f(x_{n-1}) + \\frac{1}{3}f(x_n) \\right]\\\\)" ] }, { "cell_type": "code", "execution_count": 14, "id": "ce363749", "metadata": {}, "outputs": [], "source": [ "def integralSimpson(data):\n", " return ( 1/3*data[1,0] + 4/3*np.sum(data[1,1:-1:2]) + 2/3*np.sum(data[1,2:-1:2]) + 1/3*data[1,-1] ) \\\n", " * (data[0,-1] - data[0,0]) / ( len(data[1]) - 1 )" ] }, { "cell_type": "code", "execution_count": 15, "id": "32b032ed", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.599999999999999" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "integralSimpson(data)" ] }, { "cell_type": "markdown", "id": "a48aebdc", "metadata": {}, "source": [ "We get the perfect result, because our data has been generated using a 3rd-order polynomial. This is also what the Simpson rule considers for the interpolation." ] }, { "cell_type": "markdown", "id": "44ecdd9f", "metadata": {}, "source": [ "### Higher-order Newton-Cortes equations " ] }, { "cell_type": "markdown", "id": "3d64d27b", "metadata": {}, "source": [ "There exist many more optimized methods, e.g. Boole's rule:\n", "\n", "\\\\( A = \\int_{a}^b f(x)\\,\\mathrm{d}x \\\\\\approx \\frac{b-a}{n-1} \\left[ \\frac{14}{45}f(x_1) + \\frac{64}{45}f(x_2) + \\frac{24}{45}f(x_3) + \\frac{64}{45}f(x_4) + \\frac{28}{45}f(x_5) + \\frac{64}{45}f(x_6) + \\frac{24}{45}f(x_7) \\dots + \\frac{64}{45}f(x_{n-3}) + \\frac{24}{45}f(x_{n-2}) + \\frac{64}{45}f(x_{n-1}) + \\frac{14}{45}f(x_n) \\right]\\\\)\n", "\n", "For this rule, the number of data points has to be a multiple of 5.\n", "\n", "More information: https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }