procédure Papa FLAMMY relation de récurrence et diagonalisation de matrice 2X2 FIBONACCI
J'utilise la même procédure que Papa FLAMMY, mais pour des relations de récurrence différentes de celle de Papa FLAMMY. c'est ok pour p = 3, mais pour p = 2, il semble que l'algèbre soit trop compliquée pour la dernière étape, qui est le calcul de la limite
Lorsque j'essaie de calculer cette limite, j'ai reçu le message suivant: RuntimeError: ensuite, mon navigateur Web Edge semble très lent avec p = 2, le problème vient de p = 2 et des 4 dernières ligne de code. si vous voulez expérimenter le pb, changez `si p <> 2:`
par `si p <> 127`
recurrence relation et matrice diagonalisation matrice 4X4
c'est magnifique L'algebre lineaire !
je recopie aussi mon code ci dessous, mais attention il est surement pas optimal, je pense que l'on peut faire plus simple en utilisant toutes les possibilités de SAGEMATH.
(peut etre meme qu'il suffit de 20 lignes de code SageMath optimal!), si quelqu'un peut me montrer une solution plus simple avec Sagemath , je suis preneur. surtout si quelqu'un arrive a calculer la limite pour p=2.
j'ajoute ce lien sur sageCell qui est sans doute mieux pour voir le resultat du code:
SageCell Link
- Code: Tout sélectionner
p=3
Nmax=2+p^2
# Ian Stewart Formula, ok for p=x
def IanS(n,p):
if n<=2*p:
# first formula below count the words with all sequences of contiguous >=p bits
# also count if all zero no duplicates possible
# count all sequences in word n when ones sequences >=p or all zero
# example :
#0000 n=4 p=3 2+1+1 =(n-p+1)/2 *(n-p+2) +1 =4
#1110
#0111
#1111
return (n-p+1)/2 *(n-p+2) +1
else:
# Ian Stewart formula below
return 2*IanS(n-1,p) - IanS(n-2,p) + IanS(n-p-1,p)
sIanTab=[]
nTab=[]
for n in range (p,Nmax):
sIanTab.append(IanS(n,p))
nTab.append(n)
print(' Ian Stewart Calculation, all ones packets >= '+str(p)+' in the word n or all zeros')
show(" ones packets >= "+str(p)+" List : " ,sIanTab)
show(" word of length n List : ",nTab)
" below this recurrence relation is automaticaly build \
it expresses the number of sequences of >= to p '1' consecutive plus 1 (the sequence all the bits to zero) "
if p>2 :
zeroSeq='+ 0*a_{n-'+str(2)+'} '
for i in (3..(p-1)) :
zeroSeq=zeroSeq + '+ 0*a_{n-'+str(i)+'} '
show(zeroSeq)
else :
zeroSeq=''
recurrenceEq='1*a_{n+1} = 2*a_{n+0} -1*a_{n-1} '+ zeroSeq + '+1*a_{n-'+str(p)+'}'
# you can write below your own recurrence relation (one blank between each term)
#recurrenceEq='1*a_{n+1} = 2*a_{n+0} -1*a_{n-1} +0*a_{n-2} +1*a_{n-3}'
" example : below, in a number of n bits length.(all ones packets >= 3 in the word n or all zeros) "
#recurrenceEq='1*a_{n+1} = 2*a_{n+0} -1*a_{n-1} +0*a_{n-2} +1*a_{n-3}'
# example for Fibonaci sequence, uncomment the 2 lines below
#recurrenceEq='1*a_{n+1} = 1*a_{n+0} +1*a_{n-1} '
#sIanTab=[1,1]
show(recurrenceEq)
show(LatexExpr(r"\text{Your recurrence EQ : } "+recurrenceEq+" "))
nL=[]
import re
nL=re.findall("[-+]?[.]?[\d]+(?:,\d\d\d)*[\.]?\d*(?:[eE][-+]?\d+)?", recurrenceEq)
show("nL : ",nL)
eQfactors=[]
eQexponants=[]
for i in range(0,len(nL)) :
if (i % 2) == 0:
eQfactors.append(Integer(float(nL[i])))
else :
eQexponants.append(Integer(-float(nL[i])+1))
"reverse coeff list in order to put the progression factors in the same direction \
as the progression of the polynomial with the variable n :\
like the serie [n^0, n^1 , n^2 .....]"
eQfactors.reverse()
" the -1 multiplication is present only for displaying the rhs polynomial in usual lhs way"
" stop to -2 because it does not apply to a_n which is already with good lhs signe"
fnL=[]
mR=len(eQfactors)-1
for i in (0..mR-1):
eQfactors[i]= -eQfactors[i]
fnL.append(Integer(Integer(-eQfactors[i])))
show("fnL : ",fnL)
#eQfactors[len(eQfactors)-1]=eQfactors[len(eQfactors)-1]
show("Eq factors List : ",eQfactors)
show("Eq Exponants List : ",eQexponants)
#R = PolynomialRing(QQ, 'x')
print " caracteristic polynomial p(x)"
var('x')
x = PolynomialRing(SR, 'x').gen()
p_x=0
for i in (0..(len(eQfactors)-1)):
p_x=p_x+eQfactors[i]*x^(eQexponants[i])
show("p(x) : ",p_x)
p_xRootsTuple=p_x.roots()
p_xRootsL=[]
index=0
indexT=0
# expanding the roots tuples (if multiple roots)
for t in p_xRootsTuple:
for r in range(0,t[1]) :
p_xRootsL.append(t[0])
index+=1
show("root ",str(index)," : ",p_xRootsTuple[indexT][0])
#show("root simplified",str(index)," : ",(LatexExpr(latex((p_xRootsTuple[indexT][0])._sympy_().simplify()))))
indexT+=1
# Build the matrix M by replacing a_00 by its equation
startSeqL=(sIanTab[0:eQexponants[-1]])[::-1]
show("start Seq : ",startSeqL)
startSeqM=matrix(SR,startSeqL).transpose()
aL=[fnL[::-1]] # fnl reversed to get the most significant digit up
sr=([1]+[0]*(mR-1))
for r in range(0,mR-1) :
aL.append(sr)
sr=sr[-1:]+sr[:-1] # shift right
A=matrix(QQ,aL)
show("start Seq matrix: ",startSeqM)
show(" A : ",A)
var('lambd')
lambd = PolynomialRing(SR, 'lambd').gen()
IM = identity_matrix(SR, A.dimensions()[0]);
AIM=(A-lambd*IM)
lambdaListTuple=(AIM.determinant()).roots()
lambdaList=[]
for t in lambdaListTuple:
for r in range(0,t[1]) :
lambdaList.append(t[0])
for e in lambdaList:
e=e._sympy_().simplify()
show(" lambda List : ",lambdaList)
A_lambdaI_list=[]
for l in lambdaList :
A_lambdaI_list.append(A-l*IM)
"An interval is represented as a pair of floating-point numbers a and b (where a≤b)\
and is printed as a standard floating-point number with a question mark (for instance, 3.1416?).\
The question mark indicates that the preceding digit may have an error of ±1. \
These floating-point numbers are implemented using MPFR (the same as the RealNumber elements of RealField_class)."
#show("matrices List of (A-lambda*I) : ")
#for m in A_lambdaI_list :
# m.change_ring(CC)
# show("m : ",m," m.det : ",m.det())
# P,L,U=(m.change_ring(QQbar)).LU()
# show("m : ",m," m.det : ",m.det())
# show("P : ",P)
# show("L : ",L)
# show("U : ",U)
#show(m.numerical_approx(digits=5))
#show(S)
lambdaDiag = diagonal_matrix(lambdaList)
#show(" matrix (diagonal_matrix)",Lambdad.numerical_approx(digits=5))
show("lambda diagonal matrix",lambdaDiag)
##########################################################################
print "solution variables list (replace solution (sols) internal free variables r_i by arbitrary value 1 or 0) \
here this function gives vectors solution for the nullspace matrix (A-lambda_i*I)*X=0 \
in other words, this function gives a solution to the equation (A-lambda_i*I)*X=0 other than the trivial solution \
X=zero vector, this solution is the eigen vector corresponding to the eigen value lambda_i"
def solutionMatricesList(sols,X):
XcList=[]
varList=[]
varStr=""
for eqs in sols :
for eq in eqs :
#show(" eq : ",eq," eq.rhs().variables() : ",eq.rhs().variables())
for v in (eq.rhs().variables()) :
if v not in varList :
varList.append(v)
show(" varList : ",varList," len(varList) : ",len(varList))
if len(varList)==0 :
return XcList,varList
varListValue=[0] * len(varList)
Xc=copy(X)
for i in range(0,len(varList)) :
Xc=(Xc.subs({varList[i]:0 for i in range(0,len(varList))}))
XcList.append(Xc)
varListValue[0] =1
for j in range(0,len(varList)) :
Xc=copy(X)
#show(" varListValue : ",varListValue)
for i in range(0,len(varList)) :
Xc=(Xc.subs({varList[i]:varListValue[i] for i in range(0,len(varList))}))
varListValue.insert(0,varListValue.pop())
XcList.append(Xc)
return XcList,varList
def SolveNUorUNeqM(N,M,nu) :
# solve matrix N*U=M, with U Unknowns matrix
nR=N.nrows()
nC=N.ncols()
mR=M.nrows()
mC=M.ncols()
# generate unknown matrix variable list
#################################
uR=N.ncols()
uC=M.ncols()
uL=[]
for r in [0..uR-1] :
for c in [0..uC-1] :
uL.append("u_%d%d"%(r,c))
# workAound bug sagemath 8.8
# https://ask.sagemath.org/question/47470/var-list-with-only-one-variable-error/
# uV=var(uL) does not work if len(uL)=1
uV=[var(uL[0])] if len(uL) == 1 else var(uL)
#show("uV : ",uV)
#show(" U row : ",uR, " U col : ",uC)
Ul=[[0 for c in [0..uC-1]] for r in [0..uR-1]]
index=0
for r in [0..uR-1] :
for c in [0..uC-1] :
Ul[r][c]=uV[index]
#print "r : " ,r," c : " ,c ," index : ", index ," uV[index] : ", uV[index]
index+=1
U=(matrix(Ul))
#show("U : ",U)
#show("N : ",N)
#show("M : ",M)
#
if nu :
NU=N*U
else :
NU=U*N
eqT=[]
for r in [0..mR-1] :
for c in [0..mC-1] :
eqT.append(NU[r][c]==M[r][c])
# adding suplementary constraints if needed
#eqT.append(u_30==0)
########## end adding suplementary constraints
#for eq in eqT :
# show(eq)
S=solve(eqT,uV)
if len(S)<>0 :
show("Solutions : ",S)
if len(S[0])==U.nrows() * U.ncols() :
#for s in S[0]:
# show(s)
UnumL=[]
# fill the matrix Unum with its numerical values
UnumL=[[0 for c in [0..uC-1]] for r in [0..uR-1]]
index=0
for r in [0..uR-1] :
for c in [0..uC-1] :
UnumL[r][c]=S[0][index].rhs()
index+=1
Unum=matrix(UnumL)
# verify
#show("U : ",U," Unum : ",Unum)
#show("N : ",N)
#show("M : ",M)
#show("N*U :",N*U)
return U,Unum,S
else :
print "number of solutions do not match number of unknowns!"
return U,U
else :
print "No solutions !"
return U,U
zeroL=[0]*mR
Z=matrix(QQ,zeroL).transpose()
M=Z
eigenVectorList=[]
#show(M)
nu=True
M=Z
print "remember that A_lambdaI_list is the list of (A-lambda_i*I) matrices for all the lambda_i"
print "now we will find the non-trivial (not zero vectors) solution for (A-lambda_i*I)*X=0"
index=1
for e in A_lambdaI_list:
# solve matrix N*U=M
U,Unum,S=SolveNUorUNeqM(e,M,nu)
#show(" N : ",N," U : ", U, " Unum : ",Unum ," = M : ",M)
#if nu :
# show(" N*Unum : ",N*Unum)
#else :
# show(" Unum*U : ",Unum*N)
XcList,varList=solutionMatricesList(S,Unum)
#show(" XcList : ", XcList)
if varList <>0 :
for ee in XcList :
if ee <> Z :# remove the trivial one
eigenVectorList.append(ee)
show("A-lambda_i*I matrix"+str(index)+" :",e," eigen vector x_" + str(index)+" : ",ee)
index+=1
#show("verify the result is zero vector : ",e*ee)
#show("eigenVectorList : ")
#show(eigenVectorList)
# concatenate the eigen vectors in a matrix
X=eigenVectorList[0]
for i in range(1,len(eigenVectorList)) :
X=X.augment(eigenVectorList[i])
show("eigen vectors matrix X : ",X)
Xinv=X.inverse().simplify_full()
show("X inverse : ",Xinv)
show(LatexExpr(r" \text{the greek letter } \, \Lambda \,\, \text{represents the matrix which is called lambdaDiag in the code (matrix of eigen values in diagonal)} \\ \
\text{ and matrix of eigen vectors X is also called X in the code } \\ \
\text{diagonalization: We know that} \, X^{-1} \, A \, X = \Lambda \, \text{so} \, A= X \, \Lambda \, X^{-1} "))
show(LatexExpr(r"\text{Eigenvalues of } \,A^k = \, (X \, \Lambda \, X^{-1})^k \,= \, (X \, \Lambda \, X^{-1}) \, (X \, \Lambda \, X^{-1}) \, \text{...} _, (X \, \Lambda \, X^{-1}) \
= \, (X \, \Lambda \, (X^{-1} \, X) \, \Lambda \, \text{...} \, ( X^{-1} \, X) \, \Lambda \, X) = (X \, \Lambda^k \, X^{-1})"))
show("Verify that A=X*lambdaDiag*Xinv : ", (X*lambdaDiag*Xinv).simplify_full())
print "here below using power of the matrix lambdaDiag^k, which is not optimal, "
print "X*(lambdaDiag^4)*Xinv*startSeqM compute the k=4 next values following the starter Sequence startSeqM:"
show((X*(lambdaDiag^4)*Xinv*startSeqM).simplify_full())
k=var('k')
lambdaPowKlist=[]
for i in range(0,len(lambdaList)) :
lambdaPowKlist.append(lambdaList[i]^k)
show("lambda^k list : ",lambdaPowKlist)
lambdaDiagPowK = diagonal_matrix(lambdaPowKlist)
show("lambda^t Diagonal matrix : ",lambdaDiagPowK)
ApowK=(X*(lambdaDiagPowK)*Xinv*startSeqM).simplify_full()
#show(ApowT.dimensions(),ApowT)
print"below using the matrix lambdaDiag^k with all the diagonal eingenvalues to the exponent k,\
(which is the interest of this procedure with a chalck)"
show("X*(lambdaDiag^4)*Xinv*startSeqM= ",(X*(lambdaDiagPowK.subs(k=4))*Xinv*startSeqM).simplify_full())
show("below we compute (for example here, we choose case p=3 p+1=4)" ,LatexExpr(r" \lim_{n\to\infty} \frac{a_{n+1}}{a_n} ")," using the matrices product : ",\
LatexExpr(r"\, X \Lambda^k X^{-1} \,"), " which gives the resultant vector",\
LatexExpr(r"\, \left(\begin{array}{r}\
a_{n+4+k} \\ \
a_{n+4+k-1} \\ \
a_{n+4+k-2} \\ \
a_{n+4+k-3} \
\end{array}\right)"))
show(LatexExpr(r" X= \text{eigen vector Matrix and } \,\Lambda \,\text{ = } \lambda \text{ diagonal matrix}"))
show("and now we are looking for ",LatexExpr(r"\lim_{n\to\infty} \, \frac{a_{n+4+k}}{a_{n+4+k-1}}")," which is ",
LatexExpr(r" \lim_{n\to\infty} \frac{a_{n+1}}{a_n} ") )
print "in case of p=2,we jump this Limit computation as it produces a PC's memory overflow, \
but the limit is well computed for p=3"
show(LatexExpr(r" \text{for p=2 the limit would be equals to } \phi_2 =\frac{1}{3} \, \left(\frac{1}{2}\right)^{\frac{1}{3}} {\left(3 \, \sqrt{23} \sqrt{3} + \
25\right)}^{\frac{1}{3}} + \frac{2 \, \left(\frac{1}{2}\right)^{\frac{2}{3}}}{3 \, {\left(3 \, \sqrt{23} \sqrt{3} + \
25\right)}^{\frac{1}{3}}} + \frac{2}{3} \approx 1.754877666246693 "))
show(LatexExpr(r"\phi_2 \,\text{is the real root of } \, p(x)=x^3 -2x^2 + x -1"))
print "it seems the Algebra Limit is too much complicated for p=2"
print "see page 205 of THE BOOK OF NUMBERS, John H. Conway • Richard K. Guy "
print " http://www.blackwire.com/~bjordan/The-Book-of-Numbers.pdf "
if p <> 2 :
print "but here as p="+str(p)+"we can find the limit without memory overflow"
ApowK=(X*(lambdaDiagPowK)*Xinv*startSeqM)
anPlusOne_an=(sum(ApowK[0,:][0]).real_part()/sum(ApowK[1,:][0]).real_part())
L=limit(anPlusOne_an,k=infinity)
show(LatexExpr(r" \lim_{n\to\infty} \frac{a_{n+1}}{a_n} = "), L._sympy_().simplify())