from sympy import Function,S,symbols,And,Piecewise,lambdify,Matrix,pi,latex
from sympy.parsing.sympy_parser import parse_expr
from sympy.functions.special.polynomials import jacobi
from numpy import array,float32,linspace
import src.lagrange as lg
import matplotlib.pyplot as plt
from src.dbg import dbg
'''
p=Piecewise((f,x<=0),(g,And(x>0,x<1)),(0,True))
p.subs(x,0)
p.integrate()
p.diff()
h=f*p
h.integrate()
s=lambdify(x,h,"numpy")
t=p.integrate((x,0,1))
t.evalf()
'''
class LeastSquares():
'''
find aproximation use least squares method
func must be set while initialisation
basis and points added directly or
use pointsSet for point
and createBasisH or createBasisJ for Head or Jacoby basis
see tests for examples
used as
l=LeastSquares(func='cos(x)')
l.points=list of points or
l.points=lg.choosingPoints(x0=0,n=10)
l.createBasisH()
l.leastSquares(weightFunction=1)
l.plot()
'''
def __init__(self,func='sin(x)'):
self.x=symbols('x')
self.points=lg.choosingPoints(x0=0,n=10)
self.basis=[]
self.texBasis=[]
self.texF=0
self.f=func
self.fLam=0
self.iF=0
def forEps(self,n):
x=self.x
W=(1-x)*(1+x)
self.points=[-1,1]
self.createBasisJ(n=n,alpha=1,beta=1)
self.leastSquares(weightFunction=W)
def eps(self):
a=linspace(self.points[0],self.points[-1],30)
ff=array([self.fLam(p) for p in a])
iFF=array([self.iF(p) for p in a])
return(max(ff-iFF))
def showInNotebook(self):
self.texBasis=[latex(b) for b in self.basis]
self.plot()
def test1H(self):
'''
find interplated function iF for function sin(x) (if you create LeastSquares object
as default)
in Chebyshev points
used Head function basis
'''
self.points=lg.choosingPoints(x0=0,n=10)
self.createBasisH()
self.leastSquares(weightFunction=1)
self.plot()
def test2H(self):
'''
find interplated function iF for function sin(x) (if you create LeastSquares object
as default)
in points from interval 0,\pi
used Head function basis
'''
self.points=[0,float(pi/6.),float(pi/4.),float(pi/3.),float(pi/2.),float(2*pi/3.),float(3*pi/4.),float(5*pi/6.),float(pi)]
self.createBasisH()
self.leastSquares(weightFunction=1)
self.plot()
def test1J(self):
'''
find interplated function iF for function sin(x) (if you create LeastSquares object
as default)
in Chebyshev points
used Jacoby basis
'''
x=self.x
W=(1-x)*(1+x)
self.points=[-1,1]
self.createBasisJ(n=5,alpha=1,beta=1)
self.leastSquares(weightFunction=W)
self.plot()
def plot(self):
x=self.x
points=array(self.points)
f=parse_expr(self.f)
f1=lambdify(x,f)
f=lambdify(x,f,"numpy")
a=linspace(points[0],points[-1],30)
ff=array([f1(p) for p in a])
iFF=array([self.iF(p) for p in a])
plt.plot(a,ff,a,iFF,'o')
plt.legend(('f(x)','interp f(x)'),loc='upper left')
plt.show()
def pointsSet(self,points):
self.points=points
def createBasisJ(self,n=5,alpha=0,beta=0):
'''
DESCRIPTION:
create basis from Jacobi polynomials
'''
x=self.x
basis=[jacobi(i,alpha,beta,x) for i in range(n)]
self.basis=basis
return(basis)
def createBasisH(self):
'''
DESCRIPTION:
create Head functions for points:
(x-p_{i-1})/(p_{i}-p_{i-1}) ,x \in [p_{i-1},p_{i}]
H_{i}(x)=(p_{i+1}-x)/(p_{i+1}-p_{i}) ,x \in [p_{i},p_{i+1}]
0,otherwise
it use sympy.Piecewise for conditions
every conditions must be rigorus (<,>), otherwise uninfinite integrals
will be wrong
'''
points=array(self.points)
H=[]
x=self.x
deltaPoints=(points[1:]-points[:-1])[:-1]
ddeltaPoints=points[2:]-points[1:-1]
for i in range(len(points))[1:-1]:
cond=[]
cond.append(((x-points[i-1])/(points[i]-points[i-1]),And(x>points[i-1],x<points[i])))
cond.append(((points[i+1]-x)/(points[i+1]-points[i]),And(x>points[i],x<points[i+1])))
cond.append((0,True))
H.append(Piecewise(*cond))
self.basis=H
return(H)
def leastSquares(self,weightFunction=1):
'''
DECRIPTIONS:
find a_{i} that satisfy
\integral(g(x)-\sum a_{j}*H_{j}(x))^2)->min
i.e:
find best aproximation of function func in form
L_{2} g= \sum a_{j}*H_{j}(x)
INPUT:
weight function:
for Head functions =1
for Jacobi functions=(1-x)**alpha*(1+x)**beta
self.points:
for Head functions same as for createBasisH
for Jacobi functions points=[-1,1]
self.func:
is function (as string) whom will be aproximated
when you useing Jacoby func mast be from -1 to 1
(for that use subset: f(x)=F((a+b)/2+(b-a)*x/2))
EXAMPLES:
forHead:
see self.test1H and test2H methods
l=LeastSquares(func='cos(x)')
l.points=lg.choocingPoints(x0=0,n=10)
l.createBasisH()
l.leastSquares()
l.plot()
for Jacoby:
see self.test1J methods
x=symbols('x')
W=(1-x)*(1+x)
l.points=[-1,1]
l.createBasisJ(n=5,alpha=1,beta=1)
l.leastSquares(weightFunction=W)
'''
x=self.x
points=self.points
basis=self.basis
func=self.f
f=parse_expr(func)
self.fLam=lambdify(x,f)
dbg(**{"leastSquares|part1":(points[0],points[-1])})
W=weightFunction
HH=Matrix([[(W*Hi*Hj).integrate((x,points[0],points[-1])) for Hj in basis] for Hi in basis])
b=Matrix([(W*f*Hi).integrate((x,points[0],points[-1])) for Hi in basis])
a=(HH**-1)*b
expr=sum([a[i]*basis[i] for i in range(len(basis))])
dbg(**{"leastSquares|part2":(HH,b,a,expr)})
self.texF=latex(expr)
expr1=lambdify(x,expr)
self.iF=expr1
return(expr)
class my_funct(Function):
nargs = 1
@classmethod
def eval(cls, x):
print(cls)
print(type(cls))
if x.is_Number:
if x is S.Zero:
return S.One
elif x is S.Infinity:
return S.Zero
