Ask Your Question
0

impossible to ask my question [closed]

asked 2019-08-31 06:02:25 -0500

ortollj gravatar image

updated 2019-08-31 13:24:49 -0500

Hi

I got this error:Internal server error system error log is recorded, error will be fixed as soon as possible please report the error to the site administrators I think that the code size is too large ? (clik on more comments below the last displayed comment on this page) if you want to see my question click the fourth comment google doc link. the code is ok for p=3 (I got the good limit) and ok for p=2 except the computation of the Limit at the end of the code

it is worse now I can't even post a comment on this forum !!!

Please help me: do you think all this mess comes from my PC ?. I can only write in the first message and in the response.

edit retag flag offensive reopen merge delete

Closed for the following reason not a real question by Iguananaut
close date 2019-09-04 09:23:56.830437

Comments

I cant also (without the code) put it inside a comment ! maybe I m infected by a virus ?

ortollj gravatar imageortollj ( 2019-08-31 06:08:53 -0500 )edit

strange my web browser Edge became blank when I try to ask my question , I'm gonna try with another web browser, Chrome

ortollj gravatar imageortollj ( 2019-08-31 06:17:31 -0500 )edit

no it is impossible also with Chrome ! Crazy !!!

ortollj gravatar imageortollj ( 2019-08-31 06:20:27 -0500 )edit

I put my question in google doc: My question

ortollj gravatar imageortollj ( 2019-08-31 06:28:11 -0500 )edit
1

Could you try without the link ?

tmonteil gravatar imagetmonteil ( 2019-08-31 06:33:23 -0500 )edit

1 answer

Sort by » oldest newest most voted
0

answered 2019-08-31 12:27:32 -0500

ortollj gravatar image
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())
edit flag offensive delete link more

Comments

I finally got to put the code, for that I put it in the location for the answer !! which is not opportune !

ortollj gravatar imageortollj ( 2019-08-31 12:31:00 -0500 )edit

I m gonna erase this messy question, and try to ask my original question more clearly. because it seems that my comment strange problem has disappeared now.

ortollj gravatar imageortollj ( 2019-09-01 02:35:24 -0500 )edit

no it is always not possible for me to ask my question every attempt I get:

Internal server error system error log is recorded, error will be fixed as soon as possible please report the error to the site administrators

so for the moment I do not erase this messy question ;-(

ortollj gravatar imageortollj ( 2019-09-01 03:26:41 -0500 )edit

Question Tools

1 follower

Stats

Asked: 2019-08-31 06:02:25 -0500

Seen: 92 times

Last updated: Aug 31 '19