{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Trapezoidal and Simpson rule" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next function implements Trapezoid and Simpson methods." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "# Performs Trapezoid and Simpson quadrature\n", "def integrate(a,b,n,f):\n", " h = (b-a)/n\n", " x = np.linspace(a,b,n+1)\n", " y = f(x)\n", " res1 = 0.5*y[0] + np.sum(y[1:n]) + 0.5*y[n]\n", " res2 = 4.0*np.sum(y[1:n:2]) + 2.0*np.sum(y[2:n-1:2]) + y[0] + y[n]\n", " return h*res1, (h/3.0)*res2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next function performs convergence test." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "def test(a,b,f,Ie,n,N):\n", " e1,e2 = np.zeros(N), np.zeros(N)\n", " for i in range(N):\n", " I1,I2 = integrate(a,b,n,f)\n", " e1[i],e2[i] = Ie - I1, Ie - I2\n", " if i > 0:\n", " print('%6d %18.8e %14.5g %18.8e %14.5g'%\n", " (n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))\n", " else:\n", " print('%6d %18.8e %14.5g %18.8e %14.5g'%(n,e1[i],0,e2[i],0))\n", " n = 2*n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example\n", "$$\n", "f(x) = x^3 \\sqrt{x}, \\qquad x \\in [0,1]\n", "$$\n", "The exact integral is $2/9$." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2 -7.19719516e-02 0 -3.37000954e-03 0\n", " 4 -1.81666065e-02 3.9618 -2.31491460e-04 14.558\n", " 8 -4.55322411e-03 3.9898 -1.54299774e-05 15.003\n", " 16 -1.13906170e-03 3.9973 -1.00755946e-06 15.314\n", " 32 -2.84814093e-04 3.9993 -6.48919993e-08 15.527\n", " 64 -7.12066288e-05 3.9998 -4.14075488e-09 15.672\n", " 128 -1.78018541e-05 4 -2.62556588e-10 15.771\n", " 256 -4.45047596e-06 4 -1.65759628e-11 15.84\n", " 512 -1.11261977e-06 4 -1.04333209e-12 15.888\n", " 1024 -2.78154993e-07 4 -6.55309140e-14 15.921\n" ] } ], "source": [ "f = lambda x: x**3 * np.sqrt(x)\n", "a, b = 0.0, 1.0\n", "Ie = 2.0/9.0 # Exact integral\n", "n, N = 2, 10\n", "test(a,b,f,Ie,n,N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\frac{1}{1 + (x-\\pi)^2}, \\qquad x \\in [0,5]\n", "$$\n", "The exact integral is $\\arctan(5-\\pi) + \\arctan(\\pi)$." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2 1.73111440e-01 0 -2.85329178e-01 0\n", " 4 7.10986397e-02 2.4348 3.70943729e-02 -7.692\n", " 8 7.49581780e-03 9.4851 -1.37051228e-02 -2.7066\n", " 16 1.95340266e-03 3.8373 1.05930941e-04 -129.38\n", " 32 4.89160655e-04 3.9934 1.07998844e-06 98.085\n", " 64 1.22340738e-04 3.9983 6.74323792e-08 16.016\n", " 128 3.05883472e-05 3.9996 4.21691793e-09 15.991\n", " 256 7.64728450e-06 3.9999 2.63594035e-10 15.998\n", " 512 1.91183348e-06 4 1.64748215e-11 16\n", " 1024 4.77959142e-07 4 1.02939879e-12 16.004\n" ] } ], "source": [ "f = lambda x: 1.0/(1.0 + (x-np.pi)**2)\n", "a, b = 0.0, 5.0\n", "Ie = np.arctan(b-np.pi) + np.arctan(np.pi)\n", "n, N = 2, 10\n", "test(a,b,f,Ie,n,N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\sqrt{x}, \\qquad x \\in [0,1]\n", "$$\n", "The exact integral is $2/3$." ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2 6.31132761e-02 0 2.85954792e-02 0\n", " 4 2.33836204e-02 2.699 1.01404019e-02 2.82\n", " 8 8.53644504e-03 2.7393 3.58738658e-03 2.8267\n", " 16 3.08546979e-03 2.7667 1.26847804e-03 2.8281\n", " 32 1.10773039e-03 2.7854 4.48483920e-04 2.8284\n", " 64 3.95855288e-04 2.7983 1.58563588e-04 2.8284\n", " 128 1.41009370e-04 2.8073 5.60607304e-05 2.8284\n", " 256 5.01176901e-05 2.8136 1.98204636e-05 2.8284\n", " 512 1.77851167e-05 2.818 7.00759224e-06 2.8284\n", " 1024 6.30444768e-06 2.821 2.47755801e-06 2.8284\n" ] } ], "source": [ "f = lambda x: np.sqrt(x)\n", "a, b = 0.0, 1.0\n", "Ie = 2.0/3.0\n", "n, N = 2, 10\n", "test(a,b,f,Ie,n,N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $f'(0)$ is not finite, we do not get the optimal convergence rates. But the errors still decrease and both methods converge at same rate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x) = \\exp(\\cos x), \\qquad x \\in [0,2\\pi]\n", "$$\n", "The exact integral is 7.95492652101284." ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2 -1.74053505e+00 0 7.20800573e-01 0\n", " 4 -3.43969188e-02 50.601 5.34315792e-01 1.349\n", " 8 -1.25168894e-06 27480 1.14639707e-02 46.608\n", " 16 -3.55271368e-15 3.5232e+08 4.17229640e-07 27476\n", " 32 -3.55271368e-15 1 -1.77635684e-15 -2.3488e+08\n" ] } ], "source": [ "f = lambda x: np.exp(np.cos(x))\n", "a, b = 0.0, 2*np.pi\n", "Ie = 7.95492652101284\n", "n, N = 2, 5\n", "test(a,b,f,Ie,n,N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The integrand is periodic and the error formula suggests that we should expect fast convergence." ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }