from numpy import array,ones,pi,cos,concatenate,linspace,sin,NaN,shape
from sympy import symbols,lambdify,latex,diff,Function,expand,sin
from sympy.utilities.lambdify import implemented_function
from sympy.parsing.sympy_parser import parse_expr
from src.dbg import dbg
import matplotlib.pyplot as plt
from pandas import DataFrame
class Lagrange():
def __init__(self,points=[-1,-1,1,1]):
self.points=points
self.x= symbols('x')
self.f=0
self.iF=0
self.A=0
self.d=[]
self.dd=[]
self.dataFrameD=0
self.dataFrameDD=0
self.x0=0
def showInNotebook(self):
if len(self.dd)!=0:
self.dataFrameDD=DataFrame(self.dd)
if len(self.d)!=0:
columns=[]
index=[]
m,n=shape(self.d)
for i in range(n):
if i==0:
columns.append(r'$'+latex("f(x_{i})")+'$')
elif i==1:
columns.append(r'$'+latex("f(x_{i};x_{j})")+'$')
elif i==n-1:
columns.append(r'$'+latex("f(x_{i};...;x_{j})")+'$')
else:
columns.append("...")
index=self.points
self.dataFrameD=DataFrame(self.d,index=index,columns=columns)
def plot(self):
x=self.x
iF=lambdify(x,self.iF,"numpy")
f=self.f
a=linspace(self.points[0],self.points[-1],30)
plt.plot(a,f(a),'o',a,iF(a))
plt.legend(('f(x)','interp f(x)'),loc='upper left')
plt.stem(array(self.points),iF(array(self.points)))
plt.show()
def lagrangeInterp(self,points=[],func="sin(pi*x)",y=0):
'''
INTUITION PHASE in solving Lagrange interpolation problem:
STEP 0:
a[1:]-a[:-1]
a[::2][1:]-a[::2][:-1]
a[1::2][1:]-a[1::2][:-1]
c[::3][:-1]=a[::3][1:]-a[::3][:-1]
c[1::3][:-1]=a[1::3][1:]-a[1::3][:-1]
c[2::3][:-1]=a[2::3][1:]-a[2::3][:-1]
c[3::3][:-1]=a[3::3][1:]-a[3::3][:-1]
c[4::3][:-1]=a[4::3][1:]-a[4::3][:-1]
c[5::3][:-1]=a[5::3][1:]-a[5::3][:-1]
c[6::3][:-1]=a[6::3][1:]-a[6::3][:-1]
6=9-3
for i in range(len(a)-1-3):
c[i::3]=a[i::3][1:]-a[i::3][:-1]
STEP 1:
c[::9]
c=d=ones((n+1,n))
c[0]=a
c[1]=f
QUESTIONs:
+c=ones((n+1),n)#i.e lines below
+c=ones((n+1),n)
+j from 1 or 2 ?
+if from 1 what menans i
STEP 5:
c[0]=a[1:]-a[:-1]
for j in range(len(a)-1+1)[2:]:
for i in range(len(a)-1-j+1):
c[j-1][i::j][:-1]=a[i::j][1:]-a[i::j][:-1]
d[j+2]=(d[j+1][1:]-d[j+1][:-1])/float(c[j+2])
SOME USEFUL REMARKs:
f = implemented_function(Function('f'), lambda x: x+1)
func = lambdify(x, f(x))
func(4)
expr=sp.parsing.sympy_parser.parse_expr("sin(x)+x**2")
sp.diff(expr)
ff=sp.lambdify(x,sp.diff(expr),"numpy")
TESTS:
import src.lagrange as lg
lg.lagrangeInterp(listX=[pi,2*pi,3*pi,2*pi])
lg.lagrangeInterp(listX=[0,-2,0.5,2/3.],func='x*sin(pi*x)')
lg.lagrangeInterp(listX=[-1,-1,1,1],func="sin(pi*x)")
'''
if len(points)==0:
points=self.points
else:
self.points=list(points)
a=array(points)
n=len(points)
c=ones((n,n))*NaN
d=ones((n,n))*NaN
c[0][:-1]=a[1:]-a[:-1]
x=self.x
expr=parse_expr(func)
f=lambdify(x,expr,"numpy")
self.f=f
df=lambdify(x,diff(expr),"numpy")
dbg(**{"lagrangeInterp|part1":expr})
d[0][:-1]=(f(a[1:])-f(a[:-1]))
firstElementF=f(a[0])
if y:
y=array(y)
d[0][:-1]=y[1:]-y[:-1]
firstElementF=y[0]
dbg(**{"lagrangeInterp|part2":(len(c[0][:-1]),a,d[0])})
c,d=self.ifEqPoints(a,c,d,j=1,df=df)
for j in range(len(a)-1+1)[2:]:
for i in range(j):
if(len(a[i::j])>1):
c[j-1][i::j][:-1]=a[i::j][1:]-a[i::j][:-1]
dbg(**{"lagrangeInterp|part3":(i,j,c[0][i::j][:-1],a[i::j],a[i::j][1:]-a[i::j][:-1])})
d[j-1][:-(j)]=(d[j-2][1:-(j-1)]-d[j-2][:-j])
dbg(**{"lagrangeInterp|part4":(j,d[j-1][:],d[j-2][1:-(j-1)],d[j-2][:-j])})
c,d=self.ifEqPoints(a,c,d,j=j,df=df)
dbg(**{"lagrangeInterp|part5":(a,c,d)})
self.d=concatenate((d[::-1],array([f(a)])),axis=0)[::-1][:-1].T
A=d.T[0][:]
X=a[:]
dbg(**{"lagrangeInterp|part6":(A,X)})
exprSum=firstElementF
exprMult=1
for ai in d.T[0][:-1]:
exprMult=exprMult*(x-X[0])
X=X[1:]
exprSum=exprSum+ai*exprMult
dbg(**{"lagrangeInterp|part7":(f(a[0]),exprSum)})
g=lambdify(x,exprSum,"numpy")
self.iF=exprSum
self.A=concatenate([array([firstElementF]),A[:-1]])
return(exprSum)
def ifEqPoints(self,a,c,d,j=1,df=0):
if(0 not in c[j-1][:-j]):
d[j-1][:-j]=d[j-1][:-j]/c[j-1][:-j]
dbg(**{"ifEqPoints|part1":(j,c[j-1][:],d[j-1][:])})
else:
for i in range(len(c[j-1][:-j])):
if(c[j-1][:-j][i]!=0):
d[j-1][:-j][i]=d[j-1][:-j][i]/c[j-1][:-j][i]
dbg(**{"ifEqPoints|part21":(j,i,d[j-1][:-j][i])})
else:
if(df==0):print("FATAL ERROR!!!")
d[j-1][:-j][i]=df(a[:-1][i])
dbg(**{"ifEqPoints|part22":(j,i,d[j-1][:-j][i])})
return(c,d)
def diffSocket(self,func="sin(x)",x0=0,points=[]):
'''
DECRCIPTION:
find len(points)-1 derivation of function func
step1 (find b):
b[n]=a[n]
for k=n-1,...,1 do
b[k]=a[k]+(x0-x[k])b[k+1]
b[1] is f'(x0)
step2:
a=b
and try agen step1
after k+1 steps
b[1]=f'(x0),...,b[k]=f^{k}(x0)
INPUT:
func-whom derivations we want to find
x0 is point to find derivations in
points will be generated in small area near x0 as Chebyshev points
questions:
+x0 from left in points?
+check array.copy's
'''
self.x0=x0
points=array(points)
if not len(points):
points=choosingPoints(x0=x0,n=5)
self.lagrangeInterp(points=points,func=func)
a=self.A
b=ones((len(a),len(a)))
dbg(**{"diffSocket|part1":(a,points)})
for i in range(len(a)):
b[i][-1]=a[-1]
for k in range(len(a))[1:]:
b[i][:-1][-k]=a[:-1][-k]+(x0-points[:-1][-k])*b[i][-(k)]
dbg(**{"diffSocket|part21":((len(b)-k-1,len(b)-1-k,len(b)-1-k,len(b)-k,b[i][:-1][-k]))})
a=b[i].copy()
points=concatenate([array([x0]),points[:-1]])
dbg(**{"diffSocket|part22":points})
deriv=[self.fact(i)*b[-1][i] for i in range(len(b[-1]))[1:]]
dbg(**{"diffSocket|part3":(b,deriv)})
self.dd=concatenate((b[::-1],self.d),axis=0)
return(deriv)
def fact(self,n):
if n==1:
return(1)
else:
return(n*self.fact(n-1))
def choosingPoints(x0=0,n=10):
'''Choose different points near x0
'''
x=linspace(x0+0.1,x0+2,n)
x=array(range(n))
a=x0+0.1
b=x0+2
x=array([(a+b-(a-b)*cos((2*j-1)*pi/float(2*n))/float(cos(pi/float(2*n))))/2. for j in range(n+1)[1:]])
if x[0]>x[-1]:
x=x[x.argsort()]
return(x)
def dividedDifferences(l,f=lambda x:x):
if len(l)==2:
return((f(l[1])-f(l[0]))/float(l[1]-l[0]))
else:
dD=dividedDifferences
return((dD(l[1:],f)-dD(l[:-1],f))/float(l[-1]-l[0]))
