## Numeric Interpolation, Differentiation and Integratio

### Interpolation

# Suppose you know that a function takes the values given in the following table
# 
#     x                y=f(x)
#     0                 0
#     1                 1
#     2                 8
#     3                27
#     4                64
#     5               125
#     
# *Problem :* Find a polynomial $p$ such that $f(x) = p(x)$ for the given values of $x$.
# 
# There are infinitely many such polynomials. One such polynomial can be constructed using
# Langrange's interpolation formula. Suppose we are given $n$ points $x_1, \dotsc, x_n$ and
# the value of $f$ at these points $y_1, \dotsc, y_n$.
# 
# Define
# \begin{equation}
# l_i(x) = \frac{\prod_{j \neq i}(x - x_j)}{\prod_{j \neq i}(x_i - x_j)}
# \end{equation}
# These are definitely polynomials which take the value $1$ at $x_i$ and $0$ at all other $x_j$'s.
# We define
# 
# $p(x) = \sum_{i=1}^n y_i l_i(x)$ to be the required polynomail.
# 
# This is called the Lagrange interpolation formula.

# In[1]:

# We create a class called interpolation.

class LagrangeInterpolation :
    """Creates a class for lagrange interpolation. Creates a callable polynomial.
    
    Attributes : 
        l : a list of (x, f(x)) pairs.
        
    Special methods define :
        __init__
        __call__
        __str__
        
    Other methods :
        add_pair(p) : Add pair to the list l
        poly() : Returns a list containing the coefficients of the polynomial.
    """
    
    def __init__(self, l) :
        self.l = l
    
    def add_pair(self, p) :
        self.l.append(p)
    
    def poly_add(self, l1, l2) :
        """Given two polynomial coeff lists, it adds them."""
        len1 = len(l1)
        len2 = len(l2)
        maxl = max(len1, len2)
        rlist = []
        for i in range(maxl) :
            if i < len1 and i < len2 :
                rlist.append(l1[i] + l2[i])
            elif i < len1 :
                rlist.append(l1[i])
            elif i < len2 :
                rlist.append(l2[i])
            else :
                print "Serious error. Execution should not reach here."
                rlist = None
        return rlist
    
    def poly_mult(self, l1, l2) :
        """Product of two poly coeff lists."""
        len1 = len(l1)
        len2 = len(l2)
        prodlen = len1 + len2 - 1 # Why?
        #print len1, len2, prodlen
        
        prod= []
        for i in range(prodlen) :
            i_th_coeff = 0
            for j in range(i+1) :
                #print i, j
                if j < len1 and i - j < len2 and j >= 0 and i >= j:
                    #print j, i-j, len1, len2
                    i_th_coeff += l1[j] * l2[i - j]
            prod.append(i_th_coeff)
        return prod
    
    def li(self, i) :
        l = self.l
        lenlist = len(l)
        if i >= lenlist :
            print "i = %d is out of range." % i
            li = None
        else :
            li = [1.0]
            xi = float(l[i][0])
            for j in range(lenlist) :
                if j != i : 
                    xj = float(l[j][0])
                    #print i, j, xi, xj, li
                    li = self.poly_mult(li, [-(xj/(xi - xj)), (1.0/(xi - xj))])
        return li
    
    def poly(self) :
        l = self.l
        lenl = len(l)
        p = [0.0]
        for i in range(lenl) :
            yi = [float(l[i][1])]
            ith_term = self.poly_mult(yi, self.li(i))
            p = self.poly_add(p, ith_term)
            #print ith_term, p
        return p
    
    def __call__(self, t) :
        l = self.l
        ll = len(l)
        valatx = 0
        x = []
        y = []
        for k in range(ll) :
            x.append(float(l[k][0]))
            y.append(float(l[k][1]))
        for i in range(ll) :
            ithterm = y[i]
            for j in range(ll) :
                if j != i :
                    ithterm *= (t - x[j]) / (x[i] - x[j])
            valatx += ithterm
        return valatx
    
    def __str__(self) :
        pl = self.poly()
        lpl = len(pl)
        strval = ""
        for i in range(lpl) :
            if i == 0 :
                strval += "%g" % pl[i]
            else :
                strval += " + %g x^%d" % (pl[i], i)
        return strval
    
    def __repr__(self) :
        return "LangrangeInterpolation(" + str(self.l) + ")"
            


# In[2]:

if __name__ == '__main__' :
    cl = LagrangeInterpolation([(1, 0), (2, 15), (3, 80)])
    print "Testing Interpolation :"
    print "Polynomial : "
    print cl.poly()
    cl.add_pair((-1, 0))
    cl.add_pair((10, 9999))
    print "__________________________"
    print "Polynomial after adding two points."
    print cl.poly()
    print "Value of the polynomial at .01 :", cl(.01)
    print "Testing __str__ :",
    print cl
    print "Testing __repr__ :",
    print cl.__repr__()


# Out[2]:

#     Testing Interpolation :
#     Polynomial : 
#     [35.0, -60.0, 25.0]
#     __________________________
#     Polynomial after adding two points.
#     [-0.9999999999999982, -5.329070518200751e-15, -3.552713678800501e-15, 3.552713678800501e-15, 1.0000000000000002]
#     Value of the polynomial at .01 : -0.99999999
#     Testing __str__ : -1 + -5.32907e-15 x^1 + -3.55271e-15 x^2 + 3.55271e-15 x^3 + 1 x^4
#     Testing __repr__ : LangrangeInterpolation([(1, 0), (2, 15), (3, 80), (-1, 0), (10, 9999)])
# 

### Differentiation

# We shall be differentiating with help of the approximation
# 
# \begin{equation}
# f'(x) \simeq \frac{f(x + h) - f(x - h)}{2h}
# \end{equation}
# 
# We implement that too as a class.

# In[3]:

class Derivative2 :
    """New class for derivatives
    
    Attributes : f, h, err, max_iter, show_steps
    Methods : seth, seterr, setiter, switch_show_steps
    """
    def __init__ (self, f) :
        self.f = f
        self.h = 1
        self.err = 1E-5
        self.max_iter = 10000
        self.show_steps = False
        
    def seth(self, h) :
        self.h = h
        return None
        
    def seterr(self, e) :
        self.err = e
        return None
    
    def setiter(self, N) :
        self.max_iter = N
        return None
    
    def switch_show_steps(self) :
        if self.show_steps :
            self.show_steps = False
        else :
            self.show_steps = True
        return None
            
    def _basic_der(self, y, d) :
        f = self.f
        return (f(y + d) - f(y - d))/(2*d)
            
    def __call__(self, x) :
        d = self.h
        f = self.f
        e = self.err
        
        df = self._basic_der(x, d)
        cont_loop = True
        no_iter = 0
        while cont_loop :
            no_iter += 1
            d /= 2.0
            dfnew = self._basic_der(x, d)
            if abs(df - dfnew) < e :
                retval = dfnew
                cont_loop = False
            elif no_iter > self.max_iter :
                print "Exceeded the maximum number of iterations."
                retval = None
                cont_loop = False
            else :
                df = dfnew
                if self.show_steps :
                    print "%06d. df = %10g" % (no_iter, df)
        return retval
                
    def __str__(self) :
        return "Cannot display the deriviative."


# In[4]:

if __name__ == '__main__' :
    der = Derivative2(lambda x : (1.0/3)*x*x*x + x)
    print "der(0) =", der(0), ", der(1) =", der(1)
    der.switch_show_steps()
    print "der(10) with steps:"
    print der(10)
    der.seterr(1E-50)
    der.seth(1E8)
    der.setiter(10)
    der.switch_show_steps()
    print "der(1.5) with err = 1E-50, h = 1E8, max_iter = 10 :", der(1.5)
    der.switch_show_steps()
    print "Same thing with steps :"
    print der(1.5)


# Out[4]:

#     der(0) = 1.00000127157 , der(1) = 2.00000127157
#     der(10) with steps:
#     000001. df =    101.083
#     000002. df =    101.021
#     000003. df =    101.005
#     000004. df =    101.001
#     000005. df =        101
#     000006. df =        101
#     000007. df =        101
#     000008. df =        101
#     101.000001272
#     der(1.5) with err = 1E-50, h = 1E8, max_iter = 10 : Exceeded the maximum number of iterations.
#     None
#     Same thing with steps :
#     000001. df = 8.33333e+14
#     000002. df = 2.08333e+14
#     000003. df = 5.20833e+13
#     000004. df = 1.30208e+13
#     000005. df = 3.25521e+12
#     000006. df = 8.13802e+11
#     000007. df = 2.03451e+11
#     000008. df = 5.08626e+10
#     000009. df = 1.27157e+10
#     000010. df = 3.17891e+09
#     Exceeded the maximum number of iterations.
#     None
# 

### Integration

# 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.
# 
# 
# Thus we get the equation
# 
# \begin{equation}
# \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)
# \end{equation}
# 
# The right hand side can be rewritten as
# \begin{equation}
# 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})
# \end{equation}

# #### Saving computations.
# Note that if we call the expression on the right hand side to be $S_N$, we have (check!) :
# 
# \begin{equation}
# 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}
# \end{equation}
# 
# We shall vary n starting from n = 1 and stop when two consecutiive $S_{2^n}$'s are close enough.

# In[5]:

class Integrate :
    """Given a function, and one end point it computes the integral to a point x"
    Attributes :
        f : function
        a : starting point of integration
        b : end point of integration
        max_iter : maximum number of iterations
        err : Error below which everything is considered zero.
        showstep : show steps
        
    Methods :
        integral(a, b) : Integral of f from a to b
        seta : sets a
        setb : Sets b
        setiter : sets max_iter
        seterr : sets error
        switch_showstep : Switches showstep
        __call__(x) : returns the integral from a preset a to x
    """
    
    def __init__(self, f, a) :
        self.f = f
        self.a = a
        self.b = 0
        self.max_iter = 10000
        self.err = 1E-5
        self.showstep = False
        
    def seta(self, val) :
        self.a = val
        return None
    
    def setb(self, val) :
        self.b = val
        return None
    
    def setiter(self, N) :
        self.max_iter = N
        return None
    
    def seterr(self, e) :
        self.err = e
        return None
    
    def switch_showstep(self) :
        if self.showstep :
            self.showstep = False
        else :
            self.showstep = True
    
    def integral(self, p, q) :
        f = self.f
        n = 1
        Sn = (q-p)*(f(p) + f(q))/4.0 + (q-p)*f((p+q)/2.0)/2.0
        halfpowern = 0.5
        twopowernm1 = 1
        cont_loop = True
        while cont_loop :
            n += 1
            halfpowern /= 2.0
            twopowernm1 *= 2
            Snm1 = Sn # S_{n-1}
            Tn = 0
            for k in range(twopowernm1) :
                Tn += f(p + (2*k + 1) * (q - p) * halfpowern)
            Sn = Snm1/2.0 + halfpowern * Tn * (q-p)
            
            if abs(Sn - Snm1) < self.err :
                retval = Sn
                cont_loop = False
            elif n >= self.max_iter :
                print "Maximum number of iterations exceeded."
                retval = None
                cont_loop = False
            else  :
                if self.showstep :
                    print "%07d. N = %20d, S_N = %25g." % (n, twopowernm1 * 2, Sn)
        return retval
    
    def __call__(self, x) :
        self.setb(x)
        return self.integral(self.a, self.b)
    
    def __str__(self) :
        return "Cannot display the integral."


# In[6]:

if __name__ == '__main__' :
    i = Integrate(lambda x : x*x, 0)
    i.switch_showstep()
    i.seterr(1E-10)
    i.setiter(10)
    print "Integral of x^2 between 0 and 100 with steps :"
    print i.integral(0, 100)
    i.setiter(10000)
    i.seterr(1E-5)
    print "i(3) ="
    print i(3)
    print "i =", i
    i.switch_showstep()
    print "i(6) =", i(6)


# Out[6]:

#     Integral of x^2 between 0 and 100 with steps :
#     0000002. N =                    4, S_N =                    343750.
#     0000003. N =                    8, S_N =                    335938.
#     0000004. N =                   16, S_N =                    333984.
#     0000005. N =                   32, S_N =                    333496.
#     0000006. N =                   64, S_N =                    333374.
#     0000007. N =                  128, S_N =                    333344.
#     0000008. N =                  256, S_N =                    333336.
#     0000009. N =                  512, S_N =                    333334.
#     Maximum number of iterations exceeded.
#     None
#     i(3) =
#     0000002. N =                    4, S_N =                   9.28125.
#     0000003. N =                    8, S_N =                   9.07031.
#     0000004. N =                   16, S_N =                   9.01758.
#     0000005. N =                   32, S_N =                   9.00439.
#     0000006. N =                   64, S_N =                    9.0011.
#     0000007. N =                  128, S_N =                   9.00027.
#     0000008. N =                  256, S_N =                   9.00007.
#     0000009. N =                  512, S_N =                   9.00002.
#     0000010. N =                 1024, S_N =                         9.
#     9.00000107288
#     i = Cannot display the integral.
#     i(6) = 72.0000021458
#