{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### In class demo from Spring 2018\n",
"\n",
"\n",
"*(Note: If you are an SFU student you can download this notebook and open it in https://sfu.syzygy.ca, then you can execute the code and play around with it. You'll need to login with your SFU Computing ID.)*\n",
"\n",
"\n",
"**Goals:**\n",
"* On assignment 5 we saw that to approximate the integral\n",
"$$\\int_{0}^{\\pi/2}\\frac{\\sin{x}}{x}~dx$$\n",
"to $15$ decimal places using the midpoint rule, it would take a value of $n$ around the order of $40$ million. Suggesting the Midpoint Rule won't be useful. This demo will show a creative way to get the Midpoint Rule to work for us, with relatively small $n$.\n",
"\n",
"* We'll see graphically the rate of convergence for Midpoint and Simpons's rule as $n$ increases, by plotting errors-vs-n on a log-log scale.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementation of Approximation Methods\n",
"\n",
"First we implement the approximation methods that we will use:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Midpoint Rule: ##\n",
"\n",
"The function `midpoint_rule` outputs the midpoint rule approximation of $\\int_a^b f(x) dx$ using n partitions of the interval:\n",
"\n",
"$$\\int_a^bf(x)dx \\approx \\sum_{i=1}^n f(\\overline{x}_{i-1})\\Delta x,$$\n",
"where $\\overline{x}_i$ is the midpoint of inteval $[x_{i-1},x_i]$, which is $\\overline{x}_i=\\frac{x_{i-1}+x_i}{2}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.45988929071851814\n"
]
}
],
"source": [
"import math\n",
"\n",
"def midpoint_rule(fcn,a,b,n):\n",
"# output: midpoint rule approx of int_a^b fcn(x)dx using n steps\n",
" Deltax = (b-a)*1.0/n\n",
" xs=[a+Deltax*i for i in range(n+1)]\n",
" ysmid=[fcn((xs[i]+xs[i+1])/2) for i in range(n)]\n",
" return Deltax*sum(ysmid)\n",
"\n",
"#Example\n",
"n=10;\n",
"a=0\n",
"b=1\n",
"def f(x):\n",
" return(math.sin(x))\n",
"print(midpoint_rule(f,a,b,n));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Trapezoid Rule: ##\n",
"The function `trapezoid_rule` outputs the trapezoid rule approximation of $\\int_a^b f(x) dx$ using n partitions of the interval:\n",
"\n",
"$$\\int_a^bf(x)dx \\approx \\frac{\\Delta x}{2}(f(x_0)+2f(x_1)+2f(x_2)+ \\cdots + 2f(x_{n-1})+f(x_n)).$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.4593145488579764\n"
]
}
],
"source": [
"import math\n",
"\n",
"def trapezoid_rule(fcn,a,b,n):\n",
"# output: trapezoid rule approx of int_a^b fcn(x)dx using n steps\n",
" Deltax = (b-a)*1.0/n\n",
" coeffs = [2]*(n-1)\n",
" coeffs = [1]+coeffs+[1]\n",
" valsf = [fcn(a+Deltax*i) for i in range(n+1)]\n",
" return (Deltax/2)*sum([coeffs[i]*valsf[i] for i in range(n+1)])\n",
"\n",
"#Example\n",
"n=10;\n",
"a=0\n",
"b=1\n",
"def f(x):\n",
" return(math.sin(x))\n",
"print(trapezoid_rule(f,a,b,n));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simpsons Rule: ##\n",
"\n",
"The function `simpsons_rule` outputs the Simpson's Rule approximation of $\\int_a^b f(x) dx$ using n partitions\n",
"of the interval:\n",
"\t\t\t\t\n",
"$$\\int_a^bf(x)dx \\approx \\frac{\\Delta x}{3}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+ \\cdots + 2f(x_{n-2})+4f(x_{n-1})+f(x_n)).$$\n",
"\t\t\t\t\n",
"**Note:** n must be even."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.4596979498238206\n"
]
}
],
"source": [
"import math\n",
"\n",
"def simpsons_rule(fcn,a,b,n):\n",
"# output: simpsons rule approx of int_a^b fcn(x)dx using n steps (n must be an even integer)\n",
" Deltax = (b-a)*1.0/n\n",
" n2=int(n/2)\n",
" coeffs = [4,2]*n2\n",
" coeffs = [1] +coeffs[:n-1]+[1]\n",
" valsf = [fcn(a+Deltax*i) for i in range(n+1)]\n",
" return (Deltax/3)*sum([coeffs[i]*valsf[i] for i in range(n+1)])\n",
"\n",
"#Example\n",
"n=10;\n",
"a=0\n",
"b=1\n",
"def f(x):\n",
" return(math.sin(x))\n",
"print(simpsons_rule(f,a,b,n));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Piecewise Defined Functions: ###\n",
"\n",
"Sometimes you may find the need to define $f$ as a piecewise defined function. For example, suppose you are trying to approximate the integral $\\int_{0}^{10} f(x)~dx$ where\n",
"$$f(x)=\\begin{cases}\n",
"\\displaystyle x,&x\\le 5\\\\\n",
"x^2,&x > 5\\\\\n",
"\\end{cases}\n",
"$$ (not that you would not need to approximate in this case since you can easily find antiderivates, but we'll just use it as a simple example).\n",
"\n",
"\n",
"We can replace the line `f(x) = ...` in the code above with a python function defining $f$:\n",
"\n",
"```python\n",
"def f(x):\n",
" if x<=5:\n",
" return x;\n",
" elif x>5: # elif means \"else if\"\n",
" return x^2;\n",
"```\n",
"\n",
"Try to cut-and-paste this code to replace the function $f$ in Simpson's Rule above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Class Demonstration:\n",
"\n",
"Let's investigate a 'smarter' midpoint rule for $$f(x)=\\frac{\\sin{x}}{x}.$$\n",
"\n",
"That is, we want a method that will give us a very good approximation with a relatively small value for $n$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# we'll need to import a few mathematical libraries\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from scipy import integrate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Richardson Extrapolation with Midpoint Rule:\n",
"\n",
"Using the midpoint rule to approximate the following integral\n",
"$$\\int_{0}^{\\pi/2} \\frac{\\sin(x)}{x} ~dx$$\n",
"to $15$ decimal places requires $n$ to be significantly large. Error bound estimates suggest $n$ to be at least 20 million.\n",
"\n",
"However, using Richardson Extrapolation we can take a weighted average of $M_n$ and $M_{2n}$, for $n$ relatively small (i.e. $n=2000$) and get $15$ decimal places.\n",
"\n",
"Let $I_{exact}$ denote the exact value of the integral, and $Err_n$ the error in the n-th approximation. Then\n",
"\n",
"\\begin{align*}\n",
"I_{exact} &= M_{n} + Err_n\\\\\n",
"I_{exact} &= M_{2n} + Err_{2n} \\simeq M_{2n} + \\frac{Err_{n}}{4}\n",
"\\end{align*}\n",
"\n",
"Therefore,\n",
"\\begin{align*}\n",
"4I_{exact} - I_{exact} & = 4\\left(M_{2n} + \\frac{Err_{n}}{4}\\right) - (M_{n} + Err_{n})\\\\\n",
"3I_{exact} & = 4M_{2n} - M_{n}\n",
"\\end{align*}\n",
"and so\n",
"$$I_{exact} \\simeq \\left(\\frac{4}{3}\\right)M_{2n} + \\left(\\frac{-1}{3}\\right)M_{n}.$$\n",
"\n",
"This gives a better approximation of $I_x$ than just using $n$ and $2n$ alone. \n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def f(x):\n",
" if x==0:\n",
" return(1)\n",
" else:\n",
" return(np.sin(x)/x)\n",
"a,b = 0, np.pi/2"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3707621681544886"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=2000\n",
"(4/3)*midpoint_rule(f,a,b,2*n)-(1/3)*midpoint_rule(f,a,b,n)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.370762168154486"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"simpsons_rule(f,a,b,2000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Comparing with the actual value below we see we have obtained $16$ decimal places."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3707621681544881"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate.quad(f,a,b)[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Richardson Extrapolation with Simpson's Rule:\n",
"We could do the same analysis using Simpson's rule. In this case we use the fact that doubling $n$ reduces that error by a factor of $2^4 = 16$. This gives\n",
"\n",
"$$I_{exact} \\simeq \\left(\\frac{16}{15}\\right)S_{2n} + \\left(\\frac{-1}{15}\\right)S_{n}.$$\n",
"\n",
"The next line shows that we can more than 16 decimal places of accuracy using $n$ and $2n$, where $n=100$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3707621681544886"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=100\n",
"(16/15)*simpsons_rule(f,a,b,2*n)-(1/15)*simpsons_rule(f,a,b,n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Convergence rate for Approximation Methods\n",
"\n",
"As $n$ increases we expect $S_n$ and $M_n$ to get closer to the actual values. Can we see the rate at which these approximations converge to the actual values?\n",
"\n",
"From the error bound estimate for Simpon's rule \n",
"\n",
"$$|E_S| \\le \\frac{k(b-a)^5}{180n^4} $$\n",
"\n",
"we infer that $E_S$ is proportional to $n^{-4}$, therefore plotting $\\log{E_S}$ versus $\\log{n}$ should produce something close to a line of slope $-4$. We now construct this plot to see this."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def f(x): # define our function f\n",
" return(np.exp(np.cos(x))) # use exp and cos from numpy's math library \n",
"a,b = 0,1 # set the limits of integration"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"Ival=integrate.quad(f,a,b)[0] # this is what we'll use for the *exact* value of the integral"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEOCAYAAAC0BAELAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+YZFdZoN9v0nSwY8UOFdJJm85U6KSbdFJJrT3Lj0VBAkFEA0JAwzNiNGZno7LrorgxO+6K4qwCCq7rDxyDjwFaAcGQEIKQmE4iSuJOb5qpSYdupkklHTo2UDPNlDaTtsm3f9xbTaXTP6rm9q1b55vvfZ56quqcc+/93jn9zFf33nPPEVXFcRzHcbabHVkH4DiO49jEE4zjOI6TCp5gHMdxnFTwBOM4juOkgicYx3EcJxU8wTiO4zip4AnGcRzHSQVPMI7jOE4qeIJxHMdxUqEr6wCy5Mwzz9RCoXDC2z/55JOceuqp2xdQB2HVzaoXuFuohOg2MTHxDVV97lbtTuoEUygUOHDgwAlvf+zYMU4//fRtjKhzsOpm1QvcLVRCdBORR5tp55fIEjAxMZF1CKlh1c2qF7hbqFh28wSTgJ6enqxDSA2rbla9wN1CxbKbJ5gEJLl/0+lYdbPqBe4WKpbdPMEkYGpqKusQUsOqm1UvcLdQsezmCSYBln95WHWz6gXuFiqW3TzBtMrYGBQKsGMH3/uSl0TfDVKr1bIOIRWseoG7hYplt5N6mHLLjI3Bnj2wtARA9xNPRN8Bdu/OMLDtp1qtZh1CKlj1AncLFctucjIvmbxr1y5t6TmYQgEeXWf4986dUKlsV1gdQYhj85vBqhe4W6iE6CYiE6q6a6t2fomsFR57rLXygLE6Nt+qF7hbqFh28wTTCued11p5wORyuaxDSAWrXuBuoWLZzRNMK+zbB2sfiurpicqN0d/fn3UIqWDVC9wtVCy7eYJphd27Yf/+6J6LCN8666zou7Eb/ADT09NZh5AKVr3A3ULFspuPImuV3btXE8rXH3uM8wxeHgMYHBzMOoRUsOoF7hYqlt38DCYBlocXWnWz6gXuFiqW3TzBJGBxcTHrEFLDqptVL3C3ULHs5gkmAaOjo1mHkBpW3ax6gbuFimU3TzAJsDx+3aqbVS9wt1Cx7NYxCUZE3iQiD4nIUyKya03djSJyWESmReSHNtj+fBF5QES+LCIfFZHutGPu7e1N+xCZYdXNqhe4W6hYduuYBAMcAt4A3NdYKCIjwNXAxcCrgT8WkVPW2f5dwPtU9ULgKPCz6YYL+Xw+7UNkhlU3q17gbqFi2a1jEoyqPqyq6w0Ifx3wEVV9UlUfAQ4DL2hsICICXA58PC66GfixNOMFmJ2dTfsQmWHVzaoXuFuoWHbrmASzCd8LzDV8fzwuayQPLKrqyiZttp3h4eG0D5EZ2+bWsLwBhULmyxt4n4WJu4VJWx+0FJG7gLPXqdqrqrdutNk6ZWungG6mTT2GPcAeiKZoGB8fZ2RkhEqlwtLSEqOjo0xMTNDX10d3dzdzc3MUi0VmZmZYWVmhWCwyOTlJf38/hw8fZnp6mlKpRLlcpquri6GhIcrlMgMDAywvL7OwsLC6z56eHgqFAlNTUxQKBWq1GtVqdbU+l8vR39/P9PQ0g4ODVKtVFhcXV+t7e3vJ5/PMzs4yPDzM/Pw8tVpttT6fz5PL5ahUKifsBDA/P88pp5zC4cOHEznNvetdPP+97+WU48ejf/xHH+Xb113Hl6amGLjhhrY7lUol7r33Xnbu3Gmmnxr/9paXl5menjblVO+nb33rWywuLppyqvfTgw8+yJVXXhmUU9Ooake9gHuAXQ3fbwRubPj+WeDFa7YR4BtAV/z9xcBntzrW6OioJuHuu+9OtH0nsy1uO3eqwjNfO3cm3/cJ4n0WJu7WWQAHtIn/z0O4RHYbcLWInCoi5wMXAv/U2CAWHgfeGBddA2x0RrRtWB6/vi1uHbi8gfdZmLhbmHRMghGR14vI40RnH58Wkc8CqOpDwMeAKeBvgV9Q1W/H29whIvVzthuAXxKRw0T3ZD6QdsyWx69vi1sHLm/gfRYm7hYmHTPZpareAtyyQd0+4Blz4qvqaxo+f4U1o8vSxvLwwm1x27fvaUtMA5kvb+B9FibuFiYdcwYTIpYXCtoWtzXLG7BzZ+bLG3ifhYm7hYknmARUKpWsQ0iNbXPbvRsqFXjqqeg947VzvM/CxN3CxBNMAkZGRrIOITWsuln1AncLFctunmASYPmXh1U3q17gbqFi2c0TTAKWGm9eG8Oqm1UvcLdQsezmCSYBlsevW3Wz6gXuFiqW3TzBJMDy+HWrbla9wN1CxbKbJ5gE9PX1ZR1Calh1s+oF7hYqlt08wSSguzv1Nc0yw6qbVS9wt1Cx7OYJJgFzc3NbNwoUq24d75VgeYOOd0uAu4VJx0wVEyLFYjHrEFLDqltHe42NPX1qnUcfjb5DUw+odrRbQtwtTPwMJgEzMzNZh5AaVt062mvv3qfP2wbR9717m9q8o90S4m5h4gkmASsrK1s3ChSrbh3tlXB5g452S4i7hYknmARYPrW16tbRXgmXN+hot4S4W5h4gknA5ORk1iGkhlW3jvbaty9azqCRFpY36Gi3hLhbmHiCSUDL61MHhFW3jvZKuLxBR7slxN3CxEeROU4nsXt35ksaOM524WcwCZifn886hNSw6mbVC9wtVCy7dUSCEZE3ichDIvKUiOxqKL9CRCZEpBy/X77B9u8Qka+KyGT8es167babUqnUjsNkglU3q17gbqFi2a0jEgxwCHgDcN+a8m8AV6pqEbgG+NAm+3ifqpbi1x0pxfk0yuVyOw6TCVbdrHqBu4WKZbeOuAejqg8DiMja8gcbvj4EPFtETlXVJ9sY3oZ0dXXEP18qWHWz6gXuFiqW3TrlDKYZrgIe3CS5vFVEDorIn4vIGe0IaGhoqB2HyQSrbla9wN1CxbJb21KniNwFnL1O1V5VvXWLbS8G3gW8aoMmfwK8E9D4/feAazfY1x5gD0TDA8fHxxkZGaFSqbC0tMTo6CgTExP09fXR3d3N3NwcxWKRmZkZVlZWKBaLTE5O0t/fz8GDBznzzDMplUqUy2W6uroYGhqiXC4zMDDA8vIyCwsLq/vs6emhUCgwNTVFoVCgVqtRrVZX63O5HP39/UxPTzM4OEi1WmVxcXG1vre3l3w+z+zsLMPDw8zPz1Or1Vbr8/k8uVyOSqVywk4Q3XQ8fvw4p512mimnUqnE7bffztDQkCmn+t/e0aNH6enpMeVU76cjR45w0UUXmXKq99MDDzzAVVddFZRTs4iqtrRBmojIPcDbVfVAQ9m5wN3Az6jqPzSxjwJwu6peslXbXbt26YEDB7ZqtiGHDx/mggsuOOHtOxmrbla9wN1CJUQ3EZlQ1V1btevoS2Qi0gt8Grhxs+QiIuc0fH090aCB1FleXm7HYTLBqptVLzhJ3RIsb9ApWO63jkgwIvJ6EXkceDHwaRH5bFz1VuAC4H80DEE+K97mpoYhze+OhzIfBF4OvK0dcS8sLLTjMJlg1c2qF5yEbvXlDR59FFS/s7xBYEnGcr911CWydpP0EtmxY8c4/fTTtzGizsGqm1UvOAndCoUoqaxl506oVNoR1rYQYr+ZuETW6UxMTGQdQmpYdbPqBSehW8LlDToFy/3mCSYBPWtnvjWEVTerXnASuiVc3qBTsNxvnmASUCgUsg4hNay6WfWCk9At4fIGnYLlfvMEk4CpqamsQ0gNq25WveAkdEu4vEGnYLnf7M5R0AYs//Kw6mbVC05SNwPLG1juNz+DSUCtVss6hNSw6mbVC9wtVCy7eYJJQLVazTqE1LDqZtUL3C1ULLt5gknA6Oho1iGkhlU3q17gbqFi2c0TTAIsj1+36mbVC9wtVCy7eYJJQC6XyzqE1LDqZtUL3C1ULLt5gklAq1NXh4RVN6te4G6hYtnNE0wCpqensw4hNay6WfUCdwsVy26eYBIwODiYdQipYdXNqhe4W6i03a2NSxx4gkmA5eGFVt2seoG7hUpb3dq8xIEnmAQsLi5mHUJqWHWz6gXuFiptddu7F5aWnl62tBSVp4CvB+PrwayLVTerXuBuodJWtx07ojOXtYjAU081vRtfD6YNWB6/btXNqhe4W6i01a3NSxxsOdmliDR75EVVPZYwnqDo7e3NOoTUsOpm1QvcLVTa6rZvX3TPpfEyWYpLHDQzm/LNgAKySRsF/gL44IkEISJvAt4BXAS8QFUPxOUF4GGgPo7vflW9fp3tnwN8FCgAFeDHVfXoicTSCvl8Pu1DZIZVN6te4G6h0la3+szTe/dGK3+ed16UXFKakXrLBKOqL0/lyE/nEPAG4E/XqZtV1dIW2/8q8Heq+jsi8qvx9xu2OcZnBjY7y3mBrZ7XLFbdrHqBu4VK293auMRBR9yDUdWHVTXJ00avIzrTIn7/seRRbc3w8HA7DpMJVt2seoG7hYpltxNKMCLyERH5UPx693YHtYbzReRBEblXRH5ggzZ9qvoEQPx+VsoxATA/P9+Ow2SCVTerXuBuoWLZ7URXtPyCqv5vABFp6gKiiNwFnL1O1V5VvXWDzZ4AzlPVqoiMAp8UkYuTDCYQkT3AHojmABofH2dkZIRKpcLS0hKjo6NMTEzQ19dHd3c3c3NzFItFZmZmWFlZoVgsMjk5SX9/P4888gi1Wo1SqUS5XKarq4uhoSHK5TIDAwMsLy+zsLCwus+enh4KhQJTU1MUCgVqtRrVanW1PpfL0d/fz/T0NIODg1SrVRYXF1fre3t7yefzzM7OMjw8zPz8PLVabbU+n8+Ty+WoVCon7ATRH/zx48e57777TDmVSiUOHTrE8vKyKaf6397Ro0cZHx835VTvpyNHjtDT02PKqd5P5XKZoaGhoJya/v/2RJ6DEZG7gVuAz6rqTMs72Hi/9wBvr9/kb7ZeRKaBH1TVJ0TkHOAeVd3yvNOfg9kYq25WvcDdQiVEt7Sfg3kLMAtcJSI3neA+tkREnisip8SfnwdcCHxlnaa3AdfEn68BNjoj2lZ8bH54WPUCdwsVy25NJxgR+X0REQBV/aqq3qGqv62q1yUNQkReLyKPAy8GPi0in42rXgocFJEvAh8HrlfVI/E2N4lIPYP+DnCFiHwZuCL+njo+dDI8rHqBu4WKZbdW7sH8C3CbiFytqv8qIq8Cfl1VX5I0CFW9heiS29ryTwCf2GCb6xo+V4FXJI2jVSwvFGTVzaoXuFuoWHZr+gxGVX8N+CvgHhH5PPDLRM+bnLRUKpWsQ0gNq25WvcDdQsWyW9NnMCLyCuA/Av8KnAP8bMJnV4JnZGQk6xBSw6qbVS9wt1Cx7NbKTf69wP9U1R8E3gh8VEQuTyWqQLD8y8Oqm1UvcLdQsezW9BmMql7e8LksIj9MdH/kP6QRWAgsrV1XwRBW3ax6gbuFimW3VkaRPUdE3ikifyoivwgcJ4Mb61nTuNroNde8LM3VRjNldHQ06xBSwaoXuFuoWHZr5RLZR4Aa8CmgB/g8UEwjqE5l7Wqjc3M70lxtNFOsjs236gXuFiqW3VpJMOeo6rtV9XZV/W3gSuAPUoqrI2nzaqOZ0tfXl3UIqWDVC9wtVCy7tZJgjojIpfUvqvoVojOZk4bHHmutPGS6u7uzDiEVrHqBu4WKZbdWEswe4C9F5E9E5OdF5A+Jpos5aWjzaqOZMjc3l3UIqWDVC9wtVCy7NZVgRGQHcBXwfcA40XT4XwTenF5once+fdHqoo2kuNpophSLNm+vWfUCdwsVy25NJRhVfQp4paouq+rHVPUdqvpnqno85fg6it27Yf9+2LkTRKCv7zj797dtcbi2MjOzbZNkdxRWvcDdQsWyWyuXyB4UkV+vT3h5srJ7N1Qq8NRT8LGP/ZPJ5AKwsrKSdQipYNUL3C1ULLu1MtnlANGw5J8TkQeAg8BBVf3rVCILAMuntlbdrHqBu4WKZbdWJrv8cVW9CNgJ/AZwGHhhWoGFwOTkZNYhpIZVN6te4G6hYtmtlSf5f1hE7ic6c/lVYEZV355aZAHQ6vKhIWHVzaoXuFuoWHZr5R7MHxNN0f8iYD/wHhE5qUaROY7jOM3TSoJZUNV/UNWjqnoX8ENEMyyftMzPz2cdQmpYdbPqBe4WKpbdWkkwFRH5LRGpP3b6b0Rzk520lEqlrENIDatuVr3A3ULFslsrCUaBNwBz8YqWh4lWt7wwaRAi8iYReUhEnhKRXQ3lu0VksuH1lIg8ozdE5B0i8tWGdq9JGlMzlMvldhwmE6y6WfUCdwsVy26trAfzZgAReTZwCXBZ/LpJRJ6nqgMJ4jhElLz+dM0xx4Cx+LhF4FZV3WjIxftU9XcTxNAyXV2tjPIOC6tuVr3A3ULFtFurG8RP7x+IX9uCqj4MsMUznG8G/mq7jrkdDA0NZR1Calh1s+oF7hYqlt2aTjDxCpbvAHqJ5iF7r6ren1Jc6/ETwOs2qX+riPwUUeL7ZVU9ul4jEdlDNHEn/f39jI+PMzIyQqVSYWlpidHRUSYmJujr66O7u5u5uTmKxSIzMzOsrKxQLBaZnJykv7+fgwcPcuaZZ1IqlSiXy3R1dTE0NES5XGZgYIDl5WUWFhZW99nT00OhUGBqaopCoUCtVqNara7W53I5+vv7mZ6eZnBwkGq1yuLi4mp9b28v+Xye2dlZhoeHmZ+fp1arrdbn83lyuRyVSuWEnSC66Xj8+HFOO+00U06lUonbb7+doaEhU071v72jR4/S09NjyqneT0eOHOGiiy4y5VTvpwceeICrrroqKKdmEVVtrqHII8BPAlPAKFGy+SNVbeqsQkTuAs5ep2qvqt4at7kHeLuqHliz7QuBm1R13UdeRaQP+AbRfaJ3Eq1dc+1WMe3atUsPHDjxE7HDhw9zwQUXnPD2nYxVN6te4G6hEqKbiEyo6q6t2rVyiWxBVf8h/nyXiHwBeIAmL1up6itbONZart7sOKq6UP8sIn8G3J7gWE2zvLzcjsNkglU3q17gbqFi2a3jhynHSwW8iWjJ5o3anNPw9fVEgwZSZ2FhYetGgWLVzaoXuFuoWHbrlGHKrxeRx4EXA58Wkc82VL8UeDxeQbNxm5sahjS/W0TKInIQeDnwtqQxNcPo6Gg7DpMJVt2seoG7hYplt1Ymu3yzqo4QTXb5X4kmvDyNaJhyoiXZVPUWVT1XVU9V1T5V/aGGuntU9UXrbHNd/V6Nqr5FVYuqeqmqvlZVn0gST7NMTEy04zCZYNXNqhe4W6hYdmtlFFke+HHgOPAQ8Jeq+oG0AguBnrXLWxrCqptVL3C3ULHs1solsluA5wL/C3gP8E0ReTiVqAKhUChkHUJqWHWz6gXuFiqW3VpJMDlV/U2i0WQvI3rw8S9SiSoQpqamsg4hNay6WfUCdwsVy26tJJjj8fuTIvJdqvoJoC1zfnUqln95WHWz6gXuFiqW3Vp5DuZ3ReQ5wEeBPxeRfwS+N52wwqBWszuZtFU3q17gbqFi2a2VUWSfUNUjqvpe4A5ggM2nbjFPtVrNOoTUsOpm1QvcLVQsuzU9VYxFkk4Vc+zYMU4//fRtjKhzsOpm1QvcLVRCdGt2qphW7sE4a7A8fn273MbGoFCAHTui97GxbdntCeN9FibuFiYtJxgRuTKNQEIkl8tlHUJqbIfb2Bjs2QOPPgqq0fuePdkmGe+zMHG3MDmRM5h92x5FoLQ6dXVIbIfb3r2wtPT0sqWlqDwrvM/CxN3C5EQSzKargp1MTE9PZx1CamyH22OPtVbeDrzPwsTdwuREEszJOypgDYODg1mHkBrb4Xbeea2VtwPvszBxtzDxm/wJsDy8cDvc9u2DtdMs9fRE5VnhfRYm7hYmnmASsLi4mHUIqbEdbrt3w/79sHMniETv+/dH5VnhfRYm7hYmLT8HIyJ3quoVKcXTVvw5mI2x6mbVC9wtVEJ0S+05GCvJZTuwPH7dqptVL3C3ULHs5pfIEtDb25t1CKlh1c2qF7hbqFh28wSTgHw+n3UIqWHVzaoXuFuoWHbbMsGIyB+LyM+JyEtEJLULhSLyHhH5kogcFJFbRKS3oe5GETksItMi8kMbbH++iDwgIl8WkY+KSHdasdaZnZ1N+xCZYdXNqhe4W6hYdmvmDGYSuJhoJctHRKQiIp8SkX0icvU2xnIncImqXgrMADcCiMgIcHUcw6uBPxaRU9bZ/l3A+1T1QuAo8LPbGNu6DA8Pp32IzLDqZtUL3C1ULLttmWBUdb+qvlVVX6aqeeD7gT8BjgE/sl2BqOrnVHUl/no/cG78+XXAR1T1SVV9BDgMvKBxWxER4HLg43HRzcCPbVdsGzE/P5/2ITLDqptVL3C3ULHstuWCYyKy3nPXh+JXY/2iqh7bpriuJVrYDKJFze5vqHucZy50lo+Pv7JJGwBEZA+wB6I5gMbHxxkZGaFSqbC0tMTo6CgTExP09fXR3d3N3NwcxWKRmZkZVlZWKBaLTE5O0t/fzyOPPEKtVqNUKlEul+nq6mJoaIhyuczAwADLy8ssLCys7rOnp4dCocDU1BSFQoFarUa1Wl2tz+Vy9Pf3Mz09zeDgINVqlcXFxdX63t5e8vk8s7OzDA8PMz8/T61WW63P5/PkcjkqlcoJO0H0B3/8+HHuu+8+U06lUolDhw6xvLxsyqn+t3f06FHGx8dNOdX76ciRI/T09JhyqvdTuVxmaGgoKKdm2fI5GBEZ36CqvqHEn/9CVT+4xb7uAs5ep2qvqt4at9kL7ALeoKoqIn8EfEFVPxzXfwC4I16yub7f58ZtLoi/D8RtipvF48/BbIxVN6te4G6hEqLbtj0Ho6ov3+B1efyqf940ucT7eqWqXrLOq55crgF+FNit38l8jxOtnlnnXGDtOeU3gF4R6dqkzbZjefy6VTerXuBuoWLZrWOGKYvIq4EbgNeqauMk77cBV4vIqSJyPnAh8E+N28bJaBx4Y1x0DXBr2jFbHl5o1c2qF7hbqFh2azrBxMN//0ZEfl1EXicihW2O5Q+BHHCniEyKyPsBVPUh4GPAFPC3wC+o6rfjmO4QkfpFwRuAXxKRw0T3ZD6wzfE9A8sLBVl1s+oF7hYqlt1aOYP5U+CfgSrww8AhESmLyG+KyLOSBqKqF6jqgKqW4tf1DXX7VHVQVYdV9TMN5a9R1fn481dU9QXxft6kqk8mjWkrKpVK2ofIDKtuVr3A3ULFstuWo8ga+ElVLdW/xGcYP0M0XPm9wH/e5tg6npGRkaxDSA2rbla9wN1CxbJbK2cw3xSRS+tfVHUSeJGq/i7wkm2PLAAs//Kw6mbVC9wtVCy7tXIGcz3wYRGZJHq6fxh4Kq5LfVqWTmRp7YLzhrDqZtUL3C1ULLs1fQajqg8TPUH/t8BZRE/U/6iInAZ8JJ3wOpvR0dGsQ0gNq25WvcDdQsWyWyujyJ4DvAN4JdFzJzeralVV/1VVfyul+Doay+PXrbpZ9QJ3CxXLbq3cg/kIUAM+BfQAnxeRF2y+iW36+vqyDiE1rLpZ9QJ3CxXLbq3cgzlHVd8df75dRD4K/CXwou0PKwy6u+3eerLqZtUL3C1ULLu1cgZzZM0osq8QncmctMzNzWUdQmpYdbPqBe4WKpbdWjmD2QN8QkT+HigTrc9id6WcJigWN51LM2isuln1AncLFctuzaxo+UER+SWi6e8vJ5rz67nAg8Cb0w2vs5mZmck6hNSw6tbpXmNjUCjAjh3R+9hY89t2ulsS3C1MmjmDuRm4jGgCycuA04nmBXsWcCXw16lF1+GsrKxs3ShQrLp1stfYGOzZA/XHIh59NPoOsHv31tt3sltS3C1MtlwP5hkbRFPijxAlm8tU9e1pBNYOkq4Hc/ToUc4444xtjKhzsOrWyV6FQpRU1rJzJzTzsHcnuyXF3TqLbVsPZi2quqKqB1X1QyEnl+1gcnIy6xBSw6pbJ3s99lhr5WvpZLekuFuYdMx6MCHS6vKhIWHVrZO9zltvcfJNytfSyW5Jcbcw8QTjOB3Cvn3Qs2bgf09PVO44IeIJJgHz86mvypwZVt062Wv3bti/P7rnIhK979/f3A1+6Gy3pLhbmLR8k98SfpN/Y6y6WfUCdwuVEN1Su8mfBiLyHhH5kogcFJFbRKQ3Lr9CRCbilTMnROTyDbZ/h4h8NV5qeVJEXtOOuMvlcjsOkwlW3ax6gbuFimW3jkgwwJ3AJap6KTAD3BiXfwO4UlWLRM/hfGiTfbyvYbnlO9INN6Krq5WJEMLCqptVL3C3ULHs1hEJRlU/p6r1p43uB86Nyx9U1foFyoeAZ4vIqVnEuB5DQ0NZh5AaVt2seoG7hYplt45IMGu4FvjMOuVXAQ+q6pMbbPfW+BLbn4tIWy5oWj61tepm1QvcLVQsu7XtJr+I3AWcvU7VXlW9NW6zF9gFvEEbAhORi4HbgFep6jMm2BSRPqLLaQq8k2hpgWs3iGMP0cSd9Pf3j374wx9mZGSESqXC0tISo6OjTExM0NfXR3d3N3NzcxSLRWZmZlhZWaFYLDI5OUl/fz9PPPEEqkqpVKJcLtPV1cXQ0BDlcpmBgQGWl5dZWFhY3WdPTw+FQoGpqSkKhQK1Wo1qtbpan8vl6O/vZ3p6msHBQarVKouLi6v1vb295PN5ZmdnGR4eZn5+nlqttlqfz+fJ5XJUKpUTdoJoVMsZZ5zBsWPHTDmVSiXuvfdezjrrLFNO9b+9U089lWPHjplyqvfTKaecwplnnmnKqd5PDz/8MFdccUVQTs9//vObusnfMaPIROQa4HrgFaq61FB+LnA38DOq+g9N7KcA3K6ql2zVNukosqmpKUZGRk54+07GqptVL3C3UAnRLbRRZK8GbgBeuya59AKfBm7cLLmIyDkNX18PHEor1kYWFhbacZhMsOpm1QvcLVQsu3XEGYyIHAZOBapx0f2qer2I/BrRiLIvNzR/lap+TURuAt6vqgdE5ENAiegSWQX4T6r6xFbHTXoGc+zYMU4//fQT3r6Tsepm1QvcLVRCdAvqDEZVL1DVgYZhxtfH5b+lqqc1lJdU9Wtx3XWqeiD+/BZVLarqpar62maSy3YwMTHRjsNkglU3q17gbqFi2a0jEkyo9KydOMpnaA/DAAAO8klEQVQQVt2seoG7hYplN08wCSgUClmHkBpW3ax6gbuFimU3TzAJmJqayjqE1LDqZtUL3C1ULLt5gkmA5V8eVt2seoG7hYplN08wCajValmHkBpW3ax6gbuFimU3TzAJqFarWzcKFKtuVr3A3ULFspsnmASMjo5mHUJqWHWz6gXuFiqW3TzBJMDy+HWrbla9wN1CxbKbJ5gE5HK5rENIDatuVr3A3ULFspsnmATUZxa1iFU3q17gbqFi2c0TTAKmp6ezDiE1rLpZ9QJ3CxXLbp5gEjA4OJh1CKlh1c2qF7hbqFh28wSTAMvDC626WfUCdwsVy26eYBKwuLiYdQipYdXNqhecnG5jY1AowI4d0fvYWFvD2hYs91tHrAeTFb4ezMZYdbPqBSef29gY7NkDS0vfKevpgf37YffuNgeYgBD7Laj1YELF8vh1q25WveDkc9u79+nJBaLve/e2KahtwnK/eYJJQG9vb9YhpIZVN6tecPK5PfbY+m03Ku9ULPebJ5gE5PP5rENIDatuVr3g5HM777z1225U3qlY7reOSDAi8h4R+ZKIHBSRW0SkNy4viMi3RGQyfr1/g+2fIyJ3isiX4/cz2hH37OxsOw6TCVbdrHrByee2b190z6WRnp6oPCQs91tHJBjgTuASVb0UmAFubKibVdVS/Lp+g+1/Ffg7Vb0Q+Lv4e+oMDw+34zCZYNXNqhecfG67d0c39HfuBJHoPbQb/GC73zoiwajq51R1Jf56P3Bui7t4HXBz/Plm4Me2K7bNmJ+fb8dhMsGqm1UvODnddu+GSgWeeip6Dy25gO1+64gEs4Zrgc80fD9fRB4UkXtF5Ac22KZPVZ8AiN/PSjtIsL1QkFU3q17gbqFi2a2rXQcSkbuAs9ep2quqt8Zt9gIrQP1xqSeA81S1KiKjwCdF5GJVPZYgjj3AHogmmRsfH2dkZIRKpcLS0hKjo6NMTEzQ19dHd3c3c3NzFItFZmZmWFlZoVgsMjk5SX9/P729vYyPj1MqlSiXy3R1dTE0NES5XGZgYIDl5WUWFhZW99nT00OhUGBqaopCoUCtVqNara7W53I5+vv7mZ6eZnBwkGq1yuLi4mp9b28v+Xye2dlZhoeHmZ+fp1arrdbn83lyuRyVSuWEnSD6RXXBBRdw3333mXIqlUqsrKzwj//4j6ac6n97559/PuPj46ac6v2Uz+eZmpoy5VTvpyeffJJjx44F5dQ0qtoRL+Aa4AtAzyZt7gF2rVM+DZwTfz4HmG7mmKOjo5qEu+++O9H2nYxVN6tequ4WKiG6AQe0if9jO+ISmYi8GrgBeK2qLjWUP1dETok/Pw+4EPjKOru4jShBEb/fmm7EEZaHF1p1s+oF7hYqlt06IsEAfwjkgDvXDEd+KXBQRL4IfBy4XlWPAIjITSJSn6rgd4ArROTLwBXx99SxvFCQVTerXuBuoWLZrSMSjKpeoKoDumY4sqp+QlUvVtXLVPX7VPVTDdtcp6oH4s9VVX2Fql4Yvx9pR9yVSqUdh8kEq25WvcDdQsWyW0ckmFAZGRnJOoTUsOpm1QvcLVQsu3mCSYDlXx5W3ax6gbuFimU3TzAJWFo7lashrLpZ9QJ3CxXLbp5gEjA6Opp1CKlh1c2qF7hbqFh28wSTAMvrOFh1s+oF7hYqlt08wSSgr68v6xBSw6qbVS9wt1Cx7OYJJgHd3d1Zh5AaVt2seoG7hYplN08wCZibm8s6hNSw6mbVC9wtVCy7eYJJQLFYzDqE1LDqZtUL3C1ULLt5gknAzMxM1iGkhlU3q17gbqFi2c0TTAJWVla2bhQoVt2seoG7hYplN08wCbB8amvVzaoXuFuoWHbzBJOAycnJrENIDatuVr3A3ULFspsnmAS0vLpbQFh1s+oF7hYqlt08wTiO4zip4AkmAfPz81mHkBpW3ax6gbuFimU3TzAJKJVKWYeQGlbdrHqBu4VKu93GxqBQgB07ovexsfSO5QkmAeVyOesQUsOqm1UvcLdQaafb2Bjs2QOPPgqq0fuePeklmY5IMCLyHhH5kogcFJFbRKQ3Lt8tIpMNr6dE5BnpXkTeISJfbWj3mnbE3dXV1Y7DZIJVN6te4G6h0k63vXth7fIzS0tReRp0RIIB7gQuUdVLgRngRgBVHVPVkqqWgLcAFVXdaEzf++ptVfWOdgQ9NDTUjsNkglU3q17gbqHSTrfHHmutPCkdkWBU9XOqWn+c9X7g3HWavRn4q/ZFtTV+2h4eVr3A3UKlnW7nnddaeVI6IsGs4VrgM+uU/wSbJ5i3xpfY/lxEzkgntKczMDDQjsNkglU3q17gbqHSTrd9+6Cn5+llPT1ReRq07eKfiNwFnL1O1V5VvTVusxdYAcbWbPtCYElVD22w+z8B3glo/P57RIlqvTj2AHsgesBpfHyckZERKpUKS0tLjI6OMjExQV9fH93d3czNzVEsFpmZmWFlZYViscjk5CT9/f3Mzc0xNzdHqVSiXC7T1dXF0NAQ5XKZgYEBlpeXWVhYWN1nT08PhUKBqakpCoUCtVqNarW6Wp/L5ejv72d6eprBwUGq1SqLi4ur9b29veTzeWZnZxkeHmZ+fp5arbZan8/nyeVyVCqVE3aCaNjkd3/3d3PfffeZciqVSkxOTvK1r33NlFP9b2/Hjh3Mzc2Zcqr3k6qac6r306FDhzjrrLPa4lQs5tm372x++7dzfP3rz6avb5lrrz3MlVcOMD7evFPTqGpHvIBrgC8APevUvQ/4703upwAcaqbt6OioJuHuu+9OtH0nY9XNqpequ4VKiG7AAW3i/9iOGJohIq8GbgBepqpLa+p2AG8CXrrJ9ueo6hPx19cDG53pbCujo6PtOEwmWHWz6gXuFiqW3TrlHswfAjngzniY8fsb6l4KPK6qX2ncQERuEpFd8dd3i0hZRA4CLwfe1o6gJyYm2nGYTLDqZtUL3C1ULLt1xBmMql6wSd09wIvWKb+u4fNb0olsc3rW3i0zhFU3q17gbqFi2a1TzmCCpFAoZB1Calh1s+oF7hYqlt08wSRgamoq6xBSw6qbVS9wt1Cx7CbRgICTExH5OvBogl2cCXxjm8LpNKy6WfUCdwuVEN12qupzt2p0UieYpIjIAVXdtXXL8LDqZtUL3C1ULLv5JTLHcRwnFTzBOI7jOKngCSYZ+7MOIEWsuln1AncLFbNufg/GcRzHSQU/g3Ecx3FSwRPMJojIs0Xkn0TkiyLykIj8xjptThWRj4rIYRF5QEQK7Y+0dZp0+2kR+XrDSqHXrbevTkVEThGRB0Xk9nXqguy3Olu4BdtvIlKJp32aFJED69SLiPxB3G8HReT7sojzRGjC7QdF5JsN/fY/s4hzO+mIqWI6mCeBy1X1X0TkWcDnReQzqnp/Q5ufBY6q6gUicjXwLqK1azqdZtwAPqqqb80gvu3gF4GHgdPXqQu13+ps5gZh99vLVXWj50J+GLgwfr2QaKmOF7YrsG1gMzeAv1fVH21bNCnjZzCbEM9M/S/x12fFr7U3rV4H3Bx//jjwChGRNoV4wjTpFiwici7wI8BNGzQJst+gKTfLvA74YPz3ez/QKyLnZB2Usz6eYLYgvhQxCXwNuFNVH1jT5HuBOQCNln3+JpBvb5QnRhNuAFfFlyI+LiIhLSv4+8B/A57aoD7YfmNrNwi33xT4nIhMxIsDrmW132Iej8tCYCs3gBfHl60/IyIXtzO4NPAEswWq+m1VLQHnAi8QkUvWNFnvV28QZwJNuH0KKKjqpcBdfOcXf0cjIj8KfE1VN5sHPch+a9ItyH6LeYmqfh/RpbBfEJG160AF2W8xW7n9P6IpWC4D/g/wyXYHuN14gmkSVV0E7gFevabqcWAAQES6gO8BjrQ1uIRs5KaqVVV9Mv76Z0AoKyO9BHitiFSAjwCXi8iH17QJtd+2dAu431DV+fj9a8AtwAvWNFntt5hzgfn2RJeMrdxU9Vj9srWq3gE8S0TObHug24gnmE0QkeeKSG/8+buAVwJfWtPsNqLlngHeCNytATxc1IzbmmvbryW6qdzxqOqNqnquqhaAq4n65CfXNAuy35pxC7XfROQ0EcnVPwOv4pmr094G/FQ8muxFwDcbVrPtWJpxE5Gz6/cBReQFRP8/V9sd63bio8g25xzgZhE5haizP6aqt4vIbxKtSX0b8AHgQyJymOgX8NXZhdsSzbj9FxF5LbBC5PbTmUW7DRjpt3Ux0m99wC3x/7FdwF+q6t+KyPUAqvp+4A7gNcBhYAn4mYxibZVm3N4I/JyIrADfAq4O4UfPZviT/I7jOE4q+CUyx3EcJxU8wTiO4zip4AnGcRzHSQVPMI7jOE4qeIJxHMdxUsETjOM4jpMKnmAcx3GcVPAE4zgpIyKvFJEPneC23yUi98YPxG7UpltE7ounvHGcjsETjOOkz2XAF09w22uBv1HVb2/UQFWXgb8jrPVsnJMATzCOkz6XAZMi8vz4TOMhEbmrPpGhiFwUlx8UkV+Jp6+psxu4tf5FRG4Rkd8Skb8XkX8WkVfGVZ+M2zpOx+AJxnHS5zKgDHwC+EVVvRi4E3hbfFlrLC6/FHge8SSIItINPE9VKw37ugRYVNUfAH6e7ySVQ8C/b4OL4zSNJxjHSZF4OerTgR8EPq+qD8ZVU8BZwBuAL64pr19OOxNYbNhXD9GyAu+Li7rq9fEltOX6jL2O0wl4gnGcdBkhmi5/hOgspk6RKJlcCkw2lF/S8P1bwLMb6i4GJhrux1zK06d8PxU4vm2RO05CPME4TrpcRpQwvkqUZBCR5wFvAT5ItN7HUFxeAn6S+AxGVY8Cp4hIPck0Jh+IEszBeNs88HVV/beUfRynaTzBOE661EeQfQjoF5Ey0UqU16pqNS7fJSL/l2jEWEVVv9Kw/eeA748/F3nm2U79DOblRGulOE7H4OvBOE6GiMh315fJFZFfAb5HVX+tof7fAb+kqm/ZYj9/A9yoqtOpBuw4LeBnMI6TLW+Lhy1PAgXgnY2V8c3/8a0etAQ+6cnF6TT8DMZxHMdJBT+DcRzHcVLBE4zjOI6TCp5gHMdxnFTwBOM4juOkgicYx3EcJxU8wTiO4zip4AnGcRzHSQVPMI7jOE4q/H+AVwfhpVHPZwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"slopes: [-4.000515590129417, -4.000127253447741, -4.000058127895366, -4.000547200856506]\n"
]
}
],
"source": [
"n_vals=[20,40,80,160,320] # a list of sample n values for using in Simpon's rule\n",
"Err_vals = [abs(Ival-simpsons_rule(f,a,b,n)) for n in n_vals] # list of errors: |Ix - S_n| for each n\n",
"\n",
"#we'll also look at the errors in the midpoint rule\n",
"Err_vals_mid = [abs(Ival-midpoint_rule(f,a,b,n)) for n in n_vals]\n",
"\n",
"# building a plot\n",
"plt.plot(np.log(n_vals),np.log(Err_vals),'bo') # blue points are from Simpson's rule\n",
"plt.plot(np.log(n_vals),np.log(Err_vals_mid),'ro') # red points are from midpoint rule\n",
"plt.xlabel('$log(n)$') # label for x-axis\n",
"plt.ylabel('$log|I-approx_n|$') # label for y-axis\n",
"plt.grid(color='tab:gray', linestyle='--', linewidth=0.5)\n",
"plt.show()\n",
"\n",
"# blue point look like they live on a line, lets look at slopes between consecutive points\n",
"n_vals_l = np.log(n_vals) \n",
"Err_vals_l = np.log(Err_vals)\n",
"slopes=[(Err_vals_l[i+1]-Err_vals_l[i])/(n_vals_l[i+1]-n_vals_l[i]) for i in range(len(n_vals)-1)]\n",
"print(\"slopes: \", slopes)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The slopes between consecutive points using Simpson's rule is $-4$. The illustrates that the rate at which is $S_n$ is converging to $I_x$ is of order $n^{-4}$, which is what the errors bound would suggest."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}