In [4]:
from itertools import combinations
In [2]:
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
In [3]:
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
In [6]:
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}$
In [7]:
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 ) )
Out[7]:
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}$
In [61]:
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 ) )
Out[61]:
In [63]:
( integrate( f_minush, (x,x_0-h,x_0 ) ) + integrate( f_h, (x,x_0,x_0+h ) ) ).simplify()
Out[63]:
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)$,
In [8]:
xpts = symbols('x0:'+str(2+1))
ypts = symbols('y0:'+str(2+1))
p_2x = lagrange_interp(2,x,xpts,ypts)
In [10]:
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) })
In [12]:
integrate( f2_h,(x,x_0-h,x_0+h)).simplify()
Out[12]:
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.
In [1]:
from sympy.polys.orthopolys import legendre_poly
In [15]:
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)
In [13]:
sympy.latex(legendre_poly(2,x))
Out[13]:
In [6]:
sympy.N( sympy.integrate(1/(2+x**2),(x,0,3)) )
Out[6]:
In [ ]: