# Numerical Integration


```python
from itertools import combinations
```


```python
import sympy
from sympy import Function, integrate, Product, Sum, Symbol, symbols
from sympy.abc import a,b,h,i,k,m,n,x
from sympy import Rational as Rat
```


```python
def lagrange_basis_polys(N,x,xpts=None):
    """
    lagrange_basis_polynomials(N,x,xpts)
    returns the Lagrange basis polynomials as a list
    
    INPUTS/PARAMETERS
    -----------------
    <int> N - N > 0.  Note that there are N+1 points total
    <sympy.Symbol> x
    <list> xpts
    """
    assert N > 0
    if xpts != None:
        assert len(xpts) == N + 1

    if xpts == None:
        print "I'll generate symbolic sympy symbols for you for xpts"
        xpts = symbols('x0:'+str(N+1))
        
    basis_polys = []
    for i in range(N+1):
        tmpprod = Rat(1)
        for k in [k for k in range(N+1) if k != i]:
            tmpprod = tmpprod * (x - xpts[k])/(xpts[i]-xpts[k])
        basis_polys.append(tmpprod)
    return basis_polys

def lagrange_interp(N,x,xpts=None,ypts=None):
    """
    lagrange_interp(N,x,xpts,ypts)
    Lagrange interpolation formula
    """
    if xpts != None and ypts != None:
        assert len(xpts) == len(ypts)
    
    if xpts == None:
        print "I'll generate symbolic sympy symbols for you for xpts"
        xpts = symbols('x0:'+str(N+1))
    if ypts == None:
        print "I'll generate symbolic sympy symbols for you for xpts"
        ypts = symbols('y0:'+str(N+1))
        
    basis = lagrange_basis_polys(N,x,xpts)
    p_N = sum( [ypts[i]*basis[i] for i in range(N+1)] )
    return p_N
```


```python
xpts = symbols('x0:'+str(1+1))
ypts = symbols('y0:'+str(1+1))
p_1x = lagrange_interp(1,x,xpts,ypts)
```

Below is, mathematically, $f_{-h} := p_1(x)$ with $(x_0,y_0) = (x_0-h, f(x_0-h)), (x_1,y_1) = (x_0,f(x_0))$ and  
$\int_{x_0-h}^{x_0} f_{-h}$


```python
x_0 = Symbol('x_0',real=True)
f = Function('f')
f_minush = p_1x.subs({xpts[0]:x_0-h,xpts[1]:x_0, ypts[0]:f(x_0-h), ypts[1]:f(x_0) })
integrate( f_minush, (x,x_0-h,x_0 ) )
```




    -x_0**2*(-f(x_0) + f(-h + x_0))/(2*h) + x_0*(h*f(x_0) - x_0*f(x_0) + x_0*f(-h + x_0))/h + (-h + x_0)**2*(-f(x_0) + f(-h + x_0))/(2*h) - (-h + x_0)*(h*f(x_0) - x_0*f(x_0) + x_0*f(-h + x_0))/h



Then, we can use `sympy` to calculate, symbolically, $f_{h} := p_1(x)$ with $(x_0,y_0) = (x_0, f(x_0)), (x_1,y_1) = (x_0+h,f(x_0+h))$ and  
$\int_{x_0}^{x_0+h} f_{h}$


```python
f_h = p_1x.subs({xpts[0]:x_0,xpts[1]:x_0+h, ypts[0]:f(x_0), ypts[1]:f(x_0+h) })
integrate( f_h, (x,x_0,x_0+h ) )
```




    -x_0**2*(-f(x_0) + f(h + x_0))/(2*h) - x_0*(h*f(x_0) + x_0*f(x_0) - x_0*f(h + x_0))/h + (h + x_0)**2*(-f(x_0) + f(h + x_0))/(2*h) + (h + x_0)*(h*f(x_0) + x_0*f(x_0) - x_0*f(h + x_0))/h




```python
( integrate( f_minush, (x,x_0-h,x_0 ) ) + integrate( f_h, (x,x_0,x_0+h ) )  ).simplify()
```




    h*(2*f(x_0) + f(-h + x_0) + f(h + x_0))/2



Success! Trapezoid rule was rederived (stop using pen/pencil and paper or chalkboard; computers can do computations faster and without mistakes)

For a second order polynomial, $p_{N=2}(x)$,


```python
xpts = symbols('x0:'+str(2+1))
ypts = symbols('y0:'+str(2+1))
p_2x = lagrange_interp(2,x,xpts,ypts)
```


```python
f2_h = p_2x.subs({xpts[0]:x_0-h,xpts[1]:x_0,xpts[2]:x_0+h,ypts[0]:f(x_0-h), ypts[1]:f(x_0),ypts[2]:f(x_0+h) })
```


```python
integrate( f2_h,(x,x_0-h,x_0+h)).simplify()
```




    h*(4*f(x_0) + f(-h + x_0) + f(h + x_0))/3



# Legendre Polynomials

I don't find the `sympy` documentation very satisfying (other than listing the argument syntax, no examples of usage, nor further explanation, beyond the barebones argument syntax, is given).  So what I've done here is to try to show what I've done.   


```python
from sympy.polys.orthopolys import legendre_poly
```


```python
print "n \t \t \t \t P_n(x) \n"
for i in range(11):
    print str(i) + "\t \t \t \t " , legendre_poly(i,x) 
```

    n 	 	 	 	 P_n(x) 
    
    0	 	 	 	  1
    1	 	 	 	  x
    2	 	 	 	  3*x**2/2 - 1/2
    3	 	 	 	  5*x**3/2 - 3*x/2
    4	 	 	 	  35*x**4/8 - 15*x**2/4 + 3/8
    5	 	 	 	  63*x**5/8 - 35*x**3/4 + 15*x/8
    6	 	 	 	  231*x**6/16 - 315*x**4/16 + 105*x**2/16 - 5/16
    7	 	 	 	  429*x**7/16 - 693*x**5/16 + 315*x**3/16 - 35*x/16
    8	 	 	 	  6435*x**8/128 - 3003*x**6/32 + 3465*x**4/64 - 315*x**2/32 + 35/128
    9	 	 	 	  12155*x**9/128 - 6435*x**7/32 + 9009*x**5/64 - 1155*x**3/32 + 315*x/128
    10	 	 	 	  46189*x**10/256 - 109395*x**8/256 + 45045*x**6/128 - 15015*x**4/128 + 3465*x**2/256 - 63/256



```python
sympy.latex(legendre_poly(2,x))
```




    '\\frac{3 x^{2}}{2} - \\frac{1}{2}'




```python
sympy.N( sympy.integrate(1/(2+x**2),(x,0,3)) )
```




    0.799232657543987




```python

```

d7c1a43ecc906ac987b682ffad272ef8171aa077
ipynb
Jupyter Notebook
Cpp/Integrate.ipynb
ernestyalumni/CompPhys
1f5d7559146a14a21182653b77fd35e6d6829855
Apache-2.0
2017-07-24T04:09:27.000Z
2021-12-24T16:00:41.000Z
Cpp/Integrate.ipynb
ernestyalumni/CompPhys
1f5d7559146a14a21182653b77fd35e6d6829855
Apache-2.0
2018-01-16T22:34:47.000Z
2019-01-29T22:37:10.000Z
Cpp/Integrate.ipynb
ernestyalumni/CompPhys
1f5d7559146a14a21182653b77fd35e6d6829855
Apache-2.0
2017-01-24T19:18:42.000Z
2021-03-01T07:13:35.000Z
Qwen/Qwen-72B
1. YES
2. YES