{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Numeric Interpolation, Differentiation and Integratio" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose you know that a function takes the values given in the following table\n", "\n", " x y=f(x)\n", " 0 0\n", " 1 1\n", " 2 8\n", " 3 27\n", " 4 64\n", " 5 125\n", " \n", "*Problem :* Find a polynomial $p$ such that $f(x) = p(x)$ for the given values of $x$.\n", "\n", "There are infinitely many such polynomials. One such polynomial can be constructed using\n", "Langrange's interpolation formula. Suppose we are given $n$ points $x_1, \\dotsc, x_n$ and\n", "the value of $f$ at these points $y_1, \\dotsc, y_n$.\n", "\n", "Define\n", "\\begin{equation}\n", "l_i(x) = \\frac{\\prod_{j \\neq i}(x - x_j)}{\\prod_{j \\neq i}(x_i - x_j)}\n", "\\end{equation}\n", "These are definitely polynomials which take the value $1$ at $x_i$ and $0$ at all other $x_j$'s.\n", "We define\n", "\n", "$p(x) = \\sum_{i=1}^n y_i l_i(x)$ to be the required polynomail.\n", "\n", "This is called the Lagrange interpolation formula." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# We create a class called interpolation.\n", "\n", "class LagrangeInterpolation :\n", " \"\"\"Creates a class for lagrange interpolation. Creates a callable polynomial.\n", " \n", " Attributes : \n", " l : a list of (x, f(x)) pairs.\n", " \n", " Special methods define :\n", " __init__\n", " __call__\n", " __str__\n", " \n", " Other methods :\n", " add_pair(p) : Add pair to the list l\n", " poly() : Returns a list containing the coefficients of the polynomial.\n", " \"\"\"\n", " \n", " def __init__(self, l) :\n", " self.l = l\n", " \n", " def add_pair(self, p) :\n", " self.l.append(p)\n", " \n", " def poly_add(self, l1, l2) :\n", " \"\"\"Given two polynomial coeff lists, it adds them.\"\"\"\n", " len1 = len(l1)\n", " len2 = len(l2)\n", " maxl = max(len1, len2)\n", " rlist = []\n", " for i in range(maxl) :\n", " if i < len1 and i < len2 :\n", " rlist.append(l1[i] + l2[i])\n", " elif i < len1 :\n", " rlist.append(l1[i])\n", " elif i < len2 :\n", " rlist.append(l2[i])\n", " else :\n", " print \"Serious error. Execution should not reach here.\"\n", " rlist = None\n", " return rlist\n", " \n", " def poly_mult(self, l1, l2) :\n", " \"\"\"Product of two poly coeff lists.\"\"\"\n", " len1 = len(l1)\n", " len2 = len(l2)\n", " prodlen = len1 + len2 - 1 # Why?\n", " #print len1, len2, prodlen\n", " \n", " prod= []\n", " for i in range(prodlen) :\n", " i_th_coeff = 0\n", " for j in range(i+1) :\n", " #print i, j\n", " if j < len1 and i - j < len2 and j >= 0 and i >= j:\n", " #print j, i-j, len1, len2\n", " i_th_coeff += l1[j] * l2[i - j]\n", " prod.append(i_th_coeff)\n", " return prod\n", " \n", " def li(self, i) :\n", " l = self.l\n", " lenlist = len(l)\n", " if i >= lenlist :\n", " print \"i = %d is out of range.\" % i\n", " li = None\n", " else :\n", " li = [1.0]\n", " xi = float(l[i][0])\n", " for j in range(lenlist) :\n", " if j != i : \n", " xj = float(l[j][0])\n", " #print i, j, xi, xj, li\n", " li = self.poly_mult(li, [-(xj/(xi - xj)), (1.0/(xi - xj))])\n", " return li\n", " \n", " def poly(self) :\n", " l = self.l\n", " lenl = len(l)\n", " p = [0.0]\n", " for i in range(lenl) :\n", " yi = [float(l[i][1])]\n", " ith_term = self.poly_mult(yi, self.li(i))\n", " p = self.poly_add(p, ith_term)\n", " #print ith_term, p\n", " return p\n", " \n", " def __call__(self, t) :\n", " l = self.l\n", " ll = len(l)\n", " valatx = 0\n", " x = []\n", " y = []\n", " for k in range(ll) :\n", " x.append(float(l[k][0]))\n", " y.append(float(l[k][1]))\n", " for i in range(ll) :\n", " ithterm = y[i]\n", " for j in range(ll) :\n", " if j != i :\n", " ithterm *= (t - x[j]) / (x[i] - x[j])\n", " valatx += ithterm\n", " return valatx\n", " \n", " def __str__(self) :\n", " pl = self.poly()\n", " lpl = len(pl)\n", " strval = \"\"\n", " for i in range(lpl) :\n", " if i == 0 :\n", " strval += \"%g\" % pl[i]\n", " else :\n", " strval += \" + %g x^%d\" % (pl[i], i)\n", " return strval\n", " \n", " def __repr__(self) :\n", " return \"LangrangeInterpolation(\" + str(self.l) + \")\"\n", " " ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "if __name__ == '__main__' :\n", " cl = LagrangeInterpolation([(1, 0), (2, 15), (3, 80)])\n", " print \"Testing Interpolation :\"\n", " print \"Polynomial : \"\n", " print cl.poly()\n", " cl.add_pair((-1, 0))\n", " cl.add_pair((10, 9999))\n", " print \"__________________________\"\n", " print \"Polynomial after adding two points.\"\n", " print cl.poly()\n", " print \"Value of the polynomial at .01 :\", cl(.01)\n", " print \"Testing __str__ :\",\n", " print cl\n", " print \"Testing __repr__ :\",\n", " print cl.__repr__()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Testing Interpolation :\n", "Polynomial : \n", "[35.0, -60.0, 25.0]\n", "__________________________\n", "Polynomial after adding two points.\n", "[-0.9999999999999982, -5.329070518200751e-15, -3.552713678800501e-15, 3.552713678800501e-15, 1.0000000000000002]\n", "Value of the polynomial at .01 : -0.99999999\n", "Testing __str__ : -1 + -5.32907e-15 x^1 + -3.55271e-15 x^2 + 3.55271e-15 x^3 + 1 x^4\n", "Testing __repr__ : LangrangeInterpolation([(1, 0), (2, 15), (3, 80), (-1, 0), (10, 9999)])\n" ] } ], "prompt_number": 2 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Differentiation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall be differentiating with help of the approximation\n", "\n", "\\begin{equation}\n", "f'(x) \\simeq \\frac{f(x + h) - f(x - h)}{2h}\n", "\\end{equation}\n", "\n", "We implement that too as a class." ] }, { "cell_type": "code", "collapsed": false, "input": [ "class Derivative2 :\n", " \"\"\"New class for derivatives\n", " \n", " Attributes : f, h, err, max_iter, show_steps\n", " Methods : seth, seterr, setiter, switch_show_steps\n", " \"\"\"\n", " def __init__ (self, f) :\n", " self.f = f\n", " self.h = 1\n", " self.err = 1E-5\n", " self.max_iter = 10000\n", " self.show_steps = False\n", " \n", " def seth(self, h) :\n", " self.h = h\n", " return None\n", " \n", " def seterr(self, e) :\n", " self.err = e\n", " return None\n", " \n", " def setiter(self, N) :\n", " self.max_iter = N\n", " return None\n", " \n", " def switch_show_steps(self) :\n", " if self.show_steps :\n", " self.show_steps = False\n", " else :\n", " self.show_steps = True\n", " return None\n", " \n", " def _basic_der(self, y, d) :\n", " f = self.f\n", " return (f(y + d) - f(y - d))/(2*d)\n", " \n", " def __call__(self, x) :\n", " d = self.h\n", " f = self.f\n", " e = self.err\n", " \n", " df = self._basic_der(x, d)\n", " cont_loop = True\n", " no_iter = 0\n", " while cont_loop :\n", " no_iter += 1\n", " d /= 2.0\n", " dfnew = self._basic_der(x, d)\n", " if abs(df - dfnew) < e :\n", " retval = dfnew\n", " cont_loop = False\n", " elif no_iter > self.max_iter :\n", " print \"Exceeded the maximum number of iterations.\"\n", " retval = None\n", " cont_loop = False\n", " else :\n", " df = dfnew\n", " if self.show_steps :\n", " print \"%06d. df = %10g\" % (no_iter, df)\n", " return retval\n", " \n", " def __str__(self) :\n", " return \"Cannot display the deriviative.\"" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "if __name__ == '__main__' :\n", " der = Derivative2(lambda x : (1.0/3)*x*x*x + x)\n", " print \"der(0) =\", der(0), \", der(1) =\", der(1)\n", " der.switch_show_steps()\n", " print \"der(10) with steps:\"\n", " print der(10)\n", " der.seterr(1E-50)\n", " der.seth(1E8)\n", " der.setiter(10)\n", " der.switch_show_steps()\n", " print \"der(1.5) with err = 1E-50, h = 1E8, max_iter = 10 :\", der(1.5)\n", " der.switch_show_steps()\n", " print \"Same thing with steps :\"\n", " print der(1.5)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "der(0) = 1.00000127157 , der(1) = 2.00000127157\n", "der(10) with steps:\n", "000001. df = 101.083\n", "000002. df = 101.021\n", "000003. df = 101.005\n", "000004. df = 101.001\n", "000005. df = 101\n", "000006. df = 101\n", "000007. df = 101\n", "000008. df = 101\n", "101.000001272\n", "der(1.5) with err = 1E-50, h = 1E8, max_iter = 10 : Exceeded the maximum number of iterations.\n", "None\n", "Same thing with steps :\n", "000001. df = 8.33333e+14\n", "000002. df = 2.08333e+14\n", "000003. df = 5.20833e+13\n", "000004. df = 1.30208e+13\n", "000005. df = 3.25521e+12\n", "000006. df = 8.13802e+11\n", "000007. df = 2.03451e+11\n", "000008. df = 5.08626e+10\n", "000009. df = 1.27157e+10\n", "000010. df = 3.17891e+09\n", "Exceeded the maximum number of iterations.\n", "None\n" ] } ], "prompt_number": 4 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Integration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall use the *trapezoidal* method for integration. Suppose we want to integrate a function $f : \\mathbb{R} \\to \\mathbb{R}$ from $a$ to $b$. The method involves partitioning the interval $[a, b]$ into $N$ equal intervals $[a, a + (b-a)/N]$, $[a + (b-a)/N, a + 2(b-a)/N]$, ..., $[a + (N-1)(b-a)/N, b]$. Now the assumption is that when $N$ is sufficiently large, the intervals are so small that $f$ is practically linear. Thus to integrate $f$ on the interval $[l_i = a + (i-1)(b-a)/N, r_i = a + i(b-a)/N]$, one just has to find the integral of the line joining $(l_i, f(l_i))$ and $(r_i, f(r_i))$ for $i = 1, 2, \\dotsc, N$. The integral of each such line is $0.5(r_i - l_i)(f(r_i) + f(l_i))$. Summing these over all $i$ we get the integral.\n", "\n", "\n", "Thus we get the equation\n", "\n", "\\begin{equation}\n", "\\int_a^b f dx = \\sum_{i=1}^N \\int_{a + (i-1)\\frac{b-a}{N}}^{a + i\\frac{b-a}{N}} f dx \\simeq \\sum_{i=1}^N \\frac{b-a}{2N} \\left(f(a + (i-1)\\frac{b-a}{N}) + f(a + i\\frac{b-a}{N})\\right) = \\frac{b-a}{2N} \\sum_{i=1}^N \\left(f(a + (i-1)\\frac{b-a}{N}) + f(a + i\\frac{b-a}{N})\\right)\n", "\\end{equation}\n", "\n", "The right hand side can be rewritten as\n", "\\begin{equation}\n", "S_N := \\frac{b-a}{N} \\frac{f(a) + f(b)}{2} + \\frac{b-a}{N} \\sum_{i = 1}^{N-1} f(a + i \\frac{b-a}{N})\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Saving computations.\n", "Note that if we call the expression on the right hand side to be $S_N$, we have (check!) :\n", "\n", "\\begin{equation}\n", "S_{2^{n+1}} = \\frac{1}{2} S_{2^n} + \\frac{b-a}{2^{n+1}} \\sum_{i=0}^{2^n - 1} f(a + (2i + 1) \\frac{b-a}{2^{n+1}}) = \\frac{1}{2} S_{2^n} + \\frac{b-a}{2^{n+1}} T_{n+1}\n", "\\end{equation}\n", "\n", "We shall vary n starting from n = 1 and stop when two consecutiive $S_{2^n}$'s are close enough." ] }, { "cell_type": "code", "collapsed": false, "input": [ "class Integrate :\n", " \"\"\"Given a function, and one end point it computes the integral to a point x\"\n", " Attributes :\n", " f : function\n", " a : starting point of integration\n", " b : end point of integration\n", " max_iter : maximum number of iterations\n", " err : Error below which everything is considered zero.\n", " showstep : show steps\n", " \n", " Methods :\n", " integral(a, b) : Integral of f from a to b\n", " seta : sets a\n", " setb : Sets b\n", " setiter : sets max_iter\n", " seterr : sets error\n", " switch_showstep : Switches showstep\n", " __call__(x) : returns the integral from a preset a to x\n", " \"\"\"\n", " \n", " def __init__(self, f, a) :\n", " self.f = f\n", " self.a = a\n", " self.b = 0\n", " self.max_iter = 10000\n", " self.err = 1E-5\n", " self.showstep = False\n", " \n", " def seta(self, val) :\n", " self.a = val\n", " return None\n", " \n", " def setb(self, val) :\n", " self.b = val\n", " return None\n", " \n", " def setiter(self, N) :\n", " self.max_iter = N\n", " return None\n", " \n", " def seterr(self, e) :\n", " self.err = e\n", " return None\n", " \n", " def switch_showstep(self) :\n", " if self.showstep :\n", " self.showstep = False\n", " else :\n", " self.showstep = True\n", " \n", " def integral(self, p, q) :\n", " f = self.f\n", " n = 1\n", " Sn = (q-p)*(f(p) + f(q))/4.0 + (q-p)*f((p+q)/2.0)/2.0\n", " halfpowern = 0.5\n", " twopowernm1 = 1\n", " cont_loop = True\n", " while cont_loop :\n", " n += 1\n", " halfpowern /= 2.0\n", " twopowernm1 *= 2\n", " Snm1 = Sn # S_{n-1}\n", " Tn = 0\n", " for k in range(twopowernm1) :\n", " Tn += f(p + (2*k + 1) * (q - p) * halfpowern)\n", " Sn = Snm1/2.0 + halfpowern * Tn * (q-p)\n", " \n", " if abs(Sn - Snm1) < self.err :\n", " retval = Sn\n", " cont_loop = False\n", " elif n >= self.max_iter :\n", " print \"Maximum number of iterations exceeded.\"\n", " retval = None\n", " cont_loop = False\n", " else :\n", " if self.showstep :\n", " print \"%07d. N = %20d, S_N = %25g.\" % (n, twopowernm1 * 2, Sn)\n", " return retval\n", " \n", " def __call__(self, x) :\n", " self.setb(x)\n", " return self.integral(self.a, self.b)\n", " \n", " def __str__(self) :\n", " return \"Cannot display the integral.\"" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "if __name__ == '__main__' :\n", " i = Integrate(lambda x : x*x, 0)\n", " i.switch_showstep()\n", " i.seterr(1E-10)\n", " i.setiter(10)\n", " print \"Integral of x^2 between 0 and 100 with steps :\"\n", " print i.integral(0, 100)\n", " i.setiter(10000)\n", " i.seterr(1E-5)\n", " print \"i(3) =\"\n", " print i(3)\n", " print \"i =\", i\n", " i.switch_showstep()\n", " print \"i(6) =\", i(6)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Integral of x^2 between 0 and 100 with steps :\n", "0000002. N = 4, S_N = 343750.\n", "0000003. N = 8, S_N = 335938.\n", "0000004. N = 16, S_N = 333984.\n", "0000005. N = 32, S_N = 333496.\n", "0000006. N = 64, S_N = 333374.\n", "0000007. N = 128, S_N = 333344.\n", "0000008. N = 256, S_N = 333336.\n", "0000009. N = 512, S_N = 333334.\n", "Maximum number of iterations exceeded.\n", "None\n", "i(3) =\n", "0000002. N = 4, S_N = 9.28125.\n", "0000003. N = 8, S_N = 9.07031.\n", "0000004. N = 16, S_N = 9.01758.\n", "0000005. N = 32, S_N = 9.00439.\n", "0000006. N = 64, S_N = 9.0011.\n", "0000007. N = 128, S_N = 9.00027.\n", "0000008. N = 256, S_N = 9.00007.\n", "0000009. N = 512, S_N = 9.00002.\n", "0000010. N = 1024, S_N = 9.\n", "9.00000107288\n", "i = Cannot display the integral.\n", "i(6) = 72.0000021458\n" ] } ], "prompt_number": 6 } ], "metadata": {} } ] }